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}.
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?
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.
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
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 \):
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 \):
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.
b) Write out the first four terms of the Taylor expansion of \( u_j^{n \pm 1} \), and \( u_{j \pm 1}^n \).
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.
d) Find an expression for the stability limit of the scheme in eq. \eqref{eq:leapfrog} using von Neumann method.
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 \).
e) Consider now the case that \( \lvert c\rvert \leq 1 \). What is the expression for the diffusive error \( \epsilon_D \)?
Use the results from the previous sub-exercise.
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*} $$
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} \)
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 \).
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?
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*} $$