MATH761

MATH76109/28/2018

Lab05: Finite difference schemes for hyperbolic equations

1Semi-discretization

Semi-discretization of the advection equation

q_{t} + u q_{x} = 0

using a centered difference scheme leads to the ODE system

\frac{d}{d t} 𝑸 = - \frac{u}{2 h} 𝐁 𝑸

with $𝐁 = diag ([\begin{array}{lll} - 1 & 0 & 1 \end{array}])$ . We consider the solution of the above ODE system by:

Forward Euler (FTCS scheme)
$𝑸^{n + 1} = (𝑰 - \frac{ν}{2} 𝐁) 𝑸^{n}$
with $ν = u k / h$ known as the CFL number.
Midpoint (leap-frog scheme)
$𝑸^{n + 1} = 𝑸^{n - 1} - ν 𝐁 𝑸^{n}$
Lax-Friedrichs
$𝑸^{n + 1} = (𝝁 - \frac{ν}{2} 𝐁) 𝑸^{n}$
with $𝝁 = \frac{1}{2} (diag ([\begin{array}{lll} 1 & 0 & 0 \end{array}]) + diag ([\begin{array}{lll} 0 & 0 & 1 \end{array}]))$ the averaging operator
Lax-Wendroff
$𝑸^{n + 1} = (𝑰 - \frac{ν}{2} 𝐁 + \frac{ν^{2}}{2} 𝐃) 𝑸^{n}$
with $𝐃 = diag ([\begin{array}{lll} 1 & - 2 & 1 \end{array}])$

Lab05 continues Lab04 by studying the modified equations associated with each of these schemes. The literate programming implementation of Lab04 is expanded to include analysis techniques.

2Implementation

2.1Global definitions

A module is constructed with global definitions.

SELECTED_REAL_KIND(P,R) is a Fortran intrinsic function that returns the available type closest to the requested precision of P decimal digits and an exponent range R. Two, possibly different, precisions are defined for the dependent variables ( $q$ ) and the independent variables (space, time).

2.2Time stepping

	A common interface to all finite difference schemes: method: scheme to apply m: number of interior nodes nSteps: number of time steps cfl: CFL number Q0,Q1: initial conditions Q: final state after time stepping Qlft: boundary values at left Qrgt: boundary values at right
	Internal variable declarations
	Precomputation of common expresions
	Start of SELECT statement

2.2.1FTCS

Carry out time steps.

For $u > 0$ , use specified left boundary value, extrapolate at right. For $u < 0$ use specified right boundary value, extrapolate at left.

$Q_{i}^{n + 1} = Q_{i}^{n} - \frac{ν}{2} (Q_{i + 1}^{n} - Q_{i - 1}^{n})$

Stability.

Use Von Neumann analysis

Q_{j}^{n} \sim G^{n} e^{i j θ}

, with

θ = ξ h

In[45]:=

Q[n_,j_,theta_] = G[n] Exp[I j theta]

$G (n) e^{i j θ}$

In[46]:=

FTCS = Q[n+1,j,theta] - (Q[n,j,theta] - nu/2 (Q[n,j+1,theta]-Q[n,j-1,theta]))

$\frac{1}{2} ν (G (n) e^{i (j + 1) θ} - G (n) e^{i (j - 1) θ}) + G (n) (- e^{i j θ}) + G (n + 1) e^{i j θ}$

In[47]:=

expr = Expand[FTCS/G[n]] /. G[n+1]/G[n]->A

$A e^{i j θ} - \frac{1}{2} ν e^{i (j - 1) θ} + \frac{1}{2} ν e^{i (j + 1) θ} - e^{i j θ}$

In[48]:=

AmpFact[theta_] = Simplify[ComplexExpand[A /. Solve[expr == 0,A][[1,1]]]]

$1 - i ν \sin (θ)$

In[49]:=

The amplifcation factor is $| A | ⩾ 1$ , hence the FTCS method is always unstable.

Modified equation.

In[10]:=

FTCS = s[t+k,x] - (s[t,x] - nu/2 (s[t,x+h]-s[t,x-h]))

$\frac{1}{2} ν (s (t, h + x) - s (t, x - h)) + s (k + t, x) - s (t, x)$

In[14]:=

mFTCS = Simplify[Normal[Series[FTCS,{k,0,2},{h,0,2}]]]

$h ν s^{(0, 1)} (t, x) + \frac{1}{2} k^{2} s^{(2, 0)} (t, x) + k s^{(1, 0)} (t, x)$

From above series expansions find the modified equation

s_{t} + u s_{x} = - \frac{k}{2} s_{t t} + 𝒪 (k^{2}, h^{2}) = - \frac{u^{2} k}{2} s_{x x} + 𝒪 (k^{2}, h^{2}),

an advection-diffusion equation with negative diffusivity, again indicating instability of the FTCS method.

$Q_{i}^{n + 1} = Q_{i}^{n} - ν (Q_{i}^{n} - Q_{i - 1}^{n})$ for $u > 0$

$Q_{i}^{n + 1} = Q_{i}^{n} - ν (Q_{i + 1}^{n} - Q_{i}^{n})$ for $u < 0$

Stability.

Use Von Neumann analysis

Q_{j}^{n} \sim G^{n} e^{i j θ}

, with

θ = ξ h

(assume

u > 0

)

In[87]:=

Upwind = Q[n+1,j,theta] - (Q[n,j,theta] - nu (Q[n,j,theta]-Q[n,j-1,theta]))

$ν (G (n) e^{i j θ} - G (n) e^{i (j - 1) θ}) + G (n) (- e^{i j θ}) + G (n + 1) e^{i j θ}$

In[88]:=

expr = Expand[Upwind/G[n]] /. G[n+1]/G[n]->A

$A e^{i j θ} - ν e^{i (j - 1) θ} + ν e^{i j θ} - e^{i j θ}$

In[89]:=

AmpFact[theta_] = Simplify[TrigReduce[ComplexExpand[A /. Solve[expr == 0,A][[1,1]]]]]

$- i ν \sin (θ) + ν \cos (θ) - ν + 1$

In[90]:=

A[theta_,nu_]=Sqrt[TrigReduce[ComplexExpand[AmpFact[theta] Conjugate[AmpFact[theta]]]]]

$\sqrt{- 2 ν^{2} \cos (θ) + 2 ν^{2} + 2 ν \cos (θ) - 2 ν + 1}$

In[93]:=

UpwindStabPlot=Plot[Table[A[theta,nu],{nu,0.1,1.0,0.1}],{theta,0,2Pi},

GridLines->Automatic,Axes->False,Frame->True,FrameLabel->{"theta","A(theta)"}];

Export["/home/student/courses/MATH761/lab05/UpwindStabPlot.pdf",UpwindStabPlot]

$/home/student/courses/MATH761/lab05/UpwindStabPlot.pdf$

In[94]:=

Figure 1. Upwind amplification factor obtained from Von Neumann analysis

Modified equation.

In[94]:=

Upwind = s[t+k,x] - (s[t,x] - nu (s[t,x]-s[t,x-h]))

$ν (s (t, x) - s (t, x - h)) + s (k + t, x) - s (t, x)$

In[96]:=

mUpwind = Simplify[Normal[Series[Upwind,{k,0,2},{h,0,2}]]]

$- \frac{1}{2} h^{2} ν s^{(0, 2)} (t, x) + h ν s^{(0, 1)} (t, x) + \frac{1}{2} k^{2} s^{(2, 0)} (t, x) + k s^{(1, 0)} (t, x)$

In[97]:=

From above series expansions find the modified equation

s_{t} + u s_{x} = - \frac{k}{2} s_{t t} + \frac{h^{2} ν}{2 k} s_{x x} + 𝒪 (k^{2}, h^{2}) = \frac{1}{2} (- k u^{2} + \frac{h^{2} ν}{k}) s_{x x} + 𝒪 (k^{2}, h^{2}) = \frac{ν h^{2}}{2 k} (1 - ν) s_{x x} + 𝒪 (k^{2}, h^{2}),

again an advection-diffusion equation with diffusivity

α = \frac{ν h^{2}}{2 k} (1 - ν) .

Note that the diffusivity is zero at CFL $ν = 1$ . The analytical solution to

{\begin{cases} s_{t} + u s_{x} = α s_{x x} \\ s (t = 0, x) = H (- x) \end{cases} .

s (t, x) = \frac{1}{2} [1 + \erf (\frac{x - u t}{\sqrt{4 α t}})], \erf (x) = \frac{2}{\sqrt{π}} \int_{0}^{x} e^{- ξ^{2}} d ξ

In[99]:=

s[t_,x_]=(1+Erf[(x-u t)/Sqrt[4 alpha t]])/2

$\frac{1}{2} (erf (\frac{x - t u}{2 \sqrt{α t}}) + 1)$

In[100]:=

Simplify[ D[s[t,x],t] + u D[s[t,x],x] - alpha D[s[t,x],{x,2}] ]

$0$

In[101]:=

$Q_{i}^{n + 1} = Q_{i}^{n - 1} - ν (Q_{i + 1}^{n} - Q_{i - 1}^{n})$

$Q_{i}^{n + 1} = Q_{i}^{n - 1} - \frac{ν}{2} (Q_{i + 1}^{n} - Q_{i - 1}^{n}) + \frac{ν^{2}}{2} (Q_{i + 1}^{n} - 2 Q_{i}^{n} + Q_{i - 1}^{n})$

3Numerical experiments

Compilation of the above implementation leads to a Python-loadable module than can be used for numerical experiments.

Python]

from pylab import *

Python]

import os,sys

os.chdir('/home/student/courses/MATH761/lab05')

cwd=os.getcwd()

sys.path.append(cwd)

Python]

from lab05 import *

Python]

print scheme.__doc__

q = scheme(method,cfl,q0,q1,qlft,qrgt,[m,nsteps])


Wrapper for ‘‘scheme‘‘.

Parameters
----------
method : input int
cfl : input float
q0 : in/output rank-1 array('d') with bounds (m + 2)
q1 : in/output rank-1 array('d') with bounds (m + 2)
qlft : input rank-1 array('d') with bounds (nsteps + 1)
qrgt : input rank-1 array('d') with bounds (nsteps + 1)

Other Parameters
----------------
m : input int, optional
    Default: (len(q0)-2)
nsteps : input int, optional
    Default: (len(qlft)-1)

Returns
-------

q : rank-1 array('d') with bounds (m + 2)

Python]

3.1Discontinuous boundary condition

Study now the behavior for a discontinuous (shock) initial condition $q (0, x) = H (- x)$ . The upwind solution will be compared to the exact solution $q (t, x) = H (u t - x)$ for the advection equation $q_{t} + u q_{x} = 0$ , and the exact solution

s (t, x) = \frac{1}{2} [1 + \erf (\frac{u t - x}{\sqrt{4 α t}})], \erf (x) = \frac{2}{\sqrt{π}} \int_{0}^{x} e^{- ξ^{2}} d ξ

to the modified equation $s_{t} + u s_{x} = α s_{x x}$ , with artificial diffusivity

α = \frac{ν h^{2}}{2 k} (1 - ν) .

Python]

def f(x,kappa):

return zeros(size(x))

def g(t,kappa):

return (sign(t)+1)/2

u=1;

Python]

FTCS=1; Upwind=2; LaxFriedrichs=3; LeapFrog=4; LaxWendroff=5; BeamWarming=6;

Python]

kappa=4

Python]

3.1.1Upwind

Python]

method=Upwind

from math import erf

m=99; h=1./(m+1); x=arange(m+2)*h;

dcfl=0.2; tfinal=0.5;

clf();

qex=g(tfinal-x/u,kappa);

plot(x[1:m],qex[1:m],'k.'); s = zeros(size(x));

for cfl in arange(dcfl,1+dcfl,dcfl):

k=h*cfl/u; nSteps=ceil(tfinal/k); t=arange(nSteps+1)*k; kappa=4;

Q0=f(x,kappa); Q1=zeros(size(x));

Qlft=g(t,kappa); Qrgt=zeros(size(t));

Q1=scheme(method,cfl,Q0,Q1,Qlft,Qrgt);

alpha = max(cfl*h**2*(1-cfl)/(2*k),10**(-6));

y = (u*tfinal-x)/sqrt(4*alpha*tfinal);

for i in range(m+1):

s[i] = 0.5*(1+erf(y[i]));

plot(x[1:m],Q1[1:m],'b',x[1:m],s[1:m],'r.');

xlabel('x'); ylabel('q'); title('Shock propagation by upwind scheme');

savefig("Lab05UpwindShockInitCond.pdf");

Python]

Figure 2. The exact solution to the modified equation corresponds closely to the results from the upwind scheme