33hotSfer.html

Partial Differential Equations PowerTool

by Dr. Jim Herod

Section 7.5: Warm Spheres

Maple Packages needed in Section 7.5

> restart:

> with(plots):

Warning, the name changecoords has been redefined

> with(orthopoly);

>

Another important coordinate system is that of a sphere. In that coordinate system we have

is the distance from the origin,

is the angle in the x - y plane; that is, it measures longitude.

is the angle from the top; that is, it measures latitude.

Thus, = 1, 0 < < 2 , 0 < < is a sphere with radius 1.

The connection with rectangular coordinates is

x = sin( ) cos( )

y = sin( ) sin( )

z = cos( )

> plot3d(1,theta=0..2*Pi,phi=0..Pi, coords=spherical, style=wireframe,axes=NORMAL,orientation=[35,70]);

Also, = 1, = , 0 < < is a half ring running from the north pole to the south pole of the sphere. I draw the fat so you can see it. You might experiment by making r just 1 to see what I avoided.

> plot3d([r,Pi/4,phi],r = 1..1.2,phi=0..Pi, coords=spherical,axes=NORMAL,orientation=[15,75],style=patchnogrid,scaling=constrained);

Finally, = 1, = 49 degrees, defines a part of the boundary between Western Canada and Western United States.

In this coordinate system, the Laplacian Operator is

u = { + + }

We solve 0 = u, u(1, phi) = f( ), where u is independent of . From the assumption that

u( , ) = R( ) ,

we are led to this separation of variables situation:

( R '( ) ) ' / R + ( sin( ) ' ) '/ sin( ) = 0.

This suggests the ordinary differential equations

( R '( ) ) ' - R( ) = 0, 0 < < 1,

( sin( ) ' ) ' + sin( ) = 0, 0 < < .

Neither equation has a boundary condition. We have conditions of boundedness. In the second equation, we take x = cos( ) changing the equation as follows:

> x:=phi->cos(phi);

> Phi:=phi->y(x(phi));

> diff(sin(phi)*diff(Phi(phi),phi),phi);

Hence, ( sin( ) ' ) ' = y ''(x) - 2 sin( ) cos( ) y '(x).

The second differential equation above becomes

y ''(x) - 2 cos( ) y '(x) + y(x) = 0.

In terms of x alone,

(1 - ) y ''(x) - 2 x y '(x) + y(x) = 0, -1 < x < 1.

We examine solutions of this differential equation.

> x:='x':

> dsolve((1-x^2)*diff(y(x),x,x)- 2*x*diff(y(x),x)+mu^2*y(x)=0,y(x));

>

There is the suggestion that we consider Legendre Polynomials. We digress to recall these functions.

A recollection of Legendre Polynomials

Here are three ways to conceive of the Legendre Polynomials -- four, if we include Maple.

Method 1: solve the differential equation.

> dsolve((1-x^2)*diff(y(x),x,x)- 2*x*diff(y(x),x)+n*(n+1)*y(x)=0,y(x));

> plot([LegendreP(1,x),LegendreQ(1,x)],x=-1..1,color=[red,black]);

We got only the first graph. The following shows that the second graphs are not defined on the interval of interest.

> plot([LegendreQ(1,x),LegendreQ(2,x),LegendreQ(3,x)],x=0..5,y=0..1/2);

>

Here are the first three Legendre Polynomials with an easier calling.

> plot([P(0,x),P(1,x),P(2,x)],x=-1..1);

Here is a check that, for example, the 5th one satisfies the differential equation.

> y:=x->P(5,x);
(1-x^2)*diff(y(x),x,x)- 2*x*diff(y(x),x)+5*(5+1)*y(x);
simplify(%);

>

Method 2: We could apply the Gramm Schmidt Process to 1, x, , ... .

We examined these ideas earlier. Here, we illustrate that the Legendre polynomials are orthogonal.

> int(P(3,x)*P(5,x),x=-1..1);

>

Method 3. We could generate these by taking the appropriate derivatives.

We show that is the polynomial.

> f:=(n,x)->(x^2-1)^n;
Q:=(n,x)->1/(2^n*n!)*diff(f(n,x),x$n);

> for n from 1 to 3 do
expand(Q(n,x)),`and`, P(n,x);
od;

>

Method 4: Use the recursion formulas (n+1) + n = (2 n + 1) x.

With this method, we assume you know the first two.

> P3[0]:=1;
P3[1]:=x;
for n from 1 to 3 do
P3[n+1]:=expand(1/(n+1)*((2*n+1)*P3[n]*x-n*P3[n-1]));
od;
n:='n';

> for n from 0 to 4 do
P3[n]/P(n,x);
od;
n:='n';

>

Observation 1: if 0 < m < n, then = 0. We check this for a special case.

> int(x^3*P(5,x),x=-1..1);

>

Observation 2: We have a formula for the norm of P(n,x): = 2/(2 n + 1).

We check three cases by comparing the integral with 2/(2 n + 1).

> for n from 0 to 5 do
int(P(n,x)^2,x=-1..1)-2/(2*n+1);
od;
n:='n';

>

Observation 3: We can make polynomial approximations for functions on [-1, 1]. We approximate cos( x ) with the first three Legendre polynomials.

> for n from 0 to 2 do
a[n]:=int(cos(Pi/2*x)*P(n,x),x=-1..1)*(2*n+1)/2;
od;
n:='n';

> plot([a[0]+a[1]*P(1,x)+a[2]*P(2,x),cos(Pi/2*x)],x=-1..1,
color=[black,red]);

>

This completes our discourse on Legendre Polynomials. We return to the partial differential equation.

Recall where we were. The partial differential equation led to two ordinary differential equations.

The second differential equation became

y ''(x) - 2 cos( ) y '(x) + y(x) = 0.

In terms of x alone,

(1 - ) y ''(x) - 2 x y '(x) + y(x) = 0, -1 < x < 1.

This led to

(1 - ) y ''(x) - 2 x y '(x) + y(x) = 0, -1 < x < 1

and Legendre Polynomials.

Here is the first of the differential equations from above, with this :

( R '( ) ) ' - R( ) = 0, 0 < < 1,

R '' + 2 R ' - R( ) = 0, 0 < < 1.

We solve this.

> dsolve(r^2*diff(R(r),r,r)+2*r*diff(R(r),r)-n*(n+1)*R(r)=0,R(r));

>

We see that the only bounded solutions are . Hence, we are ready to make up the general solution for the partial differential equation.

General Solution:

First, we verify that products of solutions are solutions.

> n:=4;
u:=(r,phi)->r^n*P(n,cos(phi));

> diff(r^2*diff(u(r,phi),r),r)+
1/sin(phi)*diff(sin(phi)*diff(u(r,phi),phi),phi):
simplify(%);

>

I checked this with n = 4. You check it at other places. What about sums?

> n:='n';
u:=(r,phi)->sum(a[n]*r^n*P(n,cos(phi)),n=0..4);

> diff(r^2*diff(u(r,phi),r),r)+
1/sin(phi)*diff(sin(phi)*diff(u(r,phi),phi),phi):
simplify(%);

>

We are now ready to compute the coefficients for a boundary condition. We solve

0 = u with 0 < < 1, 0 < < , 0 < < 2 with boundary condition u(1, ) = f( ).

The coefficients will be

= .

Let's take the special case that f( ) = cos( ).

> for n from 0 to 4 do
a[n]:=(2*n+1)/2*int(cos(phi)*P(n,cos(phi))*sin(phi),phi=0..Pi);
od;

Thus, the solution for this problem is
u(r, ) = r* P(1,cos( )) = r cos( ).

How can we illustrate what we have?

(1) On each cross sectional plane parallel to the x-y plane, u has the same value, namely u(r, ) = z. We draw a portion of the sphere where all the points have the same temperature.

> addcoords(z_cylindrical,[z,r,theta],[r*cos(theta),r*sin(theta),z]);
dis:=plot3d(1/2,r=0..0.86,theta=0..2*Pi,coords=z_cylindrical, orientation=[35,70],axes=normal):
sph:=plot3d(1,theta=0..2*Pi,phi=0..Pi, coords=spherical, style=wireframe,axes=NORMAL,orientation=[35,70]):
display3d({dis,sph});

(2) If we hold r fixed between 0 and 1, and if we let go from 0 to , we get imbedded spheres. Here are three of these where only half of the middle one is drawn so that the innermost can be seen.

> S1:=plot3d(1,theta=0..2*Pi,phi=0..Pi, coords=spherical, style=wireframe,axes=NORMAL,orientation=[-35,70]):
S2:=plot3d(3/4,theta=0..Pi,phi=0..Pi, coords=spherical,axes=NORMAL,orientation=[-35,70]):
S3:=plot3d(1/2,theta=0..2*Pi,phi=0..Pi, coords=spherical,axes=NORMAL,orientation=[-35,70],color=black):
display({S1,S2,S3});

The temperature on each such sphere varies from the hottest point at the top and the coldest at the bottom. Here are graphs of u on the three embedded spheres as changes from 0 to ... from the north pole to the south pole.

> u:=(r,phi)->r*P(1,cos(phi));

> plot([u(1,phi),u(3/4,phi),u(1/2,phi)],phi=0..Pi,color=[black,red,green]);

In this Section, we have made analytic what we know intuitively: if we know the temperature distribution on the surface of a sphere which has no sources or sinks, then we can determine the temperature distribution on the inside of the sphere.

EMAIL: herod@math.gatech.edu or jherod@tds.net

URL: http://www.math.gatech.edu/~herod