Wave Equation: finite difference approximation via the Centered Difference Scheme

Run the following to clear all variables:

Clear["Global`*"];

Parameters

Length of the string is L, coefficient of the wave equation is (c^2), time step is Δt, number of spatial intervals is n. We assign values to these parameters here:

L=1;
c=9;
Δt = .01;
n=10;
Δx=L/n;
s=c Δt/Δx

Initial conditions

Let the initial position be f(x)=x(x-1)(x-(1/2)) and the initial velocity be g(t)=0.

f[x_]=x(x-1)(x-1/2);
g[x_]=0;
Plot[f[x],{x,0,L}]

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

U[0]=Table[ f[j*Δx],{j,1,n-1}]

Similarly we will evaluate g(x) at intervals of Δx and store the values in V[0]:

V[0]=Table[ g[j*Δx],{j,1,n-1}]

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

Matrices

The centered difference scheme uses the following finite difference equation:

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

We build the tridiagonal matrix that encodes the coefficients of the solution values at time step m. To do this, first create an (n - 1) by (n - 1) zero matrix, and then use Do loops to set values along three diagonals:

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

Compute the approximate solution

First we must find a second row of initial values so that we can launch the iteration. We’ll use option (b) from the worksheet, computing “phantom” values at time step -1. Here is the formula for these values, derived in the worksheet:

phantom[j_]:=(s^2 f[(j-1)Δx]+2(1-s^2)f[j Δx] + s^2 f[(j+1)Δx]-2Δt g[jΔx])/2;

Compute and store the values in U[-1]:

U[-1]=Table[phantom[j],{j,1,n-1}]

Now we can iterate and find the numerical solution.  Let's go from m=0 to m=100.  Since Δt=.01, this corresponds to t = 0 to t = 1. We store the solution values at time step m in U[m]. We use the equation (u^((m+1)))=A (u^((m)))-(u^((m-1))). The following loop does all the work:

Do[U[m+1]=A.U[m]-U[m-1],{m,0,100}]

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[10]

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

Do[PrependTo[U[m],0],{m,-1,100}]
Do[AppendTo[U[m],0],{m,-1,100}]

U[10]

Plot the approximate solution

Since U is the form of a list, we can use ListPlot to plot the solution at a particular time.

ListPlot[U[0],Joined->True,DataRange->{0,L}]

Now make a bunch of plots and animate them.

plots = Table[ListPlot[U[m],Joined->True,PlotRange->{-.05,.05},DataRange->{0,L}],{m,0,100}];
ListAnimate[plots]

Questions

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

Try running the numerical scheme with different choices of Δt. Does the scheme always appear to be stable?

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