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
The letter c denotes a positive constant determined by the characteristics of the string. Separation of variables leads to solutions of the following form
, 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 and
should be chosen so that the function
approximates f (x) on [0 , L ] and the function
approximates g(x). Consequently, is the Fourier sine series coefficient for f (x) and
is the Fourier sine series coefficient for g(x).
The following entries define the functions f and g, calculate and
, 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
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 ]
Set the string into motion with a finger flick at a point one quarter of the way from the left endpoint
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]=
Out[40]=
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"]
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]
A snapshot of the waveform at t = 0.2.
In[77]:=
Plot[ U[50,x,0.2], {x,0,L}, AspectRatio->1/3]
Out[77]=
Five snapshots, one every 0.2 seconds:
In[59]:=
Plot[ Evaluate[Table[U[50,x,t], {t,0,1,0.2}]], {x,0,L}]
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}]
The waveform surface
In[76]:=
Plot3D[ U[50,x,t], {x,0,L}, {t,0,2}, ViewPoint->{1.542, -2.913, 0.764} ]
Created by Mathematica (December 9, 2004)