# MAT 223 Multiple Integrals

## 1585 days ago by kcrisman

What IS an antiderivative?

var('y') @interact def _(t = [pi/6*n for n in [2..18]],start=[0,pi/6]): P = plot(1/2*sin(y),(x,start,t),fill=True,tick_formatter=pi,ticks=pi/6) Q = plot(1/2*sin(y),(x,start,3*pi),tick_formatter=pi,ticks=pi/6) show(P+Q,figsize=5) pretty_print(html("The integral up to $y=%s$ is $\\approx %s$."%(latex(t),latex(integral(1/2*sin(y),y,start,t).n()))))

## Click to the left again to hide and once more to show the dynamic interactive window

Here's the multivariate version, then.

var("y z") @interact #def _(s = 1/2,t=(3*pi/6,[pi/6*n for n in [1..18]]),show_x_plane=False,show_z_plane=False,show_all=False): def _(s = 1/2,show_x_plane=False,show_z_plane=False,show_all=False): P = plot3d(x*sin(y),(x,0,2),(y,0,3*pi)) if show_x_plane: u,v = var('u v') P += parametric_plot3d([s,v,s*u*sin(v)],(u,0,1),(v,0,3*pi),color='red',alpha=.5) if show_z_plane: P += plot3d(0,(x,0,2),(y,0,3*pi),color='orange',alpha=.5) #if show_fill: # P += parametric_plot3d((s,y,x*sin(y)*s),(y,0,t),(x,min(0,x*sin(y)*s),max(0,x*sin(y)*s)),color='green') P += parametric_plot3d((s,y,s*sin(y)),(y,0,3*pi),color='green',thickness=10) if show_all: P += sum([parametric_plot3d((u,y,u*sin(y)),(y,0,3*pi),color='green',thickness=10) for u in xsrange(0,s,.05)]) show(P)

## Click to the left again to hide and once more to show the dynamic interactive window

s=1 P = plot3d(x*sin(y),(x,0,1),(y,0,pi/2)) P += implicit_plot3d(x==s,(x,0,1),(y,0,pi/2),(z,-2,2),color='red',alpha=.5) P += plot3d(0,(x,0,1),(y,0,pi/2),color='orange',alpha=.5) P += parametric_plot3d((s,y,s*sin(y)),(y,0,pi/2),color='green',thickness=10) show(P)
var('y') integral(integral(x*sin(y),y,0,pi/2),x,0,1)
 1/2 1/2
interacts.calculus.riemann_sum()

## Riemann integral with random sampling

$f(x)=$
# divisions

Integration interval
slider:
keyboard:

List table

## Riemann integral with random sampling

$f(x)=$
# divisions

Integration interval
slider:
keyboard:

List table

The next example is by Marshall Hampton.

from sage.plot.plot3d.platonic import index_face_set def cuboid(v1,v2,**kwds): """ Cuboid defined by corner points v1 and v2. """ ptlist = [] for vi in (v1,v2): for vj in (v1,v2): for vk in (v1,v2): ptlist.append([vi[0],vj[1],vk[2]]) f_incs = [[0, 2, 6, 4], [0, 1, 3, 2], [0, 1, 5, 4], [1, 3, 7, 5], [2, 3, 7, 6], [4, 5, 7, 6]] if 'aspect_ratio' not in kwds: kwds['aspect_ratio'] = [1,1,1] return index_face_set(f_incs,ptlist,enclosed = True, **kwds) var('x,y') R16 = RealField(16) npi = RDF(pi) sin,cos = math.sin,math.cos pretty_print(html("<h1>The midpoint rule for a function of two variables</h1>")) @interact def midpoint2d(func = input_box('y*sin(x)/x+sin(y)',type=str,label='function of x and y'), nx = slider(2,20,1,3,label='x subdivisions'), ny = slider(2,20,1,3,label='y subdivisions'), x_start = slider(-10,10,.1,0), x_end = slider(-10,10,.1,3*npi), y_start= slider(-10,10,.1,0), y_end= slider(-10,10,.1,3*npi)): f = sage_eval('lambda x,y: ' + func) delx = (x_end - x_start)/nx dely = (y_end - y_start)/ny xvals = [RDF(x_start + (i+1.0/2)*delx) for i in range(nx)] yvals = [RDF(y_start + (i+1.0/2)*dely) for i in range(ny)] num_approx = 0 cubs = [] darea = delx*dely for xv in xvals: for yv in yvals: num_approx += f(xv,yv)*darea cubs.append(cuboid([xv-delx/2,yv-dely/2,0],[xv+delx/2,yv+dely/2,f(xv,yv)], opacity = .5, rgbcolor = (1,0,0))) pretty_print(html("$$\int_{"+str(R16(y_start))+"}^{"+str(R16(y_end))+"} "+ "\int_{"+str(R16(x_start))+"}^{"+str(R16(x_end))+"} "+func+"\ dx \ dy$$")) pretty_print(html('<p style="text-align: center;">Numerical approximation: ' + str(num_approx)+'</p>')) p1 = plot3d(f,(x,x_start,x_end),(y,y_start,y_end)) show(p1+sum(cubs))

## Click to the left again to hide and once more to show the dynamic interactive window

(pi/4/sqrt(2)).n()
 0.555360367269796 0.555360367269796
s=1 P = plot3d(1,(x,0,1),(y,0,1)) P += implicit_plot3d(x==s,(x,0,1),(y,0,1),(z,-2,2),color='red',alpha=.5) P += plot3d(0,(x,0,1),(y,0,pi/2),color='orange',alpha=.5) P += parametric_plot3d((s,y,s*sin(y)),(y,0,pi/2),color='green',thickness=10) show(P)

This one is by John Travis at Mississippi College.

var('x,y') # f is the top surface # g is the bottom surface global f,g # condition1 and condition2 are the inequality constraints. It would be nice # to have any number of conditions connected by  or | global condition1,condition2 @interact def _(f=input_box(default=1+0*x,label='$f(x)=$'), g=input_box(default=0+0*x,label='$g(x)=$'), condition1=input_box(default= x^3<y,label='$Constraint_1=$'), condition2=input_box(default=y<x^2,label='$Constraint_2=$'), x_slice = .5, show_slices=false, show_3d=('Stereographic',false), show_vol=('Shade volume',false), dospin = ('Spin?',false), clr = color_selector('#faff00', label='Volume Color', widget='colorpicker', hide_box=True), xx = range_slider(-5, 5, 1, default=(0,1), label='X Range'), yy = range_slider(-5, 5, 1, default=(0,1), label='Y Range'), auto_update=false): # This is the top function actually graphed by using NaN outside domain def F(x,y): if condition1(x=x,y=y): if condition2(x=x,y=y): return f(x=x,y=y) else: return -NaN else: return -NaN # This is the bottom function actually graphed by using NaN outside domain def G(x,y): if condition1(x=x,y=y): if condition2(x=x,y=y): return g(x=x,y=y) else: return -NaN else: return -NaN P = Graphics() # The graph of the top and bottom surfaces P_list = [] P_list.append(plot3d(F,(x,xx[0],xx[1]),(y,yy[0],yy[1]),color='blue',opacity=0.9,plot_points=100)) P_list.append(plot3d(G,(x,xx[0],xx[1]),(y,yy[0],yy[1]),color='gray',opacity=0.9,plot_points=100)) # Interpolate "layers" between the top and bottom if desired if show_vol: ratios = range(10) def H(x,y,r): return (1-r)*F(x=x,y=y)+r*G(x=x,y=y) P_list.extend([ plot3d(lambda x,y: H(x,y,ratios[1]/10),(x,xx[0],xx[1]),(y,yy[0],yy[1]),opacity=0.2,color=clr), plot3d(lambda x,y: H(x,y,ratios[2]/10),(x,xx[0],xx[1]),(y,yy[0],yy[1]),opacity=0.2,color=clr), plot3d(lambda x,y: H(x,y,ratios[3]/10),(x,xx[0],xx[1]),(y,yy[0],yy[1]),opacity=0.2,color=clr), plot3d(lambda x,y: H(x,y,ratios[4]/10),(x,xx[0],xx[1]),(y,yy[0],yy[1]),opacity=0.2,color=clr), plot3d(lambda x,y: H(x,y,ratios[5]/10),(x,xx[0],xx[1]),(y,yy[0],yy[1]),opacity=0.2,color=clr), plot3d(lambda x,y: H(x,y,ratios[6]/10),(x,xx[0],xx[1]),(y,yy[0],yy[1]),opacity=0.2,color=clr), plot3d(lambda x,y: H(x,y,ratios[7]/10),(x,xx[0],xx[1]),(y,yy[0],yy[1]),opacity=0.2,color=clr), plot3d(lambda x,y: H(x,y,ratios[8]/10),(x,xx[0],xx[1]),(y,yy[0],yy[1]),opacity=0.2,color=clr), plot3d(lambda x,y: H(x,y,ratios[9]/10),(x,xx[0],xx[1]),(y,yy[0],yy[1]),opacity=0.2,color=clr) ]) # P = plot3d(lambda x,y: H(x,y,ratio/10),(x,xx[0],xx[1]),(y,yy[0],yy[1]),opacity=0.1) P_list.append(parametric_plot3d((x_slice,y,x),(y,x_slice^3,x_slice^2),(x,g(x=x_slice),f(x=x_slice)),color='red')) if show_slices: for i in [.1,.2,.3,.4,.5,.6,.7,.8,.9]: x_slice=i P_list.append(parametric_plot3d((x_slice,y,x),(y,x_slice^3,x_slice^2),(x,g(x=x_slice),f(x=x_slice)),color='red')) # Now, accumulate all of the graphs into one grouped graph. P = sum(P_list[i] for i in range(len(P_list))) if show_3d: show(P,frame=true,axes=false,xmin=xx[0],xmax=xx[1],ymin=yy[0],ymax=yy[1],stereo='redcyan',figsize=(6,9),viewer='jmol',spin=dospin) else: show(P,frame=true,axes=false,xmin=xx[0],xmax=xx[1],ymin=yy[0],ymax=yy[1],figsize=(6,9),viewer='jmol',spin=dospin)

## Click to the left again to hide and once more to show the dynamic interactive window

P = plot(x,(x,0,1),fill=True)
P.plot3d()

Finally, to get multiple integrals, you'll want to do two integrals at once.

var('y') integral(integral(x-x^3-y^2*x,y,0,1-x),x,0,1)
 1/10 1/10

Or perhaps to numerically approximate them...

var('y') assume(x<1) # needed for the integration algorithm to work integral(integral(sqrt(4*x^2+4*y^2+1),y,0,1-x),x,0,1).n() # needed because no elementary antiderivative!
 0.7511563585703602 0.7511563585703602
var('y') integral(integral(sqrt(36*x^2+17),y,0,2*x),x,0,2).n()
 36.5327442292013 36.5327442292013