Illustrating Solutions to 2-D Partial Differential Equations

Matthew M. Thomas
Department of Chemical Engineering
Washington University in Saint Louis

Engineering students are usually introduced to partial differential
equations (PDEs) during their sophomore years of college. The first
PDEs they often encounter are the one- and two-dimensional wave
and heat equations [Sokolnikoff and Redheffer, 1966]. Mathematica
can illustrate the solutions to these PDEs, for a number of initial
and boundary conditions.  Such illustrations bring life to the
tedious mathematical expressions that comprise these solutions.
The illustrations assist the student in visualizing solutions to
these PDEs, thereby enhancing his understanding of them.

Mathematica functions Plot3D and ListPlot3D can create the
three-dimensional graphs that, when animated, represent solutions
to the 2-D wave and heat PDEs.  The 2-D wave PDE, wtt = c2(wxx +
wyy), describes an oscillating membrane.  The Mathematica function
W[x_, y_, t_] in Figure 1 is the solution to the 2-D wave PDE for
a freely oscillating, uniform square membrane (1 square meter) with
surface tension: density ratio c = 1 meter/second, initial conditions

	w(x, y, 0) = x y (1-x) (1-y) 
and 
	wt(x, y, 0) = - x y (1-x) (1-y), 

and boundary conditions 

    w(x, 0 m, t) = w(x, 1 m, t) = w(0 m, y, t) = w(1 m, y, t) = 0. 

In[1]:= 
W[x_, y_, t_] := 
   (64/N[Pi]^6) *
   N[Sum[(Cos[N[Pi] t Sqrt[m^2 + n^2]] -
   (Sin[N[Pi] t Sqrt[m^2 + n^2]]/
   (N[Pi] Sqrt[m^2 + n^2])))*
   Sin[m N[Pi] x] Sin[n N[Pi] y]/(m^3 n^3),
   {m, 1, 15, 2}, {n, 1, 15, 2}]];
                 
GrayPlot[f_, ranges__, n_] :=
   Plot3D[f, ranges,
      Mesh -> True,
      PlotRange -> {-0.07, 0.07},
      PlotPoints -> n];
      
WavePlot[n_, t_] := GrayPlot[W[x,y,t], {x,0,1}, {y,0,1}, n]
WavePlot[20, 0.0]

Out[1]=
   -Graphics-

Figure 2 illustrates this solution for a 20 x 20 grid at t = 0.15.
Using Do[WavePlot[20, t], {t, tstart, tend, tdelta}] creates the
sequence which, when animated, illustrates the vibrating membrane.

The 2-D heat PDE, vt = a(vxx + vyy ), describes heat flow in the
semi-infinite direction of a Rsemi-infinitelyS long solid. The
Mathematica function T[x_, y_, t_] := U[x,t]* U[y, t], as optimized
in [Sullivan and Thomas, 1990],  and presented in a Rblack-and-whiteS
version here, is the solution to the 2-D heat PDE for a semi-infinite
solid with a square cross-section and thermal diffusivity a = 1,
having initial condition v(x, y, 0) = 1, and boundary conditions
v(x, -1, t) = v(x, 1, t) = v(-1, y, t) = v(1, y, t) = 0 for t > 0.
Figure 3 illustrates this solution for a 50 x 50 grid at t = 0.10.
Using Do[HeatPlot[50, t], {t, tstart, tend, tdelta}] creates the
sequence for illustrating the heat flow.

In[3]:=
U[x_, t_] := 
     Block[{k=Ceiling[1/Sqrt[t]]},
      2 Sum[Block[{nu=(m+.5) N[Pi]},
         (-1)^m/nu Exp[-nu^2 t] Cos[nu x]],
        {m, 0, k}]];
                  
U[x_, 0.] := 1;

ListHuePlot[array_]:= ListPlot3D[
   Table[Table[0, {Length[array]+1}],
           {Length[array]+1}],
   Map[GrayLevel[.75 #]&, array, {2}],
   Boxed -> False, ViewPoint -> {0, 0, 10},
   Mesh -> False, Axes -> None];
        
HeatPlot[n_, t_] := 
Block[{u=Table[U[x, t], {x, -1, 1, 2/n}]},
    ListHuePlot[Outer[Times, u, u]]];
HeatPlot[50, 0.1]

Out[3]=
-Graphics-

When animating the 2-D wave PDE solution, we illustrate an oscillating
membrane.  When animating the 2-D heat PDE solution, we illustrate
a stationary cross section whose temperature profile varies with
time.  As can be seen from Figure 3, the  function T[x_, y_, t_]
varies from 0 (hottest) to 1 (coldest). In the GrayLevel command,
gray level varies from 0 (hottest; black) to 0.75 (coldest; light
gray). Black and light gray correspond to hot and cold temperatures,
respectively.  If we increase gray level in inverse proportion to
temperature, we intuitively represent temperature by gray scale in
Mathematica.

For educational purposes, we can create a library of Mathematica
Notebooks that illustrate solutions to 1-D and 2-D wave and heat
PDEs, for a number of initial and boundary conditions, and for a
number of membrane and solid cross section geometries.  Two-dimensional
Mathematica plots can represent a solution to a 1-D PDE for a given
time; 3-D plots can represent a solution to a 2-D PDE solution for
a given time; animation of these plots can demonstrate the
time-dependent behavior of these solutions.  This library would
provide a needed aid for the average sophomore engineering student
who often finds illustration of physical phenomena fairly difficult
to visualize.


References

Sokolnikoff, I.S., and R.M. Redheffer. 1966. Mathematics of Physics
and Modern Engineering, 2nd ed., New York, McGraw-Hill.

Sullivan, J.M., and M.M. Thomas. 1990. The Mathematica Journal,
1(1), pages 80-84.

About the Author

Matthew M. Thomas received B.S. degrees in Chemical Engineering
(1985) and Data Processing (1988) from Washington University, and
an M.S. in Chemical Engineering (1989) from the same university.
He is now enrolled in the Chemical Engineering doctoral program at
Washington University, where his research interests include
computer-aided chemical engineering and the application of artificial
neural networks to chemical process control.

Department of Chemical Engineering
Washington University in Saint Louis
Saint Louis, MO  63130-4899
email:  thomas@wuche1.wustl.edu