• The Net Charge on a Sphere
• The Mass and Center of Mass of a Rubber Sheet

Surface Integrals of Scalar Fields

Using Maple and the vec_calc Package

This worksheet shows how to compute surface integrals of scalar fields using Maple and the vec_calc package. As examples we compute

* The Net Charge on a Sphere

* The Mass and Center of Mass of a Rubber Sheet

To start the vec_calc package, execute the following commands:

> restart;

> libname:="C:/mylib/vec_calc7", libname:

> with(vec_calc): vc_aliases:

> with(linalg):with(student):with(plots):

Warning, the protected names norm and trace have been redefined and unprotected

Warning, the name changecoords has been redefined

>

The Net Charge on a Sphere

The charge distribution on the surface of a semiconductor sphere of radius is given by where is the polar angle measured down from the -axis and is the azimuthal angle measured counterclockwise from the positive -axis. We want to find the net charge on the sphere.

We first make some assumptions on the variables which will facilitate the integrations:

> interface(showassumed=0);

> assume(0<=phi,phi<=Pi,0<=theta,theta<=2*Pi);

> about(phi,theta);

Originally phi, renamed phi~:

  is assumed to be: RealRange(0,Pi)

Originally theta, renamed theta~:

  is assumed to be: RealRange(0,2*Pi)

Then we define the density:

> density:=2*sin(theta)+3*cos(phi)*cos(theta)-2*sin(phi);

Some places the charge is positive and some places the charge is negative. To visualize this we make a spherical plot of the density. We take the radius to be 20 plus the density so that the radius will be positive. (The number 20 is arbitrary.) Then the charge is positive when the radius is bigger than 20, and negative when the radius is less than 20.

> dens_plot:=sphereplot(20+density,theta=0..2*Pi,phi=0..Pi, color=density, scaling=constrained, orientation=[15,75]): dens_plot;

We also set the color equal to the density. Then the color is blue, violet or red when the density is positive, and the color is green, yellow, orange or red when the density is negative. It is still not easy to see where the charge is positive or negative. So we plot a gray spherical grid and display it with the density:

> sph_plot:=sphereplot(20,theta=0..2*Pi,phi=0..Pi, color=grey, style=wireframe):

> display({sph_plot,dens_plot}, scaling=constrained, orientation=[15,75]);

The charge is positive when the plot is above the grid and negative when the plot is below the grid. Try rotating the graph by clicking in the plot and dragging the mouse.

We are now ready to compute the net charge which is the surface integral of the density: . The Jacobian is the length of the coordinate normal which we now proceed to compute. Since the radius is 15, the spherical coordinate system on the sphere is

> Rsph:=[15*sin(phi)*cos(theta),15*sin(phi)*sin(theta),15*cos(phi)];

The tangent vectors are:

> R[theta]:=diff(Rsph,theta);

> R[phi]:=diff(Rsph,phi);

So the normal vector is

> N:= R[theta] &x R[phi];

and the length of the normal is

> lN:=len(N); simplify(%);

Maple does not simplify this to the simplest form. So we enter it by hand:

> lN:=15^2*sin(phi);

So the net charge is

> Muint(density*lN,theta = 0 .. 2*Pi,phi = 0 .. Pi); Charge:=value(%);

Another way to compute this surface integral is to use the Surface_int_scalar command (or its alias Sis ) from the vec_calc package which works directly with the parametrized surface and the scalar function. The surface must be a function of its parameters:

> Rsphf:=MF([theta,phi],Rsph);

and the scalar must be a function of position. Looking at the formula for the density, we convert to spherical coordinates:

> densityf:=MF([x,y,z], 2*y/sqrt(x^2+y^2) + 3*z/sqrt(x^2+y^2+z^2)*x/sqrt(x^2+y^2) - 2*sqrt(x^2+y^2)/sqrt(x^2+y^2+z^2));

We check the restriction to the sphere agrees with the given function:

> dens2:=simplify(densityf(op(Rsph))); density;

I had some difficulty getting Maple to show that these two are equal, but this works:

> simplify(dens2*sqrt(1-cos(phi)^2)-density*sin(phi));

We are now ready to compute the integral for the charge:

> Sis(densityf,Rsphf,theta = 0 .. 2*Pi,phi = 0 .. Pi); Charge:=value(%);

>

The Mass and Center of Mass of a Rubber Sheet

A rubber sheet is stretched across a circular ring of radius and a weight is placed in the middle which pulls it downward. Consequently the sheet has the shape of the graph of for . Since it has been stretched down, its density is which is less in the middle and more at the edge. We want to find the mass and center of mass of the rubber sheet.

We start by parametrizing the sheet in cylindrical (or polar) coordinates:

> R:=MF([r,theta],[r*cos(theta),r*sin(theta),2-2^(2-r^2/9)]);

and plot it

> plot3d(R(r,theta), r=0..3, theta=0..2*Pi, scaling=constrained, style=PATCHCONTOUR, shading=ZHUE, orientation=[0,70]);

Then we define the density:

> rho:=MF([x,y,z], 3+(x^2+y^2)/9);

The mass is . So we next evaluate the density on the curve:

> rhoR:=simplify(rho(op(R(r,theta))));

Then we compute the tangent and normal vectors:

> Ru:=diff(R(r,theta),r);

> Rtheta:=diff(R(r,theta),theta);

> N:=Ru &x Rtheta;

So the length of the normal is

> lN:=simplify(len(N));

Thus the mass is

> M:=Muint(rhoR*lN, r=0..3, theta=0..2*Pi); value(%); M:=evalf(%);

So the mass is about 130 .

By symmetry the - and -components of the center of mass are zero. The -component is the -moment divided by the mass. The -moment is found by adding a factor of (in terms of the curve) inside the integral. The center of mass is then found by dividing the first moments by the mass.

> Muint(rhoR*lN*R(r,theta)[3], r=0..3, theta=0..2*Pi); zmom:=evalf(%);

> zcm:=zmom/M;

So the center of mass is at

> CM:=[0,0,zcm];

Another way to compute this surface integral is to use the Surface_int_scalar command (or its alias Sis ) from the vec_calc package which works directly with the parametrized surface and the scalar function. The mass is:

> Sis(rho,R, r=0..3, theta = 0 .. 2*Pi); M:=evalf(%);

To find the -moment we need a function whose value is :

> z0:=MF([x,y,z],z);

So the -moment is:

> Sis(rho*z0,R, r=0..3, theta = 0 .. 2*Pi); zmom:=evalf(%);

and the -component of the center of mass is:

> zcm:=zmom/M;