L10-SurfInt.mws

Calculus IV with Maple
Copyright 2002, Dr. Jack Wagner
j.wagner@intelligentsearch.com

Lesson 10: Surface Integrals

Topics: Integration over a surface.  Integrating over the Moebius band.

Maple commands introduced:   stackmatrix, augment

In the same way that a vector field may be integrated over a curve, it may be integrated over a surface.  To each parallelogram that forms an element of the area of the surface we assign the normal component of the vector field at some interior point.  As the partition of the surface is refined  the sum of the products of the area of the parallelograms and the normal component of the vector field is the integral of the vector field over the surface, usually written, Int(F,S)   , where S is the surface, or Int(F[U].n,A) , where dA is understood to represent the "element of area", n is the unit normal and U is the  parameterization.  This integral is the flux of F through the surface S. Because the normal to the surface plays a key role it is important that we work only with surfaces over which it is possible to consistently define a unit normal.  That is, there must be a vector valued function, n (p), that assigns a unique unit normal to each point, p. Such a surface is orientable.

The archetypical non-orientable surface is the Moebius band, which may be parameterized and plotted as follows.

>    restart: with(plots):with(LinearAlgebra): with(VectorCalculus):

>    F := <(2-v*sin(u/2))*sin(u), (2-v*sin(u/2))*cos(u), v*cos(u/2)>;

F := Vector(%id = 20937060)

>    plot3d(F, u=-Pi..Pi, v=-1..1, axes=framed, style=wireframe);

[Maple Plot]

It is not orientable because at each point the "function"  n(v, u),  creates two unit normals  each pointing in the opposite direction. For a surface S(u,v) = [x(u,v), y(u,v), z(u,v)]   define n(u,v) = S[u]*X*S[v]/abs(S[u]*X*S[v])  .

If we plot n (u, v) over a surface, we expect to get a copy of the original surface one unit away.  For the Moebius band, however:

>    dF[u] := diff(F, u);

dF[u] := Vector(%id = 20940140)

>    dF[v] := diff(F, v);

dF[v] := Vector(%id = 20940380)

>    n := dF[v] &x dF[u]:

>    N := Normalize(n, 2, inplace):

>    P1 := plot3d(F+N, u=-2*Pi..2*Pi, v=
-1..1, axes=framed, style=wireframe, color=blue, grid=[40,40]):

>    P2 := plot3d(F, u=-Pi..Pi, v=-1..1, axes=framed, color=pink):

>    display({P1, P2}, scaling=constrained, orientation=[-135,75]);

[Maple Plot]

At each point on the surface u(v, u) has created two unit normals and therefore it is not well defined on this surface.  A bit further on we will look at how non-orientability interferes with integration of a vector field over a surface.  We are also going to confine ourselves to simple surfaces, that is, surfaces which do not contain singularities such as self intersections.

Let F be the vector field. The unit normal to the surface S, known by S(u,v) = [x(u,v), y(u,v), z(u,v)] , is  

S[u]*X*S[v]/abs(S[u]*X*S[v])  . We already know that dA = abs(S[u]*X*S[v])   dudv  .  Therefore, Int(F*n,A)  = Int(Int(F*S[u]*X*S[v]*abs(S[u]*X*S[v])/abs(S[u]*X*S[v]),u),v)  =

Int(Int(F*S[u]*X*S[v],u),v)  

Example 10.1

Find the flux of F(x,y,z) = [1/(x^2), 1/(y^2), 1/(z^2)]   through the surface of the paraboloid z = x^2+y^2    1 <= x , y <= 3

>    restart: with(plots): with(LinearAlgebra): with(VectorCalculus):

>    F := <1/x^2, 1/y^2, 1/z^2>;

F := Vector(%id = 18638708)

Now select a parametrizarion for S. {x=u, y=v}

>    S := <u, v, u^2+v^2>;

S := Vector(%id = 18314476)

>    F := eval(F, {x=S[1], y=S[2], z=S[3]});

F := Vector(%id = 17955144)

>    P1 := plot3d(S, u=1..3, v=1..3, color=aquamarine, style=wireframe, axes=framed):

>    P2 := fieldplot3d(F, u=1..3, v=1..3, z=1..18, arrows=THICK, grid=[5, 5, 5], color=black):

>    display(P1, P2);

[Maple Plot]

Now we compute.

>    dS[u] := diff(S, u);

dS[u] := Vector(%id = 2690920)

>    dS[v] := diff(S, v);

dS[v] := Vector(%id = 17537628)

>    n := dS[v] &x dS[u];

n := Vector(%id = 2645640)

Note that we have selected the "outward" pointing normal.

>    FdA := F.n;

FdA := 2/u+2/v-1/((u^2+v^2)^2)

The flux over the surface is

>    int(int(FdA, u=1..3), v=1..3);

-1/18*arctan(1/3)-1/2*arctan(3)+1/9+8*ln(3)+5/36*Pi

Or numerically,

>    evalf(%);

8.693943819

In the standard basis, { i, j, k }, S[u]*X*S[v] = Det(matrix([[i, j, k], [x[u], y[u], z[u]], [x[v], y[v], z[v]]]))  where the subscripts indicate the derivative.  This can be evaluated as   i Det(matrix([[y[u], z[u]], [y[v], z[v]]]))  - j Det(matrix([[x[u], z[u]], [x[v], z[v]]]))  + k   Det(matrix([[x[u], y[u]], [x[v], x[v]]]))  . Therefore,   F dA =

Det(matrix([[F[1], F[2], F[3]], [x[u], y[u], z[u]], [x[v], y[v], z[v]]]))  which yields Int(F*n,A)  = Int(Int(Det(matrix([[F[1], F[2], F[3]], [x[u], y[u], z[u]], [x[v], y[v], z[v]]])),u),v)  

Example 10.2

Integrate   F = [x^2, y^2, z^2]  over the surface defined by   z = sin(x+y) , 0 <= x , y <= Pi

>    restart: with(plots): with(LinearAlgebra): with(VectorCalculus):with(linalg):

>    F := <x^2, y^2, z^2>;

F := Vector(%id = 20784232)

>    S := <u, v, sin(u+v)>;

S := Vector(%id = 20898604)

>    F := eval(F, {x=S[1], y=S[2], z=S[3]});

F := Vector(%id = 20898424)

>    dS[u] := diff(S, u);

dS[u] := Vector(%id = 20898124)

>    dS[v] := diff(S, v);

dS[v] := Vector(%id = 20867948)

>    FdA := <<F>|<dS[u]>|<dS[v]>>;

FdA := Matrix(%id = 20403880)

>    Int(Int(Determinant(FdA), u=0..Pi), v=0..Pi);

Int(Int(-u^2*cos(u+v)-v^2*cos(u+v)+sin(u+v)^2,u = 0 .. Pi),v = 0 .. Pi)

>    value(%);

9/2*Pi^2-16

The determinant formulation leads to the following interpretation.  Consider the first term in the expansion of the determinant in minors of the first row.    F[1]*Det(matrix([[y[u], z[u]], [x[v], z[v]]]))   -   F[2]*Det(matrix([[x[u], z[u]], [x[v], z[v]]]))    +     F[3]*Det(matrix([[x[u], y[u]], [x[v], x[v]]]))  .    F[1]   is the component of F in the x-y plane,    Det(matrix([[y[u], z[u]], [y[v], z[v]]])) = y[u]*z[v]-z[u]*y[v] . is the area of the projection, onto the x-y plane,  of the element of area on the surface of S.  Similarly for the other two components, except that the x-z plane projection has a negative sign.  (Rotate the basis vector of the x axis into the z axis basis vector and, by the right hand rule, you will find that the resulting vector has an opposite sense to that of the y axis.)  The sum of these projections is precisely the determinant given.

From the basic definition of the surface integral, i.e. Int(F.n,A)  we also obtain
Int(F.n,A) = Int(Int(F.n*abs(S[u]*X*S[v]),u),v)  = Int(Int(F.n*sqrt(E*G-F^2),u),v)  

Example 10.3

Compute the integral of F(x,y) = [x, 2^y, 2^(-y)]  over the surface, S = [x, y, -2*x-3*y]  ,    0 <= x , y <= 1  

>    restart:with(plots):with(LinearAlgebra): with(VectorCalculus):

>    f := <x, 2^y, 2^(-y)>;

f := Vector(%id = 2654888)

>    S := <x, y, -2*x-3*y>;

S := Vector(%id = 18935356)

>    P1 := plot3d(S, x=0..1, y=0..1, color=blue, style=wireframe):

>    P2 := fieldplot3d(f, x=0..1, y=0..1, z=-6..1, color=red, grid=[5, 5, 5], thickness=2):

>    display(P1, P2, axes=normal );

[Maple Plot]

For a plane whose equation is 3x-2y-z=0, we know that  [3, -2, -1] is the principal part of a normal vector at every point.

>    N := <3, 2, 1>;

N := Vector(%id = 18829524)

>    n := Normalize(N, 2, inplace);

n := Vector(%id = 18829524)

>    s[x] := diff(S, x);

s[x] := Vector(%id = 18217820)

>    s[y] := diff(S, y);

s[y] := Vector(%id = 17194284)

>    E := s[x].s[x];

E := 5

>    F := s[x].s[y];

F := 6

>    G := s[y].s[y];

G := 10

>    sqrt(E*G-F^2);

14^(1/2)

>    FdA := f.n*sqrt(E*G-F^2);

FdA := (3/14*x*14^(1/2)+1/7*2^y*14^(1/2)+1/14*2^(-y)*14^(1/2))*14^(1/2)

>    Int(Int(FdA, x=0..1), y=0..1);

Int(Int((3/14*x*14^(1/2)+1/7*2^y*14^(1/2)+1/14*2^(-y)*14^(1/2))*14^(1/2),x = 0 .. 1),y = 0 .. 1)

>    value(%);

1/2*(5+3*ln(2))/ln(2)

>    evalf(%);

5.106737600

Now let's look at what happens when we try to integrate a function over the surface of a Moebius band.

Example 10.4

Calculate the integral of

S = [(2-v*sin(u/2))*sin(u), 2-v*sin(u/2)*cos(u), v*cos(u/2)]

>    restart: with(plots): with(LinearAlgebra): with(VectorCalculus):

>    S := <(2-v*sin(u/2))*sin(u), (2-v*sin(u/2))*cos(u), v*cos(u/2)>;

S := Vector(%id = 19827304)

>    F := <u^2+v, v^2+u, sin(u+v)>;

F := Vector(%id = 19827504)

>    dS[u] := diff(S, u);

dS[u] := Vector(%id = 17827228)

>    dS[v] := diff(S, v);

dS[v] := Vector(%id = 19973280)

>    FdA := <<dS[u]>|<dS[v]>|<F>>;

FdA := Matrix(%id = 19991084)

Now we integrate over the surface of the Moebius band.  First in a positive and then in a negative direction.  The results should sum to zero.  They do not!

>    int(int(Determinant(FdA), u=0..2*Pi), v=-1..1);

2272/27+1/3*Pi-Pi*sin(1)+Pi*cos(1)-64/3*Pi^2

>    int(int(Determinant(FdA), u=0..-2*Pi), v=-1..1);

2272/27-64/3*Pi^2-1/3*Pi+Pi*sin(1)-Pi*cos(1)

>    %+%%;

4544/27-128/3*Pi^2

The surface integral of F computes the volume of a solid whose base is the surface and whose height is the normal projection of F.   If F happens to be the unit normal, then   Int(n*n,A) = Int(1,A)  = Area.

Practice

1.       Compute the indicated surface integrals.
           a)  
F(x,y) = [x, sin(y), cos(y)]  , S(x,y) = [x, y*cos(y), y*sin(y)]  , x=0..1,  y=0.. 2*Pi

           b)   F(x,y) = [x, x, exp(-.1*y)]  , S(x,y) = [y*cos(x), y*sin(x), y]  ,  x=-6..6, y=0.. 2*Pi

           c)   F(x,y) = [exp(.1e-1*x), exp(.1e-1*y), exp(.1e-1*(x+y))]  , S(x, y) = {x, y, z| x+y+3*z = 1  },  { x, y }| x^2+y^2 = 4    

2.       Using Maple, show that for a vector field, F, symmetrical around the axis of revolution, the surface integral over a surface of revolution, between t[1]   and t[2] ,  is F[1]*(A[t[1]]-A[t[2]])  , where F[1]   is the vector field component along the axis of revolution and   A[j]  is the area of the surface disc at t[j]  .