Chapter 8: The wave equation

8.0 Many coupled oscillators

Remove["Global`*"]

n=3;
A=Table[0,{n},{n}];
A[[1,1]]=-2;
A[[1,2]]=1;
A[[n,n-1]]=1;
A[[n,n]]=-2;

Do[
A[[i,i-1]]=1;
A[[i,i]]=-2;
A[[i,i+1]]=1,
{i,2,n-1}]

A=A k;
A//MatrixForm

vecX=Table[(Subscript[x, i])[t],{i,1,n}]

eqn=Table[m D[vecX[[i]],{t,2}]==A[[i]].vecX,{i,1, n}]

L=1;
f[x_]:=Exp[-((((x-L/2))^2)/0.01)]
Plot[f[x],{x,0,1}]

init=Flatten[Table[{(Subscript[x, i])[0]==f[(i-1)(L/((n-1)))],(Subscript[x, i])'[0]==0},{i,1,n}]]

k=1;
m=1;
T=80;
sol=NDSolve[Flatten[{eqn,init}],vecX,{t,0,T}][[1]];

Xsol=Table[(Subscript[X, i])[t_],{i,1,n}];

Xsol=vecX/.sol;

tab[t_]=Table[{Xsol[[i]]+i-1,1},{i,1,n}];

Manipulate[ListPlot[{tab[t]},PlotRange->{{-1.5,L*n+1.5},{0,2}}],{t,0,T}]

tabv[t_]=Table[{i-1,Xsol[[i]]},{i,1,n}];

Manipulate[ListPlot[{tabv[t]},PlotRange->{{-1.5,L*n+1.5},{-1.5,1.5}}],{t,0,T}]

8.1 One-dimensional wave equation without damping

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 μ=0.5 kg/m solidly fixed at both ends (like a guitar string). The string is initially displaced from rest according to
	η(x,t=0)=A (e^(-(((((x-L/2)/λ)))^2))),
where A=0.1 m and λ=3 m.
What is η(x,t) for the first 25 seconds?

Remove["Global`*"]

(* Parameters *)
Tension=2; 
μ=0.5;
L=20.0;
A=0.1;
λ=3;

(* Simulation time *)
T=25.;

(* Phase speed *)
c=(Sqrt[(Tension/μ)]);

Let's examine the intial condition,

Plot[A (E^(-(((((x-(L/2))/λ)))^2))),{x,0,L},AxesLabel->{"x (m)","η (m)"},PlotRange->All]

Without solving the problem, can you guess qualitatively what will happen to this initial disturbance?

Here's how we can solve this problem numerically using NDSolve.

u[x_,t_]=η[x,t]/.Flatten[NDSolve[{(∂_(t,t))η[x,t]==(c^2)(∂_(x,x))η[x,t],η[x,0]==A (E^(-(((((x-L/2)/λ)))^2))),Derivative[0,1][η][x,0]==0,η[0,t]==0,η[L,t]==0},η[x,t],{x,0,L},{t,0,T}]]

And here's how the solution looks like with a 3D plot

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

or, with a contour plot which you may find clearer,

ContourPlot[u[x,t],{x,0,L},{t,0,T},PlotPoints->100,FrameLabel->{"x (m)","t(s)"},ColorFunction->Hue]

Make sure you understand this graph.

It's even more fun to see an animation of the solution

Manipulate[Plot[u[x,t],{x,0,L},PlotRange->{{0,L},{-1.1A,1.1A}}],{t,0,T}]

If the righward boundary is soft then the boundary condition will be ((∂ η)/(∂ x))(x=L,t)=0, i.e.

u[x_,t_]=η[x,t]/.Flatten[NDSolve[{(∂_(t,t))η[x,t]==(c^2)(∂_(x,x))η[x,t],η[x,0]==A (E^(-(((((x-L/2)/λ)))^2))),Derivative[0,1][η][x,0]==0,η[0,t]==0,Derivative[1,0][η][L,t]==0},η[x,t],{x,0,L},{t,0,T}]]

ContourPlot[u[x,t],{x,0,L},{t,0,T},PlotPoints->100,FrameLabel->{"x (m)","t(s)"},ColorFunction->Hue]

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 explicit centered difference schemes.  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 gridf 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, let simply mention that if everythinh is correctly coded but that the solution does not make sense it is possible that this is beacsue the time step chosent is too large. In this course, you will always be told what time step and grid size to use for a given problem. 
Ok, let's now code Equation 2.

(* Grid size and timestep *)
Δx=0.5;
Δt=0.05;
nmax=Round[T/Δt];
imax=Round[L/Δx]+1;

(* Hard boundary conditions on both ends *)
η0=0;

Now, let's define the grid:

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

and initial conditions:

t[1]=0;
Do[η[i,n]=A (E^(-(((((x[i]-L/2)/λ)))^2))),{i,1,imax},{n,0,1}];

Note how the initial displacement is defined not only at n=1 but also for n=0. This is because initially we need to set ∂η/∂t=0 at t=0.

In this example, the boundary conditions can be imposed at the outset.

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

We can then carry out the simulation once we have discretized equation (1) as shown in class. Here's how this could be done: (Note this loop can take a few tens of seconds to execute)

Do[{t[n+1]=n Δt,
η[i,n+1]=((Δ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,T}},FrameLabel->{"x (m)","t(s)"},ColorFunction->Hue]

8.2 One-dimensional wave equation with linear damping

With damping the 1D wave equation is:

(((∂^2)η)/(∂(t^2)))+γ((∂η)/(∂t))=(c^2)(((∂^2)η)/(∂(x^2)))    (3)

Remove["Global`*"]

(* Parameters *)
Tension=2; 
μ=0.5;
L=20.0;
γ=0.2;
A=0.1;
λ=3;

(* Simulation time *)
T=25.;

(* Phase speed *)
c=(Sqrt[(Tension/μ)]);

u[x_,t_]=η[x,t]/.Flatten[NDSolve[{(∂_(t,t))η[x,t]+γ (∂_t)η[x,t]==(c^2)(∂_(x,x))η[x,t],η[x,0]==A (E^(-(((((x-L/2)/λ)))^2))),Derivative[0,1][η][x,0]==0,η[0,t]==0,η[L,t]==0},η[x,t],{x,0,L},{t,0,T}]]

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

ContourPlot[u[x,t],{x,0,L},{t,0,T},PlotPoints->100,FrameLabel->{"x (m)","t(s)"},ColorFunction->Hue]

Let's now code our own algorithm. Equation (3) can be discretized in two different ways. The first is to use a backward difference scheme for the ∂η/∂t term. The second method is to use a centered scheme for ∂η/∂t term. The second method is more accurate. 
Here's discrtetized version of (3) using the first method where the term ∂η/∂t is discetized with a backward scheeme:

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

The we can solve for (Subscript[η, i]^(n+1)), i.e.,

(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))-γΔt((Subscript[η, i]^n)-(Subscript[η, i]^(n-1)))  (4)

Here's discrtetized version of (3) using the second, more accurate, method where the term ∂η/∂t is discetized with a centered scheme:

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

The we can solve for (Subscript[η, i]^(n+1)), i.e.,

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

(* Grid size and timestep *)
Δx=0.5;
Δt=0.05;
nmax=Round[T/Δt];
imax=Round[L/Δx]+1;

(* Hard boundary conditions on both ends *)
η0=0;

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

t[1]=0;
Do[η[i,n]=A (E^(-(((((x[i]-L/2)/λ)))^2))),{i,1,imax},{n,0,1}];

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

Nsol1=Table[η[i,n],{n,1,nmax},{i,1,imax}];

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

Nsol2=Table[η[i,n],{n,1,nmax},{i,1,imax}];

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