# Checking The Divergence Theorem

## Problem

Compute the flux of the vector field

> F := x^2*i+x*y*j+x^3*y^3*k;

```                                2              3  3
F := x  i + x y j + x  y  k```
 over the surface of the tetrahedron bounded by the plane x + y + z = 1 and the coordinate planes. Compute this directly as the surface integral of the projection of the vector field on the outward unit normal and as a volume integral using the divergence theorem.

### Solution:

 Recall the expression of the Divergence Theorem:

> ;

```                   /   /              /   /   /
|   |              |   |   |
|   |  F.N dS   =  |   |   |  div F dV
|   |              |   |   |
/   /              /   /   /
S                    V```
 Let's start with the right hand side tripple integral. First we compute div F.

> divF := diverge(evalm(F),[x,y,z]);

`                                  divF := 3 x`
 and we then integrate over the solid tetrahedron,

> I2 := Int(Int(Int(divF,z=0..1-x-y),y=0..1-x),x=0..1);

```                          1    1 - x    1 - x - y
/    /        /
|    |        |
I2 :=  |    |        |           3 x dz dy dx
|    |        |
/    /        /
0    0        0```
> I2 := int(int(int(divF,z=0..1-x-y),y=0..1-x),x=0..1);
`                                   I2 := 1/8`
 Now the left hand side. Let us split the surface into the four triangles: S = Sxy + Sxz + Syz + Sxyz. Where Sab is the triangle on the ab plane and Sxyz is the triangle on the plane x+y+z=1. Let us denote the surface integrals over Sxy, etc... by J1,J2,J3 and J4. We have, that the outward unit normal for Sxy is (-k) and thus,

> ;

```                                  /   /
|   |    3  3
J1 :=  |   |  -x  y  dS
|   |
/   /
Sxy```
 which evaluates to:

> J1 := int(int(dotprod(F,-k),y=0..1-x),x=0..1);

```                                         -1
J1 := ----
1120```
 Now compute the two central surface integrals J2 and J3. The outward unit normals for Sxz and Syz are j and i respectively so,

> ;

```                                    /   /
|   |
J2 :=  |   |  x y dS = 0
|   |
/   /
Sxz```
 this integral is 0 since y=0 on the xz plane. For J3 we have,

> J3 := Int(Int(dotprod(F,i),S),u);

```                                    /   /
|   |   2
J3 :=  |   |  x  dS = 0
|   |
/   /
Syz```
 again this integral is 0 since x=0 on the yz plane. For the last integral we see that the ouward unit normal is

> N := evalm(vector([1,1,1])/sqrt(3));

```                           [     1/2       1/2       1/2]
N := [1/3 3   , 1/3 3   , 1/3 3   ]```
 Notice that we need to divide by sqrt(3) (the norm of (1,1,1)) to get a UNIT normal. Hence,

> ;

```                /   /
|   |       2  1/2            1/2        3  3  1/2
J4 :=  |   |  1/3 x  3    + 1/3 x y 3    + 1/3 x  y  3    dS
|   |
/   /
Sxyz```
 Computing the dS for this last integral we find dS = sqrt(3) dx dy (Hint: parametrize the triangle with x and y and use the formula for dS).

> J4 := Int(Int(dotprod(F,[1,1,1]),y=0..1-x),x=0..1);

```                           1    1 - x
/    /
|    |        2          3  3
J4 :=  |    |       x  + x y + x  y  dy dx
|    |
/    /
0    0```
 which evaluates to:

> J4 := int(int(dotprod(F,[1,1,1]),y=0..1-x),x=0..1);

```                                        141
J4 := ----
1120```
 Finally combining the J's we get:

> I1 := J1+0+0+J4;

`                                   I1 := 1/8`
 The same as before!!

Link to the commands in this file
Carlos Rodriguez <carlos@math.albany.edu>