Maple hints for doing Homework 12 (week 11)

(The divergence theorem)

section 10.5 of Prof. Malek-Madani's book Advanced Engineering Mathematics

sm311o_hwk12a.mws,wdj,4-17-98

> restart;with(plots):
with(linalg):

Warning, new definition for norm

Warning, new definition for trace

Problems #1, 2 : These problems can be done in the same manner as the following one:

Verify the divergence theorem for F(x,y,z)=(y,x,z) for the tetrahedron with vertices (0,0,0), (1,0,0),

(0,2,0),(0,0,3).

This means "compute both sides of the divergence theorem

surface integral of F over S = triple integral of div(F) over Q

and check that they are equal".

> F:=(x,y,z)->[y,x,z];

The easiest way to do this in MAPLE is to use the formula

surface integral of F over S = int int_S F(x,y,z)*n d(sigma)

The tetrahedron S is composed over 4 faces: S1 = front (x+y/2+z/3=1), S2 = left side (y=0), S3 = back (x=0), S4 = bottom (z=0). To get a unit normal, take the gradients of these equations and divide by the length. To make them be outward unit normals, we sometimes multiply by -1. The formula for d(sigma) involves the area element factor sqrt(f_x^2+f_y^2+f_z^2)/|f_?|, where f(x,y,z)=c describes the surface and ? is the variable (x,y, or z) along whose axis you projected the surface (eg, ?=z if you projected the surface down onto the xy plane to integrate). We must compute this as well.

> v1:=convert(grad(x+y/2+z/3,[x,y,z]),list);
n1:=(x,y,z)->v1/sqrt(dotprod(v1,v1)):
n1(x,y,z);
area_element_factor1:=sqrt(v1[1]^2+v1[2]^2+v1[3]^2)/v1[3];

> v2:=convert(grad(y,[x,y,z]),list);
n2:=(x,y,z)->-v2/sqrt(dotprod(v2,v2)):
n2(x,y,z);
area_element_factor2:=sqrt(v2[1]^2+v2[2]^2+v2[3]^2)/v2[2];

> v3:=convert(grad(x,[x,y,z]),list);
n3:=(x,y,z)->-v3/sqrt(dotprod(v3,v3)):
n3(x,y,z);
area_element_factor3:=sqrt(v3[1]^2+v3[2]^2+v3[3]^2)/v3[1];

> v4:=convert(grad(z,[x,y,z]),list);
n4:=(x,y,z)->-v4/sqrt(dotprod(v4,v4)):
n4(x,y,z);
area_element_factor4:=sqrt(v4[1]^2+v4[2]^2+v4[3]^2)/v4[3];

The base for the integrals over S1 and S4 are both {0 < x < 1, 0 < y < 2(1 - x) }. On S1, z = 3(1 - x - y/2). We can check this visually in MAPLE by simply plotting the 3d object using these ranges:

NOTE: MAPLE does not plot the xyz axes in the usual right-hand orientation we use but plots the yxz axes instead .

> P1:=plot3d(3*(1-x-y/2),y=0..(2*(1-x)),x=0..1,axes=normal,style=wireframe,color=grey):
P1;

> S1_int:=Int(Int(Dotprod(F(x,y,3*(1-x-y/2)),n1(x,y,3*(1-x-y/2)))*area_element_factor1,y=0..(2*(1-x))),x=0..1);
int(int(dotprod(F(x,y,3*(1-x-y/2)),n1(x,y,3*(1-x-y/2)))*area_element_factor1,y=0..(2*(1-x))),x=0..1);

On S4, z=0 but the range of integration is the same

> S4_int:=Int(Int(Dotprod(F(x,y,0),n4(x,y,0))*area_element_factor4,y=0..(2*(1-x))),x=0..1);
int(int(dotprod(F(x,y,0),n4(x,y,0))*area_element_factor4,y=0..(2*(1-x))),x=0..1);

The range for the integral over S2 is {0 < x < 1, 0 < z < 3(1 - x) }. We can check this visually in MAPLE by simply plotting the 3d object using these ranges:

> P2:=plot3d(0,z=0..(3*(1-x)),x=0..1,axes=normal,style=wireframe,color=grey):
P2;

> S2_int:=Int(Int(Dotprod(F(x,0,z),n2(x,0,z))*area_element_factor2,z=0..(3*(1-x))),x=0..1);
int(int(dotprod(F(x,0,z),n2(x,0,z))*area_element_factor2,z=0..(3*(1-x))),x=0..1);

The range for the integral over S3 is {0 < y < 2, 0 < z < 3(1 - y/2) }. We can check this visually in MAPLE by simply plotting the 3d object using these ranges

> P3:=plot3d(0,z=0..(3*(1-y/2)),y=0..2,axes=normal,style=wireframe,color=grey):
P3;

> S3_int:=Int(Int(Dotprod(F(0,y,z),n3(0,y,z))*area_element_factor3,z=0..(3*(1-y/2))),y=0..2);
int(int(dotprod(F(0,y,z),n3(0,y,z))*area_element_factor3,z=0..(3*(1-y/2))),y=0..2);

>

> total:=S1_int+S2_int+S3_int+S4_int;
evalf(total);

Now we must compute the triple integral of the divergence of F

> triple_int:=Int(Int(Int(div(F(x,y,z),[x,y,z]),z=0..(3(1-x-y/2))),y=0..(2*(1-x))),x=0..1);
int(int(int(diverge(F(x,y,z),[x,y,z]),z=0..(3*(1-x-y/2))),y=0..(2*(1-x))),x=0..1);

Problems 5,6 :

Both require you to do an integral over a sphere or hemisphere centered about some point (x0,y0,z0). Here's the parameterization, where a is the radius and (x0,y0,z0) is the center:

> restart;with(plots):
with(linalg):

Warning, new definition for norm

Warning, new definition for trace

> r:=(u,v,a,x0,y0,z0)->[a*cos(u)*sin(v)+x0,a*sin(u)*sin(v)+y0,a*cos(v)+z0]:
r(theta,phi,a,0,0,0);

> F:=(x,y,z)->[x^3,y^2-z^2,x^2-y^2]:
F(x,y,z);

> ru:=(u0,v0)->subs({u=u0,v=v0},diff(r(u,v,1/n,1,1,1),u)):
ru(u,v);
rv:=(u0,v0)->subs({u=u0,v=v0},diff(r(u,v,1/n,1,1,1),v)):
rv(u,v);
n0:=(u,v)->crossprod(ru(u,v),rv(u,v)):
n0(u,v);

> vol:=n->4*Pi*(1/n)^3:
a:=n->subs(m=n,int(int(dotprod(F(op(r(u,v,1/m,1,1,1))),n0(u,v),orthogonal),u=0..(2*Pi)),v=0..Pi)/vol(m)):
a(n);

> t0:=time();

> seq(evalf(a(n)),n=1..50);
print(time()-t0, `seconds used to compute this sequence`);

>