The Wave Equation - Find the Errors!

Name:
Student number:

As seen in class, the 1D wave equation without damping is:
	(((∂^2) η)/(∂ (t^2)))=(c^2)(((∂^2) η)/(∂ (x^2))) (1)
where η is the displacement from equilibrium, c is the wave phase speed, t is time and x is space. This is a partial differiential equation because η depends on both time and space. 
For waves propagating on a string the phase speed depends on the string itself (its mass) and the tension in the string:
	c=(Sqrt[(T/μ)])
where T is the string tension (in N) and μ is the mass per unit length of the string (in kg/m).
Consider a string of length L=20 m, T=2N and μ=3.125 kg/m solidly fixed at both ends (like a guitar string).

You will now go over solving η(x,t) for the first 5 seconds, first using NDSolve, then writing a finite differnece routine.  The solutions presented here are peppered with errors.  You will have to find them and correct them.

Remove("Global`*")

First, let's store the important parameters of the system.  Check these carefully.

L=20;
T=2;
μ=3.125;
c=Sqrt[(T/μ)]
tf=50;

A partial differential equation like the wave equation requires initial conditions.  Here is the initial displacement of the string, i.e. a fairly sharp pluck at x=9m.

b[x_]=If[x>17 || x<1,0,0.5 Exp[ -(((x-9))^2)/3]-Exp[-(64/3)]]; (* This is fine *)

Plot[b[x],{t,0,L},PlotRange->{{0,20},{0,0.7}}]

Additionally, we need to specifiy the initial string velocity.  For this problem, it will be zero everywhere.

We have already used NDSolve for ODE's.  It can also be used to solve PDE's, as we will do now.

First, let's enter the wave equation itself,

eq=(∂_(t,t))η[x,t]==(c^2) (∂_x)η[x,t];

Second, let us specify the initial diplacement

init1=η[x,0]=b[x];

and initial velocity (we must use Derivative to specify conditions on the derivatives of η).

init2=Derivative[0,1][η][x,0]==0;

As we are dealing with a PDE, we must specify boundary conditions.  Here, they are simple: at all times, the displacements at x=0 and x=L are zero

bc1=η[t,0]==0;

bc2=η[t,L]==0;

The syntax for using NDSolve for PDE's is similar to using it for ODE's.  The only big difference is that the functions we solve for are functions of more than one variable.  Schematically, we have
NDSolve[{equation plus initial conditions and boundary conditions},{function we are solving for},{x domain},{t domain}]

sol=NDSolve[{eq,initl,init2,bc1,bc2},η[x,t],{x,0,L},{t,0,T}]//First

s[x_,t]=η[x,t]/.sol

Plot the displacement at t=2s.

Plot[s[x,2],{x,0,tf},PlotRange->{{0,L},{-0.7,0.7}}]

Now again at t=4s.

Plot[s[x,4],{x,0,L},PlotRange->{{0,L},{-0.7,0.7}}]

We can see all the times at once with a 3D plot.  One axis is x, while the other is time.  The height gives the diplacement.

Plot3D[η[x,t],{x,0,L},{t,0,tf},PlotPoints->100,AxesLabel->{"x (m)","t(s)","η(m)"}]

Alternatively, this plot can be viewed as a contour plot.  The different colours represent different elevations.

ContourPlot[η[x,t],{x,0,L},{t,0,tf},AxesLabel->{"x (m)","t(s)","η(m)"},ColorFunction->Hue]

Perhaps most fun and informative is viewing the solution as an animation.

Manipulate[Plot[s[x,i],{x,0,L},PlotRange->{{0,L},{-0.7,0.7}}],{i,0,tf,0.25}]

Let's now solve this problem using our own algorithm instead of relying on NDsolve. As shown in class we first need to discretize the equation. Let's use an explicit centered difference scheme.  Then, the only unknown term should be (η^(n+1)). Discretized Equation (1) will look like:

(((Subscript[η, i]^(n+1))-2(Subscript[η, i]^n)+(Subscript[η, i]^(n-1)))/((Δt)^2))=(c^2)(((Subscript[η, i+1]^n)-2(Subscript[η, i]^n)+(Subscript[η, i-1]^n))/((Δx)^2))

where n is the time index and i the spatial index. The only unknow in this equation is (Subscript[η, i]^(n+1)), so we can re-arrange it like this:

(Subscript[η, i]^(n+1))=((Δt)^2)(c^2)(((Subscript[η, i+1]^n)-2(Subscript[η, i]^n)+(Subscript[η, i-1]^n))/((Δx)^2))+2(Subscript[η, i]^n)-(Subscript[η, i]^(n-1))  (2)

The time step Δt and the grid size Δx have to be chosen. For this scheme if the time step is chosen to be too large numerical instabilities may result. Without details, we simply mention that if everything is correctly coded but that the solution does not make sense it is possible that this is beacsue the time step chosen is too large.  In this course, you will always be told what time step and grid size to use for a given problem.

Δx=0.1;
Δt=0.05;

nmax=Round[T/Δt]
imax=Round[L/Δx]+1

nmax is the number of time steps we will take, while imax is the number of grid points we are using to define our x domain.

Now, let's define the grid:

Do[x[i]=(i-1)Δx,{i,1,imax}];

x[1] should be zero, while x[imax] should be L.

We will only use x[i] to help us define our initial displacement, which follows

Do[η[i,1]=b[x],{i,1,imax}];

How do we specify that the initial velocity is zero, i.e. that ∂η/∂t=0 at t=0?  If the time derivative at a given point on the x domain is zero, that means the function doesn't change from one time to the next, approximately speaking.  So to accomplish ∂η/∂t=0 at t=0, we will define η at t=0-Δt to be the same as at t=0.

Do[η[i,0]=b[x],{i,1,imax}];

What about boundary conditions?  In this case, the displacement is fixed at zero at both ends (at x=0 and x=L) for all time steps, so we write

Do[η[1,n]=0,{n,1,nmax+1}]

Do[η[imax-1,n]=0,{n,1,nmax+1}]

We can now iterate the discretized equation (2). Here's how this could be done, noting that the Do loop executes the outer index first, i.e. i changes from 2 to imax-1 before n goes up by 1.  Note this loop can take a few tens of seconds to execute.

Do[
η[i,n]=((Δt)^2) (c^2) (((η[i+1,n]-2η[i,n]-η[i-1,n])/((Δx)^2)))+2 η[i,n]-η[i,n-1],
{n,1,nmax},{i,2,imax-1}]

(* Create a table from the solution *)
Nsol=Table[η[i,n],{n,1,nmax+1},{i,1,imax}];

We can then plot the solution with the following command:

ListContourPlot[Nsol,DataRange->{{0,L},{0,tf}},FrameLabel->{"x (m)","t(s)"},ColorFunction->Hue]