Introduction
The last two issues of CACHE News have
included articles demonstrating how partial differential equations can be
solved. Cutlip and Shacham (1998) used their own Polymath package to carry out
the necessary computations using the Method of Lines (MOL), while Rosen (1999)
has shown how the Excel spreadsheet can be used for other finite difference
computations. Our purpose here is to show how easy it can be solve such
problems using the computer algebra system Maple V.
The Method of Lines (MOL) is
a versatile approach to obtaining numerical solutions to partial differential
equations. The method has been made popular in science and engineering largely
through the work of W.E. Schiesser and coworkers (see, for example, Schiesser,
1991; Silebi and Schiesser, 1992).
MOL can be used to deal with a wide
variety of different types of PDE systems involving various types of PDE and/or
boundary conditions. The fundamental idea behind MOL is that the partial
derivatives in any spatial direction are replaced by finite difference
approximations, leaving only the time derivatives (or first derivatives in one
spatial direction if the equation is independent of time). Application of MOL
to a PDE may lead to a system of purely ordinary differential equations (ODEs).
Thus, all that may be needed to solve these equations is a good numerical ODE
solver. ODE solvers are widely available including in the computer algebra
system Maple V. Quite often, however, application of MOL leads to a mixed
system of differential and algebraic equations or DAEs. While DAE systems
sometimes can be solved by converting them to ODE systems, or by tricking an
ODE solver into thinking that they are ODEs, it makes sense to solve the DAE
system using a DAE solver. The availability of BESIRK, a DAE solver for Maple
(Schwalbe et al., 1996), makes implementation of MOL in Maple very easy indeed.
In what follows we demonstrate the use of
Maple by considering the same example used to illustrate the article of Cutlip
and Shacham (1998). Actual Maple code (and results) are provided for a problem
in which we find the temperature profiles in a slab as a function of space and
time.
Example
Heat transfer in a slab is described by
the following PDE:
> restart;
> PDE :=
diff(T(x,t),t) =
>
alpha*diff(T(x,t),x,x): PDE;
T is the temperature and a is the thermal diffusion coefficient. Initially, the slab
is at a uniform temperature (T0
).
> IC := t=0,T(x,t)=T[0]: IC;
At t = 0 the temperature at one
side is changed (to a value not yet specified, we shall call it (TA ).
> BC1:= x=0,T(x,t)=T[A]: BC1;
The other side is insulated and the temperature gradient is
zero.
> BC2 := x=L,diff(T(x,t),x)=0:
BC2;
We are going to use the method of lines to solve this
problem. The PDE applies to a sequence of values of T which are identified by
subscript i. The derivatives of T with
respect to x are to be replaced by second order finite difference
approximations. The conversion to appropriate form is accomplished with the `convert/fddiff` command that is part of this authors fdpack package and which we now load into Maple (along with some
other utilities).
> read
`e:/maple/numerics/pde/fdpack.mpl`:
> read
`e:/maple/utils/utils.mpl`:
> read
`e:/maple/thermo/plots/tplot.mpl`:
`convert/fddiff` can derive finite difference
approximations of any order, and with any (appropriate) number of forward and
backward steps. In this particular case the temperature at each node is to
remain a function of temperature only, and we only replace the spatial
derivatives with a second order difference approximation. This is accomplished
by including a zero in the order argument and none in the
indexletters argument to `convert/fddiff` (in the second position here since t appears in the second
position in the arguments to temperature in the original PDE).
> PDE1 := convert(PDE,fddiff,
> order=[2,0], forward=[1,0],
> indexletters=[i,none]): PDE1;
The above equation holds for all interior points, 1 < i < n, where n is the number of grid
points. It is not valid for the boundary lines (i = 1, i = n) since that
would require us to include fictitious points (i = 0, i = n+1). The next step is to specify n.
> n:=11;
The line spacing is
> hspec := h[x]=L/(n-1): hspec;
The first boundary condition may be expressed as
> BC1a := convert(BC1[2],fddiff,
> order=[2,0], forward=[0,0],
> indexletters=[1,none]): BC1a;
The insulated surface boundary condition becomes
> BC2a:=convert(BC2[2],fddiff,
> order=[2,0],forward=[0,0],
> indexletters=[n,none]): BC2a;
The initial condition becomes
> inclist := [t=0,seq(T[m]=T[0],
> m=1..n)]: inclist;
We make a list of all the equations:
> eqnlist := [BC1a,seq(PDE1,
> i=2..n-1),BC2a]: eqnlist;
… much more - similar - output omitted …
This is a differential algebraic equation (DAE) system, not
a purely ODE system. We can integrate this system using BESIRK which we now
read into Maple.
> read
`e:/maple/numerics/integ/BESIRK`:
The parameters in the model are given
numerical values:
> params := {alpha=2e-5,T[0]=100,L=1,
> T[A]=0}: params;
We now substitute the parameters into the list of
equations:
> eqnlist2 :=
subs(hspec,params,eqnlist):
and into the list of initial values:
> inclist2:=subs(params,inclist):
The integration is carried out using the following command;
> parta := BESIRK(eqnlist2,
> inclist2,0..6000):
The output from BESIRK is an array, the
first row of which lists the names of each variable. Subsequent rows provide
the numerical values of these variables. We can select some of the results with
the datatable function from the BESIRK package. The output shown below this command has the
same structure as the actual output, but the numbers have been rounded to 0.01
degrees Celcius. Also, the lines surrounding the numbers do not appear in the
actual Maple output.
Ø
datatable(parta,[t,seq(T[i],
> i=[1,3,5,7,9,11])]);
Ø
|
t |
T1 |
T3 |
T5 |
T7 |
T9 |
T11 |
|
0 |
100 |
100 |
100 |
100 |
100 |
100 |
|
6000 |
0 |
31.71 |
58.49 |
77.46 |
88.22 |
91.66 |
BESIRK uses the 3rd order semi-implict
Runge-Kutta method of Michelsen (1976), combined with a Bulirsch-Stoer
extrapolation technique. It is a very efficient integrator and found the
solution at t = 6000 seconds in a single large integration step. While these results are
in almost perfect agreement with the numerical results given by Cutlip and
Shacham it is not possible to obtain an accurate plot of the temperature profiles
with only two points. We repeat the computation with a limit on the maximum
step size is to ensure that the resulting time history appears to be relatively
smooth.
> parta := BESIRK(eqnlist2,inclist2,
> 0..6000,hmax=100):
We may plot the solution as a function of the independent
variable using the datatableplot procedure which is part of
the tplot package. The second
argument to this function is a list of the two variable names that are to be
plotted.
> datatableplot(parta,[[t,T[2]],
> [t,T[3]],[t,T[4]],[t,T[5]]],
> color=black,axes=boxed,
> view=[0..6000,0..110],thickness=2);
A 3D plot of the solution as a function of time using the
MOLplot procedure which is part of the BESIRK
package. Maple is capable of displaying surfaces in a variety of different
styles (here Maple has used it's default patch style).
Ø
MOLplot(parta,subs(params,
> rhs(hspec)),2..62,axes=boxed,
> contour,labels=[`t`,`x`,`T`],
> orientation=[-50,50]);
The actual time (in seconds on a 450 MHz Pentium PC) needed
to compute this solution can be found using the following command:
> BESIRKtime;
The second part of the example in Cutlip
and Shacham is to repeat the above calculations with 20 and 30 sections. To
complete this part of the assignment it suffices to re-run the above commands,
the sole change being to set the number of lines to 21 and 31 respectively. The
third part of the problem is to repeat the first part where heat convection is
present at the slab surface with h =
25 W/m2K and a thermal conductivity (in the solid) of k = 10 W/mK. The boundary condition for this situatation becomes:
> BC1 := x=0,h*(T[infinity]-T(x,t))
> = - k*diff(T(x,t),x): BC1;
The boundary condition at the surface of the slab may be
converted to finite difference form as follows:
> BC1a:=convert(BC1[2],fddiff,
> order=[2,0],forward=[2,0],
> indexletters=[1,none]): BC1a;
Proceeding as before, the parameters in the model are given
numerical values as follows
> params := {alpha=2e-5,T[0]=100,L=1,
> T[A]=0, T[infinity]=0,k=10,h=25}:
> params;
We form a list of all the equations (with the parameters
included),
> eqnlist2 := subs(hspec,params,
> [BC1a,seq(PDE1,i=2..n-1),BC2a]):
The initial condition is unchanged from
before.
> inclist2:=subs(params,inclist):
The numerical solution is obtained using
BESIRK,
> partc := BESIRK(eqnlist2,inclist2,
> 0..6000,hmax=100):
and the results displayed:
> datatableplot(partc,[[t,T[2]],
> [t,T[3]],[t,T[4]],[t,T[5]]],
> color=black,axes=boxed,
> view=[0..6000,50..110],thickness=2);
> MOLplot(partc,subs(params,
> rhs(hspec)),2..62,axes=boxed,
> style=patchcontour,
> labels=[`t`,`x`,`T`],
> orientation=[-50,50]);
Here we have
used the patchcontour style for the 3D surface plot.
The Maple worksheet for the example
described above and the code packages used here are available on the authors
web site (http://www.clarkson.edu/~chengweb/
faculty/taylor/maple/). Additional worksheets illustrating the application of
Maple to solve PDEs by the method of lines or by finite difference methods also
can be found there.
References
1.
Cutlip, M.B., and M.
Shacham The Numerical Method Of Lines for Partial Differential Equations, CACHE
News, No. 47, 18-21, 1999.
2.
Michelsen,
M.L. AIChE J., 22, 554, 1976.
3.
Rosen, E.M. Excel 7.0:
Partial Differential Equations, CACHE News, No. 47, 18-21, 1999.
4.
Schiesser, W.E. The
Numerical Method Of Lines, Academic Press, New York, 1992.
5.
Schwalbe, D., H.A.
Kooijman, and R. Taylor, Solving Stiff Differential Equations and Differential
Algebraic Systems with Maple V Maple Tech, 3(2),
47-53, 1996.
6.
Silebi,
C.E., and Schiesser, W.E., Dynamic Modeling of Transport Process Systems,
Academic Press, New York, 1992