Transport Equation Finite Differences

MATH 330 Differential Equations II

Forward Difference Approximations

Run the following to clear all variables:

Clear["Global`*"]

Parameters

Let the length in the x-direction be len, the number of spatial intervals be n, the transport speed be c, and time step be Δt.
Assign some values to these parameters here:

len=1;
n=200;
c=0.5; 
Δt = .005;

From the previous parameters we can compute the spatial subinterval length Δx and the parameter σ=((c Δt)/(Δx)):

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

Inspect the value of Δx:

N[Δx]

Initial condition

Let the initial wave form be f(x)=0.4(e^(-300(((x-0.5))^2)))+0.1(e^(-300(((x-0.65))^2))). Define this and make a plot:

f[x_]:=0.4Exp[-300(((x-0.5))^2)]+0.1Exp[-300(((x-0.65))^2)];
Plot[f[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. This will give the initial values for our approximate solution.

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

Matrices

We will compute the approximation at each time step from the previous time step via a matrix equation of the form (u^((j+1)))=A (u^((j))).

To build matrix A, we create a (n+1)×(n+1) zero matrix, and then use Do loops to set values along two diagonals. Fill in the missing values in the code below.

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

To see matrix A, remove the semicolon in the last line above.

Compute the approximate solution

Next we can iterate and find the numerical solution.  Let's go from m=0 to m=100. We store the solution values at time step j in u[ j ]. As mentioned before, we use the equation (u^((j)))=A (u^((j-1))). The following loop does all the work:

Do[u[j+1]=A.u[j],{j,0,99}]

Plot the approximate solution

Since u[m] returns a list , we can use ListPlot to plot the solution at a particular time.  The option Joined draws a solid line instead of plotting discrete points, and the option DataRange scales the x-axis properly. We will make a bunch of plots and animate them.

plots = Table[ListPlot[u[j],Joined->True,DataRange->{0,1},AxesLabel->{"x","displacement"},PlotRange->{0,0.5}],{j,0,100}];
ListAnimate[plots]

Questions

What do you observe in the approximate solution to the transport equation?

Experiment with wave speeds. For example, try c=-0.5, c=-1.0, c=-1.5, and other values of your choice. What do you observe?

For what wave speeds is this approximation scheme stable? Why do you think this is?