Appendix A3. Partial Differential Equations

The wave equation

The following entries show how Mathematica 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^2u_ (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) = Underoverscript[∑, n = 1, arg3] (A_ncos((c n π t)/L) + B_nsin((c n π t)/L)) sin((n π x)/L) ,  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) = Underoverscript[∑, n = 1, arg3] A_n sin((n π x)/L)

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

∂_t U_N(x, 0) = Underoverscript[∑, n = 1, arg3] (c n π B_n)/L sin((n π x)/L)

approximates g(x). Consequently, A_n is the Fourier sine series coefficient for f (x) and  (c n π 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) = 0.2 x + (0.2 (1 - x) - 0.2 x) UnitStep(x - 0.5)

See the following definitions and plot.

In[35]:=

L = 1; c = 1;
f[x_] := 0.2*x + (0.2*(1-x)-0.2*x)*UnitStep[x-0.5];
Plot[ f[x], {x,0,1}, PlotRange->{0,0.125}, AspectRatio->1/3 ]

[Graphics:HTMLFiles/Math_A3_12.gif]

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

g(x) = 0.1 DirecDelta(x - 0.25)

In[34]:=

g[x_] := 0.1*DiracDelta[x - 0.25]

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

The next entries calculate the formulas for the coefficients An and Bn.

In[39]:=

An = 2/L*Integrate[f[x]*Sin[n*Pi*x/L], {x,0,L} ]
Bn = L/(c*n*Pi)*Integrate[g[x]*Sin[n*Pi*x/L], {x,0,L} ]

Out[39]=

2 ((0. Cos[1.5708 n])/n + (0. Cos[n π])/n + (0.0405285 Sin[1.5708 n])/n^2 - (0.0202642 Sin[n π])/n^2)

Out[40]=

(0.031831 Sin[0.785398 n])/n

This is the definition of the function U as a function on N, x, t :

In[43]:=

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}]

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).

In[50]:=

Plot[ {g[x], Sum[ c*n*Pi*Bn/L*Sin[n*Pi*x/L], {n,1,30}]}, {x,0,L}, PlotRange->{-0.5,2}, PlotLabel->"Initial Velocity"]

[Graphics:HTMLFiles/Math_A3_16.gif]

This curve is a typical approximation to a Dirac delta. The area under the curve is approximately 1/10.
The fjollowing plot of U(20,x,0) shows that the An coefficients are also correct.

In[54]:=

Plot[ U[20,x,0], {x,0,L}, PlotRange->{0,0.12}, PlotLabel->"Initial Waveform", AspectRatio->1/3]

[Graphics:HTMLFiles/Math_A3_17.gif]

A snapshot of the waveform at t = 0.2.

In[77]:=

Plot[ U[50,x,0.2], {x,0,L}, AspectRatio->1/3]

[Graphics:HTMLFiles/Math_A3_18.gif]

Out[77]=

⁃Graphics⁃

Five snapshots, one every 0.2 seconds:

In[59]:=

Plot[ Evaluate[Table[U[50,x,t], {t,0,1,0.2}]], {x,0,L}]

[Graphics:HTMLFiles/Math_A3_20.gif]

A movie (see the Help Browser: A Practical Introduction to Mathematica/Graphics and Sound/Special Topic: Animated Graphics).

In[68]:=

<<Graphics`Animation`
Animate[Plot[U[50,x,t], {x,0,L}, PlotRange->{-0.1,0.1}], {t,0,2,0.05}]

[Graphics:HTMLFiles/Math_A3_62.gif]

The waveform surface

In[76]:=

Plot3D[ U[50,x,t], {x,0,L}, {t,0,2}, ViewPoint->{1.542, -2.913, 0.764} ]

[Graphics:HTMLFiles/Math_A3_63.gif]


Created by Mathematica  (December 9, 2004)