Wave Equation Finite Differences

Centered Difference Scheme

Run the following to clear all variables:

Clear["Global`*"];

Parameters

Length of the string is len, wave speed is c, time step is Δt, number of spatial intervals is n. Assign some values to these parameters here:

len=1;
c=0.5; 
Δt = .01; 
n=100;

From the previous parameters we can compute Δx and σ:

Δx=(len/n);
σ=((c Δt)/(Δx));

Inspect the value of Δx:

N[Δx]

We will also specify the time length of our simulation here, so it’s easy to change later.

maxTime=4;
jmax=Ceiling[(maxTime/(Δt))]

Initial conditions

Let the initial position be f(x)=0.4(e^(-300(((x-0.5))^2)))+0.1(e^(-300(((x-0.65))^2))), and let the initial velocity be g(x)=0.

f[x_]:=0.4Exp[-300(((x-0.5))^2)]+0.1Exp[-300(((x-0.65))^2)];
g[x_]:=0;
Plot[{f[x],g[x]},{x,0,len},PlotRange->All]

We will evaluate f(x) at intervals of Δx and store the values in u[0], with the zero denoting time zero.

u[0]=N[Table[ f[m*Δx],{m,1,n-1}]];
ListPlot[u[0],DataRange->{0,1},PlotRange->All,AxesLabel->{"x","displacement"}]

Boundary conditions

For now, we will use zero Dirichlet boundary conditions (that is, endpoints of the string fixed at zero). We don’t have to do anything with these boundary conditions.

Coefficient matrix

The centered difference scheme uses the following finite difference equation:

(Subscript[u, j+1,m])=(σ^2)(Subscript[u, j,m-1])+2(1-(σ^2))(Subscript[u, j,m])+(σ^2)(Subscript[u, j,m+1])-(Subscript[u, j-1,m])

We build the tridiagonal matrix that encodes the coefficients of the solution values at time step j. To do this, first create an (n-1)×(n-1) zero matrix, and then use Do loops to set values along three diagonals. Replace the question marks below with the correct coefficients.

B = ConstantArray[0, {n-1,n-1}];
Do[B[[j,j]] = ????,{j,1,n-1}]; (* set entries on the main diagonal *)
Do[B[[j,j+1]] = ????,{j,1, n-2}]; (* set entries just above the main diagonal *)
Do[B[[j+1,j]] =????,{j,1,n-2}]; (* set entries just below the main diagonal *)
MatrixForm[B];

Compute the approximate solution

First we must find a second row of initial values so that we can launch the iteration. We’ll use the method derived in class, using the centered difference approximation of ((∂u)/(∂t)).

u[1]=N[Table[((σ^2)/2) f[(m-1)Δx]+(1-(σ^2))f[m Δx]+((σ^2)/2) f[(m+1)Δx]+g[m Δx]Δt,{m,1,n-1}]];

Now we can iterate and find the numerical solution from time index j=2 to j=jmax. We store the solution values at time step (Subscript[t, j]) in u[ j ]. We use the equation (u^((j+1)))=B (u^((j)))-(u^((j-1))). The following loop does all the work:

Do[u[j+1]=B.u[j]-u[j-1],{j,1,jmax}]

We can append our vector u to include zero displacement at the boundaries x = 0 and x = 1:

Do[PrependTo[u[j],0],{j,-1,jmax}]
Do[AppendTo[u[j],0],{j,-1,jmax}]

Inspect some values of our approximate solution

Examine the solution at a particular time step, such as m=10, which corresponds to t=0.1:

u[2]

Plot the approximate solution

Since u[j] is the form of a list, we can use ListPlot to plot the solution at a particular time.

ListPlot[u[40],Joined->True,DataRange->{0,len},PlotRange->All]

Now make a bunch of plots and animate them.

plots = Table[ListPlot[u[j],Joined->True,PlotRange->{-.5,.5},DataRange->{0,len}],{j,0,jmax}];
ListAnimate[plots]

Questions

Try some different initial conditions. Do the solutions appear reasonable?

Try running the numerical scheme with different choices of the wave speed. For what wave speeds is this scheme stable?

How would you incorporate nonhomogeneous time-dependent boundary conditions? For example: u(t,0)=a(t) and u(t,len)=b(t) for some functions a(t) and b(t)? Try it!