Maple_A3.mw

Appendix A3. Partial Differential Equations

The wave equation

The following entries show how Maple can be used to plot approximations to solutions of the wave equation on a finite domain: A string of length L with ends clamped at x = 0 and x = L. Let u(x,t) denote the vertical displacement of the string at point x at time t. For small vibrations u satisfies the wave equation

u[t, t] = c^2*u[x, x] .

The letter c denotes a positive constant determined by the characteristics of the string. Separation of variables leads to solutions of the following form

U[N](x, t) = sum((A[n]*cos(c*n*Pi*t/L)+B[n]*sin(c*n*Pi*t/L))*sin(n*Pi*x/L), n = 1 .. N) , N a positive integer.

See Ledder, Chapter 8, Section 3.

Set the string into motion

The string is set into motion at t = 0 by giving it an initial shape  f (x)  and an initial velocity distribution, g(x). Thus the coefficients A[n] and B[n] should be chosen so that the function

U[N](x, 0) = sum(A[n]*sin(n*Pi*x/L), n = 1 .. N)

approximates f (x) on [0 , L ] and the function

eval(diff(U[N](x, t), t), t = 0) = sum(c*n*Pi*B[n]/L, n = 1 .. N) sin(n*Pi*x/L)

approximates g(x). Consequently, A[n] is the Fourier sine series coefficient for f (x) and c*n*Pi*B[n]/L is the Fourier sine series coefficient for g(x).

The following entries define the functions f and g, calculate A[n] and B[n] , then create various solution curves. We assume that L = 1, c = 1 and the string is initially stretched "tent like" over the  x axis  with the shape

f(x) = piecewise(x < Float(5, -1), Float(2, -1)*x, Float(2, -1)*(1-x))

Set it into motion with a finger flick at a point one quarter of the way from the left endpoint

g(t) = Float(1, -1)*delta(t-Float(25, -2))

You may, of course, change these to fit any situation that you would like to explore.

> L := 1: c := 1:
f := x -> piecewise(x<L/2,2/5*x,2/5*(L-x)):

g := x -> 1/10*Dirac(x - L/4):

An := 2/L*int(f(x)*sin(n*Pi*x/L), x=0..L):

Bn := L/(c*n*Pi)*2/L*int(g(x)*sin(n*Pi*x/L), x=0..L):

The following entry simplifies the formulas for An and Bn, then displays them.

> C := [An,Bn] assuming n::integer: 'An'=C[1],'Bn'=C[2];

An = -2/5*(-2*sin(1/2*n*Pi)+cos(1/2*n*Pi)*n*Pi)/(n^2*Pi^2)+2/5*(cos(1/2*n*Pi)*n*Pi+2*sin(1/2*n*Pi))/(n^2*Pi^2), Bn = 1/5*sin(1/4*n*Pi)/(n*Pi)

This it the definition of U as a function of N, x, and t.

> U := (N,x,t) -> sum((An*cos(c*n*Pi*t/L)+Bn*sin(c*n*Pi*t/L))*sin(n*Pi*x/L),
                    n=1..N);

U := proc (N, x, t) options operator, arrow; sum((An*cos(c*n*Pi*t/L)+Bn*sin(c*n*Pi*t/L))*sin(n*Pi*x/L), n = 1 .. N) end proc

The first plot checks that the coefficients are correct for the velocily function g. (A check for the shape function f  is made when we plot U at t  =0 below).

> plot( [g(x), sum(c*n*Pi*Bn/L*sin(n*Pi*x/L), n=1..30)], x=0..L,
      color=[red,blue], title="Initial Velocity");

[Plot]

This curve is a typical approximation to a Dirac delta. The area under the curve is approximately 1/10.

The plot of U(20,x,0) shows that the An coefficients are also correct.

> plot( U(20,x,0), x=0..L, 0..0.25, ytickmarks=3,
     title="Initial Waveform");

[Plot]

A snapshot of the waveform at t = 0.2:

> plot( U(50,x,0.2), x=0..L, color=red);

[Plot]

Five snapshots, one every 0.2 seconds:

> plot( [U(50,x,0.2*t)$t=0..5], x=0..L, color=red);

[Plot]

A movie (see the Help page for plots[animate]):

> plots[animate]( plot, [ U(50,x,t), x=0..L ], t=0..2, frames=40);

[Plot]

The waveform surface

> plot3d( U(50,x,t), x=0..L, t=0..2, axes=boxed, orientation=[-60,70],
       lightmodel=light2 );

[Plot]