Numerical Solutions to PDE Boundary Value Problems in Maple 8

Copyright 2002 Waterloo Maple Inc.

Maple 8 can now compute numerical solutions to linear PDE systems over rectangular regions.  One application of this feature is the solution of classical boundary-value problems from physics, such as the heat conduction equation and the wave equation.

Example 1: Heat flow in a rod

Let's compute the temperature  over a rod, where t is time and x  is distance along the rod. We assume the ends of the rod maintain a constant temperature.  For each instance of the problem, we must specify the initial heat distribution  and the thermal diffusivity  of the rod.  

We use the time-dependent, one-dimensional heat equation in cartesian coordinates:

,

subject to the boundary conditions

,
, and

The Laplacian  function in Maple 8's new VectorCalculus  package is useful for writing down 2nd-order PDEs, so we first load that package.

>	with(VectorCalculus);

Warning, the assigned names <,> and <|> now have a global binding

Warning, these protected names have been redefined and unprotected: *, +, ., Vector, diff, int, limit, series

We write down the heat equation using the Laplacian  function. We assume

>	Heat_Eqn := { diff( u(x,t), t) = Laplacian( u(x,t), 'cartesian'[x]) / 10 };

Next we write down the boundary conditions:  The ends of the rod are kept at a constant temperature (0), and we assume the initial temperature distribution is

>	BCs:={u(x,0) = (1-x)*(1-exp(-50*x)), u(0,t) = 0, u(1,t) = 0};

Now we solve the PDE boundary-value problem numerically with the pdsolve  command and numeric  option specified.  We can set the accuracy of the solution by specifying the time step and space step of the discretization over the distance-time rectangle.  

>	Heat_Solution := pdsolve( Heat_Eqn, BCs, numeric, timestep=1/100, spacestep=1/100);

Note that Maple delivered the solution as a Maple module.  The module provides functions for viewing the solution as plots ( plot  and plot3d ), animations over time ( animate ) and numerical values ( value ).

Here's an animation of the temperature distribution

>	Heat_Solution :- animate( u(x,t) ,t=0..1.5, frames=80, labels=["x", "u(x,t)"], labelfont=[TIMES,ROMAN,14]);

And a plot of the solution space.

>	Heat_Solution :- plot3d(u(x,t), t=0..2, shading=zhue, axes=boxed, labels=["x","t","u(x,t)"], labelfont=[TIMES,ROMAN,16]);

To get numerical values, we can create a procedure using the value  procedure from the solution's module.

>	RodTemperature := Heat_Solution :- value();

The temperature at ,

>	RodTemperature( .5, 3 );

Let's model a more complicated physical system.  Here the right end of the rod is more dissipative than the left, and the rod's dissipativity decreases with time.

>	Heat_Eqn := {diff(u(x,t),t)=exp(-2*t)*(0.01+5*x)/10*diff(u(x,t),x,x)};

>	BCs:={u(x,0)=sin(Pi*x), u(0,t)=0, u(1,t)=0};

>	Heat_Solution := pdsolve( Heat_Eqn, BCs, numeric, timestep=1/100, spacestep=1/100);

>	Heat_Solution :- animate( u(x,t) ,t=0..3.5, frames=50, labels=["x", "u(x,t)"], labelfont=[TIMES,ROMAN,14]);

>

Example 2: Vibration of a cord

Let's compute the vertical displacement  of a cord with length 1 and fixed endpoints, where t is time and x  is distance along the cord. For each instance of the problem, we must specify the initial displacement of the cord , the initial speed of the cord  and the horizontal wave speed c.

We use the one-dimensional wave equation in cartesian coordinates:

,

subject to the boundary conditions

,
,

, and

We write down the wave equation using the Laplacian  function with .

>	Cord_Eqn := { diff( u(x,t), t, t) = Laplacian( u(x,t), 'cartesian'[x])/10 };

We specifiy that the cord is fixed at both ends and has 0 initial speed.

>	BC1 := {u(0,t) = 0, u(1,t) = 0, D[2](u)(x,0)=0};

Case 1:  Cord is plucked from the center

>	BC2 := {u(x,0)=piecewise( x >= 0 and x <= 1/2, x, x > 1/2 and x <= 1, 1-x)};

>	sol:=pdsolve( Cord_Eqn, BC1 union BC2, numeric, spacestep=1/200, timestep=1/100):

>	sol:-animate( u(x,t) ,t=0..2*Pi, frames=30, labels=["x","u(x,t)"], labelfont=[TIMES,ROMAN,14], scaling=constrained);

Case 2:  Cord starts in the shape of a trigonometric function

>	BC2 := {u(x,0)=sin(4*Pi*x)/2};

>	sol := pdsolve( Cord_Eqn, BC1 union BC2, numeric, spacestep=1/200, timestep=1/100);

>	sol:-animate( u(x,t), t=0..Pi/2, frames=30,               labels=["x","u(x,t)"], labelfont=[TIMES,ROMAN,14],               scaling=constrained);

You can also display functions of the solution, or several functions of the solution on the same set of axes .

>	sol:-animate( [[u(x,t) + sin(120*x)/10, color=red],                [-u(x,t)/2 + sin(120*x)/20, color=blue]], t=0..Pi/2, frames=30,               labels=["x","u(x,t)"], labelfont=[TIMES,ROMAN,14],               scaling=constrained);

To learn more about Maple 8's numeric PDE solver, see the help pages at: