Application Center - Maplesoft

App Preview:

Electric Charge on a Sphere and the Mass of a Rubber Sheet

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application


 

surfintscal.mws

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 15*cm is given by density := [2*sin(theta)+3*cos(phi)*cos(theta)-2*si... coul/(cm^2) where phi is the polar angle measured down from the z -axis and theta is the azimuthal angle measured counterclockwise from the positive x -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);

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

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;

[Maple 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]);

[Maple Plot]

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: Int(density,S) = Int(Int(density*Jacobian,theta = 0... . 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)];

Rsph := [15*sin(phi)*cos(theta), 15*sin(phi)*sin(th...

The tangent vectors are:

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

R[theta] := [-15*sin(phi)*sin(theta), 15*sin(phi)*c...

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

R[phi] := [15*cos(phi)*cos(theta), 15*cos(phi)*sin(...

So the normal vector is

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

N := [-225*sin(phi)^2*cos(theta), -225*sin(phi)^2*s...

and the length of the normal is

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

lN := sqrt(50625*sin(phi)^4*cos(theta)^2+50625*sin(...

225*sqrt(1-cos(phi)^2)

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

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

lN := 225*sin(phi)

So the net charge is

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

Int(Int(225*(2*sin(theta)+3*cos(phi)*cos(theta)-2*s...

Charge := -450*Pi^2

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);

Rsphf := [proc (theta, phi) options operator, arrow...

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));

densityf := proc (x, y, z) options operator, arrow;...

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

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

dens2 := (2*sin(phi)*sin(theta)+3*cos(phi)*sin(phi)...

2*sin(theta)+3*cos(phi)*cos(theta)-2*sin(phi)

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));

0

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

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

Int(Int(450*sin(phi)*sin(theta)+675*cos(phi)*sin(ph...

Charge := -450*Pi^2

>

The Mass and Center of Mass of a Rubber Sheet

A rubber sheet is stretched across a circular ring of radius R = 3 cm and a weight is placed in the middle which pulls it downward. Consequently the sheet has the shape of the graph of z = 2-2^(2-(x^2+y^2)/9) for x^2+y^2 <= 9 . Since it has been stretched down, its density is rho = 3+(x^2+y^2)/9 gm/(cm^2) 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)]);

R := [proc (r, theta) options operator, arrow; r*co...

and plot it

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

[Maple Plot]

Then we define the density:

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

rho := proc (x, y, z) options operator, arrow; 3+1/...

The mass is M = Int(rho,S) ` ` = Int(Int(rho(R)*abs(N),r = 0 .. 3),theta = 0 .... . So we next evaluate the density on the curve:

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

rhoR := 3+1/9*r^2

Then we compute the tangent and normal vectors:

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

Ru := [cos(theta), sin(theta), 2/9*2^(2-1/9*r^2)*r*...

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

Rtheta := [-r*sin(theta), r*cos(theta), 0]

> N:=Ru &x Rtheta;

N := [-2/9*2^(2-1/9*r^2)*r^2*ln(2)*cos(theta), -2/9...

So the length of the normal is

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

lN := 1/9*sqrt(r^2*(64*r^2*2^(-2/9*r^2)*ln(2)^2+81)...

Thus the mass is

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

M := Int(Int(1/9*(3+1/9*r^2)*sqrt(r^2*(64*r^2*2^(-2...

2*int(1/9*(3+1/9*r^2)*sqrt(r^2*(64*r^2*2^(-2/9*r^2)...

M := 129.9875449

So the mass is about 130 gm .

By symmetry the x - and y -components of the center of mass are zero. The z -component is the z -moment divided by the mass. The z -moment is found by adding a factor of z (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(%);

Int(Int(1/9*(3+1/9*r^2)*sqrt(r^2*(64*r^2*2^(-2/9*r^...

zmom := -104.2681112

> zcm:=zmom/M;

zcm := -.8021392456

So the center of mass is at

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

CM := [0, 0, -.8021392456]

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(%);

Int(Int(1/9*(3+1/9*r^2)*r*sqrt(64*r^2*2^(-2/9*r^2)*...

M := 129.9875449

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

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

z0 := proc (x, y, z) options operator, arrow; z end...

So the z -moment is:

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

Int(Int(1/9*(6-12*2^(-1/9*r^2)+2/9*r^2-4/9*r^2*2^(-...

zmom := -104.2681112

and the z -component of the center of mass is:

> zcm:=zmom/M;

zcm := -.8021392456