L-34 MCS 494 Wednesday 13 November 2002

> restart;

Partial Differential Equations in Maple

The command to solve PDEs in Maple is pdsolve . Maple offers some tools to handle PDEs in the package PDEtools :

> with(PDEtools);

In this worksheet we illustrate the command pdsolve and some of these tools.

Plotting a wave

The command PDEplot from PDEtools can handle first-order PDEs.

The first PDE we encountered provided a simplistic model for air quality. It occurs in our textbook on page 218 and has applications to traffic flow.

> wave := diff(u(x,t),t) + diff(u(x,t),x) = 0;

At the beginnin of time, we consider the wave 1-x^4 over the interval [-1,1]. We plot for (x,t) in [-1,3]x[0,2].

> PDEplot(wave,[s,0,1-s^4],s=-1..1,x=-1..3,t=0..2,animate=true);

Notice that in the command above we invoked the animate flag. If we click on the picture and then on the animation toolbar above to start the animation we see the initial wave flowing over the surface. The following two instructions create an animation which shows the build up of the surface:

> frames := seq(PDEplot(wave,[s,0,1-s^4],s=-1..1,x=-1..3,t=0..k/20),k=1..20):

> plots[display](frames,insequence=true);

Unfortunately, this PDEplot only works for first-order PDEs and not for second-order PDEs like the heat equation.

Solving the heat equation

When calling pdsolve on a PDE, Maple attempts to separate the variables.

Consider the heat equation, to model the change of temperature in a rod.

> heat := diff(u(x,t),t) = diff(u(x,t),x$2);

The Maple command to solve PDEs is pdsolve:

> ans := pdsolve(heat,u(x,t));

We see that Maple proposes the separation of variables and presents the decoupled ODEs. It is not so hard to select the ODEs:

> ode1 := op([2,1,1],ans);

> ode2 := op([2,1,2],ans);

We can now solve the ODEs:

> f1 := op(2,rhs(ode1));

> sol1 := dsolve(ode1,f1);

> f2 := op(2,rhs(ode2));

> sol2 := dsolve(ode2,f2);

> sol := rhs(sol1)*rhs(sol2);

Let us use the boundary and initial conditions:

> conds := u(0,t) = 0, u(1,t) = 0, u(x,0) = 1;

> pde := [heat,conds];

> pdsolve(pde,u(x,t));

We recognize the solution we derived in class, but apparently Maple is unable to take advantage of the conditions to arrive at the series solution we obtained on Monday.

Let us try to apply the fourier transforms, abbreviated as ft and ift:

> alias(ft = inttrans[fourier]):

> alias(ift = inttrans[fourier]):

First we apply the Fourier transform with respect to the space coordinate:

> ft_heat_x := ft(heat,x,omega);

> ode := subs(fourier(u(x,t),x,omega) = U(t), ft_heat_x);

> sol_t := dsolve(ode, U(t));

We must assume time is positive, otherwise Maple cannot determine whether the integral converges or not:

> assume(t>0);

> uf := ift(rhs(sol_t),omega,x);

> u1 := subs(_C1=1,uf);

> plot3d(u1,x=0..1,t=0..20,axes=boxed);

Note that this is a solution of the heat equation as an initial value problem. We have to separate variables to solve the equation with boundary conditions.

> t := 't': # remove the assumption on t

Change of Coordinates

With dchange of PDEtools we can change coordinates easily:

> LaplacePDE := diff(F(x,y),x$2) + diff(F(x,y),y$2) = 0;

Let us change to polar coordinates:

> polar_PDE := dchange({x=r*cos(phi),y=r*sin(phi)},LaplacePDE);

taking the numerator simplifies the equation a lot:

> numer(lhs(polar_PDE)) = 0;

We encounter the Laplacian when studying the vibrations of a circular drum (see Problem 5 on page 204 of the textbook):

> drum := subs(F(phi,r) = u(r,phi,t),numer(lhs(polar_PDE)));

Just as in the textbook, we are going to assume that we strike the drum at the very center so that the vibration is radially symmetric and we can ignore the angle phi.

> drum := subs(diff(u(r,phi,t),phi$2) = 0,drum);

> drum := normal(drum/r) = diff(u(r,phi,t),t);

> ans := pdsolve(drum,u(r,phi,t));

The equation in r is interesting as it leads to Bessel's equation:

> ode_r := op([2,1,3],ans);

Careful here: if you execute the above command, the interesting ODE in r may be at a different location...

> n_ode_r := normal(r^2*ode_r);

We need to suggest the right type of constant to arrive at the Bessel equation:

> bessel := subs(_c[1] = -1,_F1(r)=y(r),n_ode_r);

> b_sol := dsolve(bessel,y(r));

> j0 := rhs(subs(_C1 = 0,_C2 = 1,b_sol));

> plot(j0,r=0..20);

So we can think of the solution to the Bessel equation as a damped cosine.

>