Theory and python exercise 6

Fredrik Eikeland Fossan


Apr 11, 2019


ODEs

Problem 1: Local truncation error and absolute stability of a simple model equation (theory)

Consider the simple model equation $$ \begin{align} y^{\prime} = \lambda \,y\,, \label{eq:simple_eq} \end{align} $$ where \( \lambda \) is a constant smaller than zero (\( \lambda < 0 \)), and where the initial condition is given by $$ \begin{align} y\left(0\right) = 1\,. \label{eq:init} \end{align} $$ The analytical solution is given by: $$ \begin{align} y\left(x\right) = e^{\lambda x}\,. \label{eq:analytical} \end{align} $$

a) Derive the local truncation error for euler's method and for Heun's method applied on Eqs. \eqref{eq:simple_eq} -\eqref{eq:init}.

Solution.

The approximate differential operator for Euler's scheme is $$ \begin{align*} L^{euler}_a = y_{n+1} - [y_n + h \lambda y_n], \end{align*} $$

whereas the LTE can be computed by inserting \( y(x) \) $$ \begin{align*} \tau_n = & \frac{1}{h}\left\{ L_a(y(x_n))\right\} = \frac{1}{h}\left\{ y(x_{n+1}) - [y(x_n) + h \lambda y(x_n)]\right\}\\ = & \frac{1}{h}\left\{ y(x_n) + h y'(x_n) + \frac{h^2}{2} y''(x_n) + \ldots + \frac{1}{p!} h^p y^{(p)}(x_n) - y(x_n) - h \lambda y(x_n) \right\} \\ = & \frac{1}{2} h y''(x_n) + \ldots + \frac{1}{p!} h^{p-1} y^{(p)}(x_n) \\ \approx & \frac{1}{2} h y''(x_n). \end{align*} $$

For Eq. \eqref{eq:simple_eq}, heun's method is given by $$ \begin{align*} y_{n + 1}^p = y_n + h \lambda y_n \,,\\ y_{n + 1} = y_n + \frac{h}{2}\left( \lambda y_n + \lambda y_{n + 1}^p\right) \\ =y_n + \frac{h}{2}\left( \lambda y_n + \lambda y_n + h \lambda y_n\right) \\ = y_n + \frac{h}{2}\left( \lambda y_n + \lambda y_n + h \lambda^2 y_n\right) \\ = y_n + h \lambda y_n + \frac{h^2}{2}\,\lambda^2\,y_n \end{align*} $$

The approximate differential operator for Heun's scheme may then be formulated as $$ \begin{align*} L^{heun}_a = y_{n+1} - [y_n + h \lambda y_n + \frac{h^2}{2}\,\lambda^2\,y_n], \end{align*} $$

whereas the LTE can be computed by inserting \( y(x) \) $$ \begin{align*} \tau_n = & \frac{1}{h}\left\{ L_a(y(x_n))\right\} = \frac{1}{h}\left\{ y(x_{n+1}) - [y(x_n) + h \lambda y(x_n) + \frac{h^2}{2}\,\lambda^2\,y(x_n)]\right\}\\ = & \frac{1}{h}\left\{ y(x_n) + h y'(x_n) + \frac{h^2}{2} y''(x_n) + \frac{h^3}{6} y'''(x_n) + \text{H.O.T} - y(x_n) - h \lambda y(x_n) - \frac{h^2}{2}\,\lambda^2\,y(x_n) \right\} \\ = & \frac{1}{6} h^2 y'''(x_n) + \ldots + \frac{1}{p!} h^{p-1} y^{(p)}(x_n) \\ \approx & \frac{1}{6} h^2 y'''(x_n). \end{align*} $$

b) What are the numerical amplification factors for euler's method and Heun's method applied on Eqs. \eqref{eq:simple_eq} -\eqref{eq:init}? What are the accompanying requirements for absolute stability?

Solution.

The numerical amplification factor \( G \) is given by $$ \begin{align*} G = \frac{y_{n +1}}{y_n} \end{align*} $$

which gives $$ \begin{align*} G^{\text{euler}} = \frac{y_n + h \lambda y_n}{y_n} = 1 + h \lambda\,. \end{align*} $$ The stability criteria \( \lvert G^{\text{euler}} \rvert \leq 1 \) then gives the following stability region \( -2 < \lambda h < 0 \).

For Heun's method we have $$ \begin{align*} G^{\text{heun}} = \frac{y_n + h \lambda y_n + \frac{h^2}{2} \lambda^2 y_n}{y_n} = 1 + h \lambda + \frac{h^2 \lambda^2}{2}\,, \end{align*} $$ The stability criteria \( \lvert G^{\text{heun}} \rvert \leq 1 \) then gives the following stability region \( -2 < \lambda h < 0 \).

Problem 2: Local truncation error and absolute stability of a simple model equation (python)

a) Set \( \lambda = -1 \), and perform a grid refinement test to compute the observed order of accuracy for Euler's method and for Heun's method. Choose a step-size \( h=\Delta x = 0.25 \) for the coarsest grid and refine 5 times with a refinement factor of 2. Compare the observed order with what is expected based on the local truncation error.

Hint.

The observed order of accuracy of a numerical solution can be estimated as: $$ \begin{align*} p = \frac{log\left(\frac{\epsilon_{n - 1}}{\epsilon_n}\right)}{log\left(r\right)} \,, \end{align*} $$ where \( \epsilon_{n - 1} \) is an error calculated on a grid with step-size \( h \) and \( \epsilon_{n} \) is an error calculated on a grid with step-size \( \frac{h}{r} \). For convenience we may choose \( r=2 \). The observed error may thus be estimated by comparing the error on grids which are successively refined by a factor of 2. The maximum absolute error, \( max(\lvert y\left(x_i\right)-y_i\rvert) \), could be used as the error metric for each grid. Here, \( y\left(x_i\right) \) refers to the analytical solution at location \( x_i \) and \( y_i \) is the numerical approximation at the same point. The procedure may be summarized as follows

  • discretize the domain with current step-size
  • compute the numerical solution \( y_n \) with given numerical scheme
  • estimate the current error, \( \epsilon_n \)
  • if the current solution is not obtained on the first/coursest grid, estimate the observed order, p
  • reduce the step-size with a factor of 2
  • repeat the above steps \( N \) times (\( N=5 \)).

Solution. 

from ODEschemes import euler, heun, rk4
import numpy as np
import matplotlib.pylab as plt

def func(y, x):
    return [l*y] 

l = -1
h = 0.25
x = np.arange(0, 1 + h, h)
y0 = [1]

y = euler(func, y0, x)[:, 0]
y_analytical = np.exp(l*x)


plt.figure()
plt.plot(x, y)
plt.plot(x, y_analytical, 'k--')
plt.show()
n_refine = 5

observedErrorList = []
error_list = []
for n in range(n_refine):

    y = euler(func, y0, x)[:, 0]
    y_analytical = np.exp(l*x)
    
    eps = np.max(np.abs(y - y_analytical))
    error_list.append(eps)
    if n > 0:
        observedErrorList.append(np.log2(error_list[-2]/error_list[-1]))
    h *= 0.5

    plt.figure()
    plt.plot(x, y)
    plt.plot(x, y_analytical, 'k--')
    plt.show()
    
    x = np.arange(0, 1 + h, h)
print observedErrorList

b) Compute the observed order of accuracy, as in the previous sub-exercise, for the rk4 method (optional).

c) Choose the euler method and set \( \lambda=-10 \) and plot the numerical solution (together with the analytical) for each of the following step-sizes, \( h=\Delta x \):

For each solution, argue if the solution is stable and/or unstable and whether it is oscillatory. Compare the observations with what you would expect based on the numerical amplification factor.

Solution. 

from ODEschemes import euler, heun, rk4
import numpy as np
import matplotlib.pylab as plt

def func(y, x):
    return [l*y] 

l = -10

y0 = [1]



h_list = [0.05, 0.1, 0.15, 0.2, 0.25, 0.33]

for h in h_list:
    x = np.arange(0, 1 + h, h)
    y = euler(func, y0, x)[:, 0]
    y_analytical = np.exp(l*x)
    plt.figure()
    plt.plot(x, y)
    plt.plot(x, y_analytical, 'k--')
    plt.title("h = {}".format(h))
    plt.show()

d) Choose Heun's method and set \( \lambda=-10 \) and plot the numerical solution (together with the analytical) for each of the following step-sizes, \( h=\Delta x \):

For each solution, argue if the solution is stable and/or unstable and whether it is oscillatory. Compare the observations with what you would expect based on the numerical amplification factor.

PDEs

Problem 3: Numerical analysis on schemes for the linear advection equation (theory)

The linear advection equation takes the form: $$ \begin{align} \frac{\partial u}{\partial t} + a \, \frac{\partial u}{\partial x} = 0, \label{eq:advection} \end{align} $$

where \( a \) is the wave speed. The analytical solution is given by $$ \begin{equation}\label{eq:hyp_analytical} u = u_0(x-a\,t ), \end{equation} $$ In which \( u_0(x) \) is an initial solution/wave (\( u(x, 0)=u_0(x) \)).

a) Discretize eq. \eqref{eq:advection} using central differences for both the temporal and spatial derivative, and show that the corresponding scheme (called the Leapfrog scheme) may be written: $$ \begin{align} u_j^{n+1} = u_j^{n-1} - c \left(u_{j+1}^n - u_{j-1}^n\right), \label{eq:leapfrog} \end{align} $$

where \( c=a\frac{\Delta t}{\Delta x} \) is the Courant number.

Solution.

Using centered finite differences for time and space derivatives in \eqref{eq:advection} yields $$ \begin{align*} \frac{u_j^{n+1} - u_j^{n-1}}{2 \Delta t} = - a \frac{u_{j+1}^n - u_{j-1}^n}{2 \Delta x}, \end{align*} $$ that can be rearranged to $$ \begin{align*} u_j^{n+1} = u_j^{n-1} - a\frac{\Delta t}{\Delta x} \left(u_{j+1}^n - u_{j-1}^n\right), \end{align*} $$ and finally yield $$ \begin{align*} u_j^{n+1} = u_j^{n-1} - c \left(u_{j+1}^n - u_{j-1}^n\right). \end{align*} $$

b) Write out the first four terms of the Taylor expansion of \( u_j^{n \pm 1} \), and \( u_{j \pm 1}^n \).

Solution.

Taylor series expansion in time read $$ \begin{align} u_j^{n \pm 1} = u_j^{n} \pm \Delta t \left. \frac{\partial u}{\partial t}\right|_j^n + \frac{{\Delta t}^2}{2} \left. \frac{\partial^2 u}{\partial t^2}\right|_j^n \pm \frac{{\Delta t}^3}{6} \left. \frac{\partial^3 u}{\partial t^3}\right|_j^n + O(h^4), \label{eq:taylor3b1} \end{align} $$ while expansions in space are given by $$ \begin{align} u_{j \pm 1}^{n} = u_j^{n} \pm \Delta x \left. \frac{\partial u}{\partial x}\right|_j^n + \frac{{\Delta x}^2}{2} \left. \frac{\partial^2 u}{\partial x^2}\right|_j^n \pm \frac{{\Delta x}^3}{6} \left. \frac{\partial^3 u}{\partial x^3}\right|_j^n + O(h^4). \label{eq:taylor3b2} \end{align} $$

c) Find an expression for the local truncation error of the scheme in eq. \eqref{eq:leapfrog}, and show that the scheme is second order in time and space.

Hint.

$$ \begin{align*} \tau_u = \frac{{\Delta t}^2}{3} \left. \frac{\partial^3 u}{\partial t^3}\right|_j^n + a \frac{{\Delta x}^2}{3} \left. \frac{\partial^3 u}{\partial x^3}\right|_j^n\,. \end{align*} $$

Solution.

The numerical operator for scheme \eqref{eq:leapfrog} is $$ \begin{align} \label{eq:LeapOperator} L_u(u_{i}) = \frac{u^{n+1}_{j} -u^{n-1}_{j}}{\Delta t} + a \frac{u^{n}_{j+1} -u^{n}_{j-1}}{ \Delta x}, \end{align} $$ and its local truncation error can be then defined as $$ \begin{align} \tau_u = L_u(u(x_j,t^n)) \label{eq:LeapLTE}, \end{align} $$ which, using Taylor series expansions, yields $$ \begin{align*} \begin{split} \tau_u = & \frac{\left( u_j^{n} + \Delta t \left. \frac{\partial u}{\partial t}\right|_j^n + \frac{{\Delta t}^2}{2} \left. \frac{\partial^2 u}{\partial t^2}\right|_j^n + \frac{{\Delta t}^3}{6} \left. \frac{\partial^3 u}{\partial t^3}\right|_j^n + O(h^4)\right) - \left( u_j^{n} - \Delta t \left. \frac{\partial u}{\partial t}\right|_j^n + \frac{{\Delta t}^2}{2} \left. \frac{\partial^2 u}{\partial t^2}\right|_j^n - \frac{{\Delta t}^3}{6} \left. \frac{\partial^3 u}{\partial t^3}\right|_j^n + O(h^4)\right)}{\Delta t} + \\ & + a \frac{\left( u_j^{n} + \Delta x \left. \frac{\partial u}{\partial x}\right|_j^n + \frac{{\Delta x}^2}{2} \left. \frac{\partial^2 u}{\partial x^2}\right|_j^n + \frac{{\Delta x}^3}{6} \left. \frac{\partial^3 u}{\partial x^3}\right|_j^n + O(h^4)\right) - \left( u_j^{n} - \Delta x \left. \frac{\partial u}{\partial x}\right|_j^n + \frac{{\Delta x}^2}{2} \left. \frac{\partial^2 u}{\partial x^2}\right|_j^n - \frac{{\Delta x}^3}{6} \left. \frac{\partial^3 u}{\partial x^3}\right|_j^n + O(h^4)\right)}{\Delta x},\\ \end{split} \end{align*} $$

in which terms with \( u_j^n, \frac{\partial u^2}{\partial t^2}, \frac{\partial u^2}{\partial x^2} \) cancel. Continuing without including higher order terms yields: $$ \begin{align*} \begin{split} \tau_u & = \frac{\left(2 \Delta t \left. \frac{\partial u}{\partial t}\right|_j^n + \frac{{\Delta t}^3}{3} \left. \frac{\partial^3 u}{\partial t^3}\right|_j^n \right) }{\Delta t} + a \frac{\left( 2 \Delta x \left. \frac{\partial u}{\partial x}\right|_j^n + \frac{{\Delta x}^3}{3} \left. \frac{\partial^3 u}{\partial x^3}\right|_j^n \right) }{\Delta x},\\ & = 2 \left( \left.\frac{\partial u}{\partial t}\right|_j^n + a \left.\frac{\partial u}{\partial x}\right|_j^n \right) + \frac{{\Delta t}^2}{3} \left. \frac{\partial^3 u}{\partial t^3}\right|_j^n + a \frac{{\Delta x}^2}{3} \left. \frac{\partial^3 u}{\partial x^3}\right|_j^n,\\ & = \frac{{\Delta t}^2}{3} \left. \frac{\partial^3 u}{\partial t^3}\right|_j^n + a \frac{{\Delta x}^2}{3} \left. \frac{\partial^3 u}{\partial x^3}\right|_j^n, \end{split} \end{align*} $$ which shows that the local truncation error has an order \( O({\Delta t}^2, {\Delta x}^2) \).

d) Find an expression for the stability limit of the scheme in eq. \eqref{eq:leapfrog} using von Neumann method.

Hint.

After some algebra one obtains $$ \begin{align*} G^{2} + 2 i \,G \, c \,sin(\delta) - 1 = 0\,, \end{align*} $$ on which you may apply the quadratic formula. Further separate between cases where \( \vert c \rvert \leq 1 \) and \( \rvert c \rvert > 1 \). For the case \( \rvert c \rvert > 1 \) you may assume \( sin \, \delta = 1 \).

Solution.

Substituting \( u_j^n \rightarrow E_j^n = G^n e^{i\beta x_j} \) in eq. \eqref{eq:leapfrog} $$ \begin{align*} G^{n+1} e^{i\beta x_j} = G^{n-1} e^{i\beta x_j} - c \left(G^n e^{i\beta x_{j+1}} - G^n e^{i\beta x_{j-1}}\right), \end{align*} $$ and further division by \( G^{n-1} e^{i\beta x_j} \), and introducing \( \delta = \beta (x_j-x_{j-1}) \) gives $$ \begin{align*} G^{2} = 1 - c \left(G e^{i \delta} - G e^{-i \delta}\right). \end{align*} $$

Noting that the following trigonometric relation holds: $$ \begin{align*} e^{ix} - e^{-ix} = i 2 sin(x), \end{align*} $$ we obtain $$ \begin{align*} G^{2} + 2 i \,G \, c \,sin(\delta) - 1 = 0. \end{align*} $$

Further by employing the quadratic formula we get the following roots for \( G \) $$ \begin{align*} G_\pm = -i \, c \,sin(\delta) \pm \sqrt{1 - c^2 \, sin^2(\delta)}. \end{align*} $$ We have two cases. For \( |c|\leq 1 \), one can see that $$ \begin{align*} |G|^2 = c^2\,sin^2(\delta) + 1 - c^2\,sin^2(\delta) = 1, \end{align*} $$ so that the scheme is conditionally stable for \( |c|\leq 1 \). On the other hand, if \( |c|>1 \), there will be some wavelengths for which \( c \, sin(\delta)>1 \), in which case the two roots are purely imaginary. Consider the case \( c>1 \) and \( sin(\delta)=1 \), then $$ \begin{align*} G_\pm = -i c \pm i \sqrt{c^2-1} = -i \left( c \mp \sqrt{c^2-1} \right) , \end{align*} $$ for which \( G_->1 \), resulting in an unstable scheme.

e) Consider now the case that \( \lvert c\rvert \leq 1 \). What is the expression for the diffusive error \( \epsilon_D \)?

Hint.

Use the results from the previous sub-exercise.

Solution.

$$ \begin{align*} \epsilon_D = \lvert G\rvert ) = 1, \end{align*} $$ so there should be no amplitude errors, however there might still be round-of errors etc.

f) Briefly explain the concept of the modified equation, and show that the modified equation for the Leapfrog scheme may be expressed as: $$ \begin{align*} \frac{\partial u}{\partial t} + a \frac{\partial u}{\partial x} = \frac{{\Delta x}^2 \cdot a}{6} \left(c^2 - 1\right) \, \frac{\partial^3 u}{\partial x^3} + \text{H.O.T.} \end{align*} $$

Hint.

By applying the Cauchy-Kowalewsky procedure we can express to express the temporal derivative, \( \frac{\partial^3 u}{\partial t^3} \) in terms of spatial derivatives \( \frac{\partial^3 u}{\partial t^3} = -a^3\frac{\partial^3 u}{\partial x^3} \)

Solution.

The Leapfrog scheme is given by $$ \begin{align*} \frac{u^{n+1}_{j} -u^{n-1}_{j}}{\Delta t} + a \frac{u^{n}_{j+1} -u^{n}_{j-1}}{ \Delta x} = 0\,. \end{align*} $$

Inserting Taylor series expansions, yields $$ \begin{align*} \begin{split} & \frac{\left( u_j^{n} + \Delta t \left. \frac{\partial u}{\partial t}\right|_j^n + \frac{{\Delta t}^2}{2} \left. \frac{\partial^2 u}{\partial t^2}\right|_j^n + \frac{{\Delta t}^3}{6} \left. \frac{\partial^3 u}{\partial t^3}\right|_j^n + O(h^4)\right) - \left( u_j^{n} - \Delta t \left. \frac{\partial u}{\partial t}\right|_j^n + \frac{{\Delta t}^2}{2} \left. \frac{\partial^2 u}{\partial t^2}\right|_j^n - \frac{{\Delta t}^3}{6} \left. \frac{\partial^3 u}{\partial t^3}\right|_j^n + O(h^4)\right)}{\Delta t} \\ & + a \frac{\left( u_j^{n} + \Delta x \left. \frac{\partial u}{\partial x}\right|_j^n + \frac{{\Delta x}^2}{2} \left. \frac{\partial^2 u}{\partial x^2}\right|_j^n + \frac{{\Delta x}^3}{6} \left. \frac{\partial^3 u}{\partial x^3}\right|_j^n + O(h^4)\right) - \left( u_j^{n} - \Delta x \left. \frac{\partial u}{\partial x}\right|_j^n + \frac{{\Delta x}^2}{2} \left. \frac{\partial^2 u}{\partial x^2}\right|_j^n - \frac{{\Delta x}^3}{6} \left. \frac{\partial^3 u}{\partial x^3}\right|_j^n + O(h^4)\right)}{\Delta x} = 0,\\ \end{split} \end{align*} $$

in which terms with \( u_j^n, \frac{\partial u^2}{\partial t^2}, \frac{\partial u^2}{\partial x^2} \) cancel. Continuing without including higher order terms yields: $$ \begin{align*} \begin{split} & \frac{\left(2 \Delta t \left. \frac{\partial u}{\partial t}\right|_j^n + \frac{{\Delta t}^3}{3} \left. \frac{\partial^3 u}{\partial t^3}\right|_j^n \right) }{\Delta t} + a \frac{\left( 2 \Delta x \left. \frac{\partial u}{\partial x}\right|_j^n + \frac{{\Delta x}^3}{3} \left. \frac{\partial^3 u}{\partial x^3}\right|_j^n \right) }{\Delta x}\\ & = 2 \left( \left.\frac{\partial u}{\partial t}\right|_j^n + a \left.\frac{\partial u}{\partial x}\right|_j^n \right) + \frac{{\Delta t}^2}{3} \left. \frac{\partial^3 u}{\partial t^3}\right|_j^n + a \frac{{\Delta x}^2}{3} \left. \frac{\partial^3 u}{\partial x^3}\right|_j^n = 0\\ & \rightarrow \left.\frac{\partial u}{\partial t}\right|_j^n + a \left.\frac{\partial u}{\partial x}\right|_j^n = - \frac{{\Delta t}^2}{6} \left. \frac{\partial^3 u}{\partial t^3}\right|_j^n - a \frac{{\Delta x}^2}{6} \left. \frac{\partial^3 u}{\partial x^3}\right|_j^n, \end{split} \end{align*} $$ By applying the Cauchy-Kowalewsky procedure we can express the temporal derivative, \( \frac{\partial^3 u}{\partial t^3} \) in terms of spatial derivatives, \( \frac{\partial^3 u}{\partial t^3} = -a^3\frac{\partial^3 u}{\partial x^3} \), such that modified equation becomes: $$ \begin{align*} \frac{\partial u}{\partial t} + a \frac{\partial u}{\partial x} = \frac{{\Delta x}^2 \cdot a}{6} \left(c^2 - 1\right) \, \frac{\partial^3 u}{\partial x^3} + \text{H.O.T.} \end{align*} $$

Problem 4: Numerical analysis on schemes for the linear advection equation (python)

Consider a smooth sine squared wave between \( x=0 \) and \( x=0.2 \) as the initial condition: $$ \begin{align}\label{eq:initialWave} u_0(x) = sin^2(\pi(x/0.2)), \quad 0 < x < 0.2, \\ u_0(x) = 0, \quad 0.2 \leq x \leq l, \label{_auto1} \end{align} $$

where \( l = 1 \). Discretize the domain with equal spacing, \( \Delta x = \frac{l}{N} \) and calculate \( \Delta t \) from the CFL condition. Choose \( N=200 \) and Simulate from \( t_{start}=0 \) to \( t_{end}=0.6 \). This way we avoid the need to introduce special boundary conditions. (The wave will not leave the domain, and \( u(0)=0 \) and \( u(l)=0 \)). Otherwise the following boundarycondtions should be applied at the right boundary: $$ \begin{equation}\label{eq:bc} u(l)^{n + 1} = u(l - a \cdot \Delta t)^n, \end{equation} $$

a) Compare the solution obtained with the leapfrog scheme with the analytical solution for the linear advection equation. choose \( a=1 \) and a cfl-constraint condition of \( c=0.5 \).

Note that the leapfrog scheme is a three time-level scheme, and thus require an initial condition at \( n=0 \), and also at \( n=1 \). Use the ftbs scheme (see python exercise 5) to obtain the solution at the first time level, \( n=1 \).

Solution. 

'''
Created on May 16, 2017

@author: fredrik
'''
#### Start import 
import numpy as np
import matplotlib.pyplot as plt
from animationHyp import myHypAnimation
from scipy import interpolate

def f(x):
    """A smooth sine^2 function between x_left and x_right"""
    f = np.zeros_like(x)
    x_left = 0
    x_right = 0.2
    xm = (x_right-x_left)/2.0
    f = np.where((x>x_left) & (x<x_right), np.sin(np.pi*(x-x_left)/(x_right-x_left))**2.0,f) 
    return f

def f2(x):
    """A squared box-function between x_left and x_right"""
    f = np.zeros_like(x)
    x_left = 0.0
    x_right = 0.2
    f = np.where((x>x_left) & (x<x_right), 1,f) 
    return f

def ftbs(U_initial, c):
    """ Forward in time backward in space scheme"""
    
    Unew = U_initial[1:-1] - c*(U_initial[1:-1] - U_initial[:-2])
    
    return Unew

def leapfrog(Uold, Uolder, c):
    """ Leapfrog scheme"""
    
    Unew = Uolder[1:-1] - c*(Uold[2:] - Uold[:-2])
    
    return Unew

def lax_wendroff(u, c): 
    u[1:-1] = c/2.0*(1+c)*u[:-2] + (1-c**2)*u[1:-1] - c/2.0*(1-c)*u[2:]
    return u[1:-1]


f = f
a = 1.0 # linear wavespeed
tmin, tmax = 0.0, 0.6 # start and stop time of simulation
xmin, xmax = 0.0, 1.0 # start and end of spatial domain
Nx = 200 # number of spatial points
c = 0.5 # courant number, 

# discretize
x = np.linspace(xmin, xmax, Nx+1) # discretization of space
dx = x[1] - x[0] # spatial step size
dt = c/a*dx # stable time step calculated from stability requirement
Nt = int((tmax-tmin)/dt) # number of time steps
time = np.linspace(tmin, tmax, Nt) # discretization of time
print a*dt/dx
# initialize variables

u_initial = f(x)
u = np.zeros_like(x)
un = np.zeros((len(time), len(x))) # holds the numerical solution
u1 = np.zeros_like(u_initial)

u_analytical = np.zeros((len(time),len(x)))
u_analytical[0, :] = f(x)
u_analytical[1, :] = f(x - a*dt)
u1[1:-1] = ftbs(u_initial, c) # Since the leapfrog-scheme is a two time-step method we need an initial condition also for the first time-step
un[0, :] = u_initial
un[1, :] = u1#f(x - a*dt)

for n in range(1, Nt - 1):
        
    u_bc = interpolate.interp1d(x[-2:], u[-2:]) # interpolate at right boundary
    u[1:-1] = leapfrog(un[n, :], un[n - 1, :], c) # calculate numerical solution of interior
    #u[1:-1] = ftbs(un[n, :], c) # uncomment to use ftbs instead of leapfrog
    u[-1] = u_bc(x[-1] - a*dt) # interpolate along a characteristic to find the boundary value
    
    un[n + 1,:] = u # storing the solution for plotting
    u_analytical[n + 1, :] = f(x - a*time[n + 1])

plt.figure()
plt.plot(x, un[-1, :])
plt.plot(x, u_analytical[-1, :], 'k--')
myHypAnimation(x, un, u_analytical, [leapfrog], len(time))

b) Solve the same problem as the previous sub-exercise with the ftbs scheme.

c) The modified equation for the ftbs is given by: $$ \begin{align*} \frac{\partial u}{\partial t} + a \frac{\partial u}{\partial x} = \frac{a \Delta x}{h}\left(1 - c\right)\frac{\partial^2 u}{\partial x^2} + \text{H.O.T.}\,. \end{align*} $$ Furthermore, the diffusion error is given by: $$ \begin{align*} \epsilon_D = \sqrt{1 - 4\,c\left(1 - c\right)\,sin\,\left(\frac{\delta}{2}\right)}\,. \end{align*} $$

Use these expressions together with the corresponding expressions for the leapfrog scheme (see theory exercise above) to explain the numerical errors you observe for the two schemes. Is it possible to say anything about the stability limit of the ftbs scheme based on the modified equation?

Solution.

The leading term of the RHS of the modified equation for the leapfrog-scheme is \( \frac{{\Delta x}^2 \cdot a}{6} \left(c^2 - 1\right) \, \frac{\partial^3 u}{\partial x^3} \), which can be recognized as a dispersive terms, and the errors should thus be reflected in terms of error in timing and phase in the wave. As one can see from the simulations, the numerical wavespeed is slightly different than the analytical, and spurious oscillations develops in regions which sharp gradients. There is however minimal amplitude-error for the peak, which reflects the fact that \( \epsilon_D = \lvert G\rvert = 1 \) for the leapfrog-scheme. For the ftbs scheme on the other hand, the leading term of the RHS of the modified equation is \( \frac{a \Delta x}{h}\left(1 - c\right)\frac{\partial^2 u}{\partial x^2} \), which can be recognized as a diffusive terms, in which the coefficient, \( \frac{a \Delta x}{h}\left(1 - c\right) \) can be seen upon as a numerical viscosity. If \( c > 1 \) the numerical viscosity is negative, and amplitudes will grow (it will be unstable). For \( c < 1 \) on the other hand we see from the simulations that the ftbs scheme is very diffusive.

d) Run the above analysis on a squared initial wave: $$ \begin{align*} u_0(x) = 1, \quad 0 < x < 0.2, \\ u_0(x) = 0, \quad 0.2 \leq x \leq l, \end{align*} $$

Solution.

dependencies 

# src-ch1/ODEschemes.py

import numpy as np
from matplotlib.pyplot import plot, show, legend, hold,rcParams,rc, figure, axhline, close,\
    xticks, title, xlabel, ylabel, savefig, axis, grid, subplots, setp

# change some default values to make plots more readable 
LNWDT=3; FNT=10
rcParams['lines.linewidth'] = LNWDT; rcParams['font.size'] = FNT
font = {'size' : 10}; rc('font', **font)


# define Euler solver
def euler(func, z0, time):
    """The Euler scheme for solution of systems of ODEs.
    z0 is a vector for the initial conditions,
    the right hand side of the system is represented by func which returns
    a vector with the same size as z0 ."""

    z = np.zeros((np.size(time), np.size(z0)))
    z[0,:] = z0

    for i in range(len(time)-1):
        dt = time[i+1] - time[i]
        z[i+1,:]=z[i,:] + np.asarray(func(z[i,:], time[i]))*dt

    return z


# define Heun solver
def heun(func, z0, time):
    """The Heun scheme for solution of systems of ODEs.
    z0 is a vector for the initial conditions,
    the right hand side of the system is represented by func which returns
    a vector with the same size as z0 ."""

    def f_np(z,t):
        """A local function to ensure that the return of func is an np array
        and to avoid lengthy code for implementation of the Heun algorithm"""
        return np.asarray(func(z,t))

    z = np.zeros((np.size(time), np.size(z0)))
    z[0,:] = z0
    zp = np.zeros_like(z0)

    for i, t in enumerate(time[0:-1]):
        dt = time[i+1] - time[i]
        zp = z[i,:] + f_np(z[i,:],t)*dt   # Predictor step
        z[i+1,:] = z[i,:] + (f_np(z[i,:],t) + f_np(zp,t+dt))*dt/2.0 # Corrector step

    return z


# define rk4 scheme
def rk4(func, z0, time):
    """The Runge-Kutta 4 scheme for solution of systems of ODEs.
    z0 is a vector for the initial conditions,
    the right hand side of the system is represented by func which returns
    a vector with the same size as z0 ."""

    z = np.zeros((np.size(time),np.size(z0)))
    z[0,:] = z0
    zp = np.zeros_like(z0)

    for i, t in enumerate(time[0:-1]):
        dt = time[i+1] - time[i]
        dt2 = dt/2.0
        k1 = np.asarray(func(z[i,:], t))                # predictor step 1
        k2 = np.asarray(func(z[i,:] + k1*dt2, t + dt2)) # predictor step 2
        k3 = np.asarray(func(z[i,:] + k2*dt2, t + dt2)) # predictor step 3
        k4 = np.asarray(func(z[i,:] + k3*dt, t + dt))   # predictor step 4
        z[i+1,:] = z[i,:] + dt/6.0*(k1 + 2.0*k2 + 2.0*k3 + k4) # Corrector step

    return z


if __name__ == '__main__':
    pass


'''
Created on Apr 25, 2018

@author: fredrik
'''
import matplotlib.pyplot as plt
from matplotlib import animation
import numpy as np
def myHypAnimation(x, un, uanalytical, solvers, Nt, xmin=0, xmax=1):
    fig = plt.figure()
    ax = plt.axes(xlim=(xmin,xmax), ylim=(np.min(un), np.max(un)*1.1))
    
    lines=[]     # list for plot lines for solvers and analytical solutions
    legends=[]   # list for legends for solvers and analytical solutions
    
    for solver in solvers:
        line, = ax.plot([], [], '-', lw=2)
        lines.append(line)
        legends.append(solver.func_name)
    
    line, = ax.plot([], [], '-', lw=1) #add extra plot line for analytical solution
    lines.append(line)
    legends.append('Analytical')
    
    plt.xlabel('x-coordinate [-]')
    plt.ylabel('Amplitude [-]')
    plt.legend(legends, loc=3, frameon=False)
     
    # initialization function: plot the background of each frame
    def init():
        for line in lines:
            line.set_data([], [])
        return lines,
    
    # animation function.  This is called sequentially
    def animate(i):
        for k, line in enumerate(lines):
            if (k==0):
                line.set_data(x, un[i,:])
            else:
                line.set_data(x, uanalytical[i,:])
        return lines,
     
    # call the animator.  blit=True means only re-draw the parts that have changed.
    anim = animation.FuncAnimation(fig, animate, init_func=init, frames=Nt, interval=10, blit=False)
    
    plt.show()