MATLAB Project 4 Due: M 12/11

Some of the following material has been adapted and copied with permission from material created by Justin Wyss-Gallifent.

The basic format of this guide is the same as the first two. The project involves creating a script m-file and put the tasks in order into that file, then publishing that file using MATLAB's publish to HTML feature and printing the resulting document. You should separate each task into cells with the title of the task preceded by a %%.

---

Fundamental Theorem of Line Integrals

We can use MATLAB to verify the fundamental theorem of line integrals. For example, suppose that we wanted to find $$ \int_C x\mathrm{d}x + y \mathrm{d}y + z \mathrm{d}z $$ over the curve $C$ defined by $\mathbf{r}(t) = t\mathbf{i} + 3t^2\mathbf{j} + (t-1)^3\mathbf{k}, 0\leq t\leq 1$. We can first check whether the vector field is conservative using the potential command.

syms x y z;
F = [x,y,z];
f = potential(F,[x,y,z]);

f=

x^2/2 + y^2/2 + z^2/2

If the vector field doesn't have a potential (i.e. it has a non-zero curl) then the function spits out NaN, for instance,

clear;
F = [x,y,x*z];
f = potential(F,[x,y,z]);

f=

NaN

In the case of the line integral above, we can use the fundamental theorem of line integrals to evaluate it.

clear;
syms x y z t;
F = [x,y,z];
f = potential(F,[x,y,z]);
r = [t,3*t^2,(t-1)^3];
P = subs(r,t,0);
Q = subs(r,t,1);
subs(f,[x,y,z],Q) - subs(f,[x,y,z],P)

ans=

9/2

We can compare this to the explicit integration

sub = subs(F,[x,y,z],r);
int(dot(sub,diff(r,t)),0,1)

ans=

9/2

---

Surface Integrals of Functions

The nice thing about having a surface in Matlab as a vector is that it’s easy to work with.

For example suppose we want to find $\iint_{\Sigma} xy \mathrm{d}S$ where $\Sigma$ is the part of $z = x^2 + y^2$ having $0\leq x\leq 2$ and $0\leq y \leq 3$.

As part of this procedure we need the magnitude of the cross products of the derivatives of r. Here is how we do just that bit:

clear all;
syms x y;
rbar = [x,y,x^2+y^2];
mylength = @(u) sqrt(u*transpose(u));
simplify(mylength(cross(diff(rbar,x),diff(rbar,y))))

ans =

(4*x^2 + 4*y^2 + 1)^(1/2)

Be aware that the first vector soln is giving you the two possible $x$ values. These pair up with the two possible y values. Thus in friendlier terms the two solutions are $(3, 1)$ and $(−1, −3)$. Also note that the solve command returns the solution in alphabetical order, meaning it returns $x$ and then $y$. This is why we assign the solution to [xsoln,ysoln]=.

So now here’s the really nice thing. To do our entire integral involves plugging the i, j, k components in for x, y, z, multiplying by that magnitude and integrate over the appropriate limits. Here it is all together:

clear all;
syms x y z;
rbar = [x,y,x^2+y^2];
f = x*y;
mylength = @(u) sqrt(u*transpose(u));
mag = simplify(mylength(cross(diff(rbar,x),diff(rbar,y))));
subresult = subs(f,[x,y,z],rbar);
int(int(subresult*mag,x,0,2),y,0,3)

ans =

(2809*53^(1/2))/240 - (1369*37^(1/2))/240 - (289*17^(1/2))/240 + 1/240

Make sure you read this carefully (after the clear all) to see what each line does. The second line sets the symbolic variables. The third line is the parametrization of the surface. The fourth line is the function to integrate. The fifth line creates a length function we’ll need later. The sixth line finds the magnitude of the cross product of the derivatives. The seventh line substitutes the components from the parametrization into the real-valued function we want to integrate. The eighth and final line does the double integral required.

---

Surface Integrals of Vector Fields

Similarly we can take the surface integral of a vector field. We only need to be careful in that Matlab can’t take care of orientation so we’ll need to do that and instead of needing the magnitude of the cross product we just need the cross product. Here is problem 6 from the 15.6 exercises.

clear
syms t x y z;
rbar = [cos(t),y,sin(t)];
F = [x,y,z];
kross = simplify(cross(diff(rbar,t),diff(rbar,y)));
sub = subs(F,[x,y,z],rbar);
int(int(dot(sub,kross),y,-2,1),t,0,2*pi)

ans =

-6*pi

Again read this carefully. After the clear the first line sets the symbolic variables. The second sets the parametrization and the third sets the vector field. The fourth finds the cross product of the derivatives. The fifth substitutes the parametrization into the vector field. The sixth does the double integral of the dot product as required for the surface integral of a vector field.

---

Stokes Theorem

We can explore and verify Stokes' Theorem using MATLAB. Suppose we wanted to compute the flux integral $$ \iint_\Sigma (xy\mathbf{i} + xz\mathbf{j} - zy\mathbf{j})\cdot \mathbf{n}\,\mathrm{d}S $$ where $\Sigma$ is the top half of the sphere of radius 1, oriented with upward facing normal. If we knew that $\mathbf{F}= \nabla \times \mathbf{A}$ for some vector field $\mathbf{A}$ then we could apply Stokes Theorem to turn this flux integral into a line integral around the circle $x^2 + y^2 =1$. Any vector field $\mathbf{A}$ satisfying $\mathbf{F} = \nabla\times \mathbf{A}$ is known as a vector potential. As is turns out, similar to scalar potentials for conservative vector fields, MATLAB has a command vectorPotential which can find a vector potential of $\mathbf{F}$ (assuming it exists). For instance, for the vector field in the above flux integral:

clear
syms x y z;
F = [x*y,x*z,-z*y];
vectorPotential(F, [x y z])

ans =

 (x*z^2)/2
    -x*y*z
         0

If there is no vector potential (i.e. if $\nabla\cdot \mathbf{F} \neq 0$), then MATLAB returns NaN

clear
syms x y z;
F = [x,y,z];
vectorPotential(F, [x y z])

ans =

 NaN
 NaN		    
 NaN

We can of course use the vector potential to find flux integrals by Stokes Theorem $$ \iint_{\Sigma} \mathbf{F}\cdot \mathbf{n}\, \mathrm{d}S = \int_C \mathbf{A}\cdot \mathrm{d}\mathbf{r}. $$ For the flux integral in the begining of this section we can compute the associated line integral

clear
syms x y z t;
F = [x*y,x*z,-z*y];
A = vectorPotential(F, [x y z]);
r = [cos(t) sin(t) 0];
int(dot(subs(A, [x y z], r),diff(r,t)),t,0,2*pi)

ans =

 0

Of course, we can verify that we get the same answer by computing the flux integral directly

clear
syms x y z t p;
F = [x*y,x*z,-z*y];
r = [sin(p)*cos(t) sin(p)*sin(t) cos(p)];
nds = simplify(cross(diff(r,t),diff(r,p)));
int(int(dot(subs(F, [x y z], r),nds),p,0,pi/2),t,0,2*pi)
		    

ans =

 0