Time Dependent Problems

Drawing Characteristic Curves - Figure 5.3

We can draw characteristics by placing straight lines with variable slope on a grid of initial values.

drawChar[v_, inx0_,xf_,dx_,Tf_] := Module[{retgraph},
retgraph = Plot[Table[(x - x0)/v[x0],{x0,inx0,xf,dx}],{x,inx0,xf},PlotRange->{{inx0,xf},{0,Tf}},PlotStyle->GrayLevel[0.5]];
Return[retgraph];
]

If we set the velocity to be a constant, we get back the parallel characteristics of the wave equation:

v[x_] := 1/2

drawChar[v,-3,3,.5,10]

But we can use more exotic velocity functions - that will lead to characteristic crossings

rho0[x_] := 20 Exp[-x^2]

rhomax = 264;
vmax = 65;

here is the velocity definition for the traffic flow model

v[x_] := vmax (1 - 2 rho0[x]/rhomax)

G0 = drawChar[v,-10,10,1,.25]

G1 = Plot[rho0[x]/rhomax,{x,-10,10}, PlotRange->{0,.25},PlotStyle->GrayLevel[0]]

we can plot the initial density, together with the curves it generates.

Show[G0,G1]

Lax Friedrichs - Implementation 5.1

This is the Lax Friedrichs solver - it takes in the initial table of spatial values u(x,0), called u0.  You must also provide a function v that computes the velocity, then give the time-step dt, number of steps to take, Nt, and spatial grid size dx (you must enforce any stability relationships between dx and dt yourself), and finally, the boundary type (a string, “periodic”, “diffuse” or “zero”).  The output is a set of tables, one for each time step, giving an approximation to u(x,t).

LFPDE[u0_, v_, dt_,Nt_,dx_, bndtype_] := Module[{retu, u,uup, Nx,i,x,j},
u = u0;
Nx = Length[u0];
uup = Table[0,{k,1,Nx}];
retu = Table[0,{k,1,Nt}];
For[i = 1, i ≤ Nt, i = i+1,
x= dx;
For[j = 2, j < Nx, j = j+1,
uup[[j]] = (1/2) (u[[j+1]] + u[[j-1]]) 
- dt/(2 dx) ( u[[j+1]] v[x+dx,u[[j+1]]] - u[[j-1]] v[x-dx,u[[j-1]]]);
x = x + dx;
];
If[bndtype=="periodic",
uup[[1]] = uup[[Nx-1]];
uup[[Nx]] = uup[[2]];
];
If[bndtype =="diffuse",
uup[[1]] = uup[[2]];
uup[[Nx]] = uup[[Nx-1]];
];
If[bndtype =="zero",
uup[[1]] = 0.0;
uup[[Nx]] = 0.0;
];
u = uup;
retu[[i]] = uup;
];
Return[retu];
]

Example 1

We’ll take the boundary points to be 0 and 1 with a periodic boundary condition, then we want our initial data to extend from 0 to 1 + dx, so that we can use periodicity to set the values of f(0) = f(1) and f(1+dx) = f(dx).

dx = .01;
u0 = Table[Sin[2 Pi x],{x,0,1+dx,dx}];
vm = 65;
dt = dx/vm;

v[x_,u_] := vm

We’ll compute the numerical solution and plot the exact solution together with it:

wave = LFPDE[u0, v, dt,100,dx,"periodic"];
exact = Table[ Table[Sin[2 Pi (x - vm t)],{x,0,1+dx,dx}],{t,dt,100 dt,dt}];
mov = Table[Show[{ListPlot[wave[[j]]],ListLinePlot[exact[[j]]]}],{j,1,Length[wave]}];
ListAnimate[mov]

Example 2 - Figure 5.5

The traffic flow example - an initially Gaussian distribution will steepen until a discontinuity forms

vm = 65;
pm = 264;
dx = 0.01;
dt = dx/vm;

Here is the velocity function appropriate to traffic flow:

v[x_,p_] := vm (1 -p/pm)

p0 = N[Table[pm Exp[-(x-5)^2],{x,0,10,.01}]];
traffic= LFPDE[p0,v,dt,300,dx,"zero"];
mov = Table[ListPlot[traffic[[j]],PlotRange->{{0.0,10/dx},{-1,270}}],{j,1,Length[traffic],3}];
ListAnimate[mov]

data = Table[Table[{(j-1) dx, traffic[[k,j]]/pm},{j,1,Length[traffic[[k]]]}],{k,1,Length[traffic]}];

Some snapshots to form Figure 5.5:

GraphicsGrid[Partition[Table[ListPlot[data[[jj]],PlotStyle->{GrayLevel[0.5], PointSize[.01]}],{jj,10,90,10}],2]]