numerics for geofluid dynamics

A solution \(u^n_{j}\) is stable if the error \(u^n_j-u(j\Delta x,n\Delta t)\)remains bounded as n increases, for fixed values of \(\Delta x\) and \(\Delta t\).

Of course, we assume that also the true solution is bounded, which is usually the case. 

A finite difference scheme is stable if this is true for any initial condition. 

• only able to test the stability of linear equations or linearized versions
• most frequently used
• main concept: every solution to a linear equation can be expressed in a Fourier series, where each harmonic component is also a solution to the equation. Thus, it is enough to study the stability of one Fourier Component on it's own! That is: \(u_j^n=u_0e^{ik(j\Delta x-C_Dn\Delta t)}\)
• We analyze stability with the amplification factor \(\lambda\). If the true solution doesn't grow over time it is sufficient fo stability that: \(|\lambda|\leq1\).
• we define: \(U^{n+1}=\lambda U^{(n)}\)

In general, centred scheme has an accuracy of order 2, whereas uncentred scheme are only of order 1 in accuracy!

equation: \(\frac{dU}{dt}=i\omega U\)

analytical solution: \(U(t)=U(0)e^{i\omega t}\)

start with the stability analysis: \(U^{(n+1)}=\lambda U^{n}\)

We define \(\lambda=|\lambda|e^{i\Theta}\)

Then the numerical solution can formally be written as:

\(U^{(n)}=|\lambda|^nU^{(0)}e^{in\Theta}\)

\(|\lambda|(<,>,=) 1\) leads to: damping, unstable or neutral

The relative phase change leads to:

\(\frac{\Theta}{\omega\Delta t}(>=<)1\) leads to: accelerating, no effect, decelerating

For accuracy we want to have the ampflification factor and the relative phase speed close to unity. Unless it is a computational mode. Then we want strong damping!

Euler scheme:

• \(\lambda =1+ip\)
• unconditionally unstable

Backward Scheme:

• \(\lambda=\frac{1}{1+\omega\Delta t}(1+i\omega\Delta t)\)
• unconditionally stable
• it is always damping and damping increases with frequency.

Trapezoidal scheme:

• \(\lambda=\frac{1}{1+\frac{1}{4}p^2}(1-\frac{1}{4}p^2+ip)\)
• \(|\lambda|=1\)
• the trapezoidal scheme is always neutral (and thus stable)

Matsuno scheme:

• stable if: \(|\omega\Delta t|\leq 1\) (so time step must be small enough!)

Heun scheme:

• very weak unstable. With small time steps often the instability can be tolerated

Leap-frog scheme:

• stable and neutral if \(|\omega \Delta t|\leq1\)

Fazit: Ich would use leap frog because it has a high accuracy and under specific conditions it is stable and even neutral! Although it is implicit what is in general more complicated to deal with, so then trapezoidal would be an alternative: unconditionally neutral. Unconditionally unstable is the euler forward and thus not useable. The Heun Scheme is also unconditionally unstable but at some situation tolerable. 

• It's through the "round of error" that one can not eliminate the computational mode completely, even when it was avoided in the initail conditions. 
• Although " round of errors" are of little effect in atmospheric models. 
• It can be shown that the Heun scheme is a good choice to compute the computional initial condition, since it will give a smaller amplitude of the computational mode. 

We might like to use the leapfrog scheme for the oscillation term but we cannot use it for the friction term. So we use two different schemes with Euler forward for the friction!

\(\frac{du}{dt}=i\omega u-\gamma u\\\frac{u_j^{n+1}-u_j^{n-1}}{2\Delta t}=i\omega u^n-\gamma u^{n-1}\)

To define time schemes, we consider the equation:

\(\frac{du}{dt}=f(u,t)\;\;\;where\;\;u=u(t)\)

1. Two-level-schemes: use timesteps n and n+1
   1. \(u^{n+1}=u^n+\int^{(n+1)\Delta t}_{n\Delta t}f(u,t)dt\)
   2. Euler or forward scheme:

      • \(u^{n+1}=u^n+\Delta tf^n\)

      • first order accurate scheme: \(O(\Delta t)\)

      • uncentred

      • explicit

To define time schemes, we consider the equation:

\(\frac{du}{dt}=f(u,t)\;\;\;where\;\;u=u(t)\)

1. Two-level-schemes: use timesteps n and n+1
   1. \(u^{n+1}=u^n+\int^{(n+1)\Delta t}_{n\Delta t}f(u,t)dt\)
   2. Euler or forward scheme:

      • \(u^{n+1}=u^n+\Delta tf^n\)

      • first order accurate scheme: \(O(\Delta t)\)

      • uncentred

      • explicit

   3. Backward scheme:

      • \(u^{n+1}=u^n+\Delta tf^{n+1}\)

      • first order accurate scheme: \(O(\Delta t)\)

      • uncentred
      • implicit 
   4. Trapezoidal scheme
      • \(u^{n+1}=u^n+\frac{1}{2}\Delta t(f^n+f^{n+1})\)

      • second order accurate scheme: \(O((\Delta t)^2)\)

      • implicit
   5. Matsuno or Euler-backward scheme
      • to increase the accuracy: first an Euler forward time step, than a backward step
      • \(1.\;u^{n+1/2}=u^n+\Delta t f^n\\ 2.\;u^{n+1}=u^n+\Delta tf^{n+1/2}\)
      • explicit

      • first order accurate scheme: \(O(\Delta t)\)

   6. Heun scheme
      • first Euler forward, then Trapezoidal
      • \(1.\;U^{n+1*}=U^{n}+\Delta t*f^n\\2.\;U^{n+1}=U^n+\frac{1}{2}\Delta t(f^n+f^{n+1})\)
      • second order of accuracy
      • explicit 

1. Three-level-schemes: use timesteps n-1, n, n+1
   1. \(u^{n+1}=u^{n-1}+\int^{(n+1)\Delta t}_{(n-1)\Delta t}f(u,t)dt\)
   2. leapfrog scheme
      • \(u^{n+1}=u^{n-1}+2\Delta tf^n\)

      • second order accurate scheme: \(O((\Delta t)^2)\)

      •  

        most widely used scheme in atmospheric and ocean models. 

To define time schemes, we consider the equation:

\(\frac{du}{dt}=f(u,t)\;\;\;where\;\;u=u(t)\)

1. Two-level-schemes: use timesteps n and n+1
   1. \(u^{n+1}=u^n+\int^{(n+1)\Delta t}_{n\Delta t}f(u,t)dt\)
   2. Euler or forward scheme:

      • \(u^{n+1}=u^n+\Delta tf^n\)

      • first order accurate scheme: \(O(\Delta t)\)

      • uncentred

      • explicit

   3. Backward scheme:

      • \(u^{n+1}=u^n+\Delta tf^{n+1}\)

      • first order accurate scheme: \(O(\Delta t)\)

      • uncentred
      • implicit 
   4. Trapezoidal scheme
      • \(u^{n+1}=u^n+\frac{1}{2}\Delta t(f^n+f^{n+1})\)

      • second order accurate scheme: \(O((\Delta t)^2)\)

      • implicit
   5. Matsuno or Euler-backward scheme
      • to increase the accuracy: first an Euler forward time step, than a backward step
      • \(1.\;u^{n+1/2}=u^n+\Delta t f^n\\ 2.\;u^{n+1}=u^n+\Delta tf^{n+1/2}\)
      • explicit

      • first order accurate scheme: \(O(\Delta t)\)

   6. Heun scheme
      • first Euler forward, then Trapezoidal
      • \(1.\;U^{n+1*}=U^{n}+\Delta t*f^n\\2.\;U^{n+1}=U^n+\frac{1}{2}\Delta t(f^n+f^{n+1})\)
      • second order of accuracy
      • explicit 

1. Three-level-schemes: use timesteps n-1, n, n+1
   1. \(u^{n+1}=u^{n-1}+\int^{(n+1)\Delta t}_{(n-1)\Delta t}f(u,t)dt\)
   2. leapfrog scheme
      • \(u^{n+1}=u^{n-1}+2\Delta tf^n\)

      • second order accurate scheme: \(O((\Delta t)^2)\)

      •  

        most widely used scheme in atmospheric and ocean models. 

To define time schemes, we consider the equation:

\(\frac{du}{dt}=f(u,t)\;\;\;where\;\;u=u(t)\)

1. Two-level-schemes: use timesteps n and n+1
   1. \(u^{n+1}=u^n+\int^{(n+1)\Delta t}_{n\Delta t}f(u,t)dt\)
   2. Euler or forward scheme:

      • \(u^{n+1}=u^n+\Delta tf^n\)

      • first order accurate scheme: \(O(\Delta t)\)

      • uncentred

      • explicit

   3. Backward scheme:

      • \(u^{n+1}=u^n+\Delta tf^{n+1}\)

      • first order accurate scheme: \(O(\Delta t)\)

      • uncentred
      • implicit 
   4. Trapezoidal scheme
      • \(u^{n+1}=u^n+\frac{1}{2}\Delta t(f^n+f^{n+1})\)

      • second order accurate scheme: \(O((\Delta t)^2)\)

      • implicit
   5. Matsuno or Euler-backward scheme
      • to increase the accuracy: first an Euler forward time step, than a backward step
      • \(1.\;u^{n+1/2}=u^n+\Delta t f^n\\ 2.\;u^{n+1}=u^n+\Delta tf^{n+1/2}\)
      • explicit

      • first order accurate scheme: \(O(\Delta t)\)

   6. Heun scheme
      • first Euler forward, then Trapezoidal
      • \(1.\;U^{n+1*}=U^{n}+\Delta t*f^n\\2.\;U^{n+1}=U^n+\frac{1}{2}\Delta t(f^n+f^{n+1})\)
      • second order of accuracy
      • explicit 

1. Three-level-schemes: use timesteps n-1, n, n+1
   1. \(u^{n+1}=u^{n-1}+\int^{(n+1)\Delta t}_{(n-1)\Delta t}f(u,t)dt\)
   2. leapfrog scheme
      • \(u^{n+1}=u^{n-1}+2\Delta tf^n\)

      • second order accurate scheme: \(O((\Delta t)^2)\)

      •  

        most widely used scheme in atmospheric and ocean models. 

To define time schemes, we consider the equation:

\(\frac{du}{dt}=f(u,t)\;\;\;where\;\;u=u(t)\)

1. Two-level-schemes: use timesteps n and n+1
   1. \(u^{n+1}=u^n+\int^{(n+1)\Delta t}_{n\Delta t}f(u,t)dt\)
   2. Euler or forward scheme:

      • \(u^{n+1}=u^n+\Delta tf^n\)

      • first order accurate scheme: \(O(\Delta t)\)

      • uncentred

      • explicit

   3. Backward scheme:

      • \(u^{n+1}=u^n+\Delta tf^{n+1}\)

      • first order accurate scheme: \(O(\Delta t)\)

      • uncentred
      • implicit 
   4. Trapezoidal scheme
      • \(u^{n+1}=u^n+\frac{1}{2}\Delta t(f^n+f^{n+1})\)

      • second order accurate scheme: \(O((\Delta t)^2)\)

      • implicit
   5. Matsuno or Euler-backward scheme
      • to increase the accuracy: first an Euler forward time step, than a backward step
      • \(1.\;u^{n+1/2}=u^n+\Delta t f^n\\ 2.\;u^{n+1}=u^n+\Delta tf^{n+1/2}\)
      • explicit

      • first order accurate scheme: \(O(\Delta t)\)

   6. Heun scheme
      • first Euler forward, then Trapezoidal
      • \(1.\;U^{n+1*}=U^{n}+\Delta t*f^n\\2.\;U^{n+1}=U^n+\frac{1}{2}\Delta t(f^n+f^{n+1})\)
      • second order of accuracy
      • explicit 

1. Three-level-schemes: use timesteps n-1, n, n+1
   1. \(u^{n+1}=u^{n-1}+\int^{(n+1)\Delta t}_{(n-1)\Delta t}f(u,t)dt\)
   2. leapfrog scheme
      • \(u^{n+1}=u^{n-1}+2\Delta tf^n\)

      • second order accurate scheme: \(O((\Delta t)^2)\)

      •  

        most widely used scheme in atmospheric and ocean models. 

To define time schemes, we consider the equation:

\(\frac{du}{dt}=f(u,t)\;\;\;where\;\;u=u(t)\)

1. Two-level-schemes: use timesteps n and n+1
   1. \(u^{n+1}=u^n+\int^{(n+1)\Delta t}_{n\Delta t}f(u,t)dt\)
   2. Euler or forward scheme:

      • \(u^{n+1}=u^n+\Delta tf^n\)

      • first order accurate scheme: \(O(\Delta t)\)

      • uncentred

      • explicit

   3. Backward scheme:

      • \(u^{n+1}=u^n+\Delta tf^{n+1}\)

      • first order accurate scheme: \(O(\Delta t)\)

      • uncentred
      • implicit 
   4. Trapezoidal scheme
      • \(u^{n+1}=u^n+\frac{1}{2}\Delta t(f^n+f^{n+1})\)

      • second order accurate scheme: \(O((\Delta t)^2)\)

      • implicit
   5. Matsuno or Euler-backward scheme
      • to increase the accuracy: first an Euler forward time step, than a backward step
      • \(1.\;u^{n+1/2}=u^n+\Delta t f^n\\ 2.\;u^{n+1}=u^n+\Delta tf^{n+1/2}\)
      • explicit

      • first order accurate scheme: \(O(\Delta t)\)

   6. Heun scheme
      • first Euler forward, then Trapezoidal
      • \(1.\;U^{n+1*}=U^{n}+\Delta t*f^n\\2.\;U^{n+1}=U^n+\frac{1}{2}\Delta t(f^n+f^{n+1})\)
      • second order of accuracy
      • explicit 

1. Three-level-schemes: use timesteps n-1, n, n+1
   1. \(u^{n+1}=u^{n-1}+\int^{(n+1)\Delta t}_{(n-1)\Delta t}f(u,t)dt\)
   2. leapfrog scheme
      • \(u^{n+1}=u^{n-1}+2\Delta tf^n\)

      • second order accurate scheme: \(O((\Delta t)^2)\)

      •  

        most widely used scheme in atmospheric and ocean models. 

The one-dimensional gravity waves: 

\(\frac{\partial u}{\partial t}=-g\frac{\partial h}{\partial x},\;\;\frac{\partial h}{\partial t}=-H\frac{\partial u}{\partial x}\)

The differential/difference equations for the one-dimensional case are: 

\(\frac{u_j}{\partial t}=-g\frac{h_{j+1}-h_{j-1}}{2\Delta x},\;\;\frac{\partial h_j}{\partial t}=-H\frac{u_{j+1}-u_{j-1}}{2\Delta x}\)

If we use for that a non-staggered grid (like the A-grid), than we will have two elementary subgrids, with the solution on one of these subgrids beeing completely decoupled from the other. Thus it would be better to calculate only one of these solutions, that is, to use a staggered grid (like the C-grid). 

This reduces computation time by a fator of two and the truncation error stays the same. £££

• We have shown, that the C-grid is the best choose for gravity waves. The error for the numerical solution can be decreased through the win in computation time. 
• Although, wheres the space derivatives are straight forward for the C-grid, the Corilis term is complicated to calculate. 
• To decide wheter which grid is to choose, we investigate the effect of the space distribution of dependent variables on the dispersive properties of the gravity-inertia waves. This will be done using the simplest centered approximations for the space derivatives, leaving the time derivatives in their differential form. That is the differential-diffrence equation.

The frequency equation for the one-dimesnional system (when wave soltions were plugged in) is: 

\((\frac{\vartheta}{f})^2=1+\frac{gH}{f^2}k^2\)

with the Rossby Radius of deformation is: 

\(R=\frac{\sqrt{gH}}{f}\)

It follows that the frquency of the gravity-inertia waves is monotonically increasing function of k. Therefore the group velocity is never zero. This is very important for the geostrophic adjustment process, as it precludes a local accumulation of wave energy. 

We now look at the effect of the finite differencing in space in this case. Obtaining again the frequency equation, we find that the non-dimensional frequency \(\vartheta/f\) depend on the two paramters: \(k\Delta x\) and \(R\Delta x\)

For the one-dimensional case C and B have also a monotinically increasing frequency as a function of wavenumber. Although the two-grid wave has a group velocity of zero. 

In the two-deminsional case, only the C-grid as no negative group velocities for all kombination of wavenumbers l and k. Although if the Rossby Radius is smaller than a half grid (\(\sqrt{gH}/f<\Delta x/2\)) than the group velocity is everywhere negative. However, this does not happen for typical grid sizes used in atmospheric models and thus the C-grid is the best choice to simulate the geostrophic adjustment process.  

Analytical:

\(\frac{\partial u}{\partial t}-c\frac{\partial u}{\partial x}\)

solution: \(u(t,x)=f(x-ct)\) 

-----------------

Euler forward: Only stable when with upstream space differencing. Introduces a strong damping. The CFL criterion must be also fulfilled. 

Leap-frog scheme: conditionally stable with the CFL criterion. 

Fazit: I would use leapfrog since it has a higher order of accuracy than the others. 

1. forward difference: \((\frac{du}{dx})_j\rightarrow \frac{u_{j+1}-u_j}{\Delta x}\)
2. centered difference: \((\frac{du}{dx})_j \rightarrow \frac{u_{j+1}-u_{j-1}}{2\Delta x}\)
3. backward difference: \((\frac{du}{dx})_j \rightarrow \frac{u_{j}-u_{j-1}}{\Delta x}\)

1. the error of the numerical solution is brought up by:
   1. inadequacy of the scheme 
   2. insufficient information about the initial consitions (only known at discrete time points) 

-> Thus, an increase of accuracy of the scheme improves only one of these two components, and the result is not too impressive. 

2.  for stability requirements, it is necessary to choose a time step significantly smaller than that required for adequate accuracy.

-> With timesteps usually chosen, other errors, for example in the space differencing, are much greater that those to the time differencing. 

NMiMaO, S.14

In mathematics, the Courant–Friedrichs–Lewy (CFL) condition is a necessary condition for convergence while solving certain partial differential equations (usually hyperbolic PDEs) numerically by the method of finite differences. It arises in the numerical analysis of explicit time integration schemes, when these are used for the numerical solution.

If a wave is moving across a discrete spatial grid and we want to compute its amplitude at discrete time steps of equal duration, then this duration must be less than the time for the wave to travel to adjacent grid points.

\(c\leq \frac{\Delta x}{\Delta t}\)

consider the equation \(\frac{du}{dt}=f(u,t)\;\;where\;u=u(t)\)

if a value of f for integration is taken at time level n+1, than this scheme is implicit. 

if this is not the case, than the scheme is implicit. 

second order central:

\((\frac{d^2u}{dx^2 })_j \rightarrow \frac{u_{j+1}-2u_j+u_{j-1}}{\Delta x^2}\)

second order forward:

\((\frac{d^2u}{dx^2})_j \rightarrow \frac{u_{j+2}-2*u_{j+1}+ u_j}{\Delta x^2}\)

second order backward:

\((\frac{d^2u}{dx^2 })_j \rightarrow \frac{u_j-2*u_{j-1}+u_{j-2}}{\Delta x^2 }\)

\(\frac{\partial u}{\partial t}+c*\frac{\partial u}{\partial x}=0\)

\(\frac{u^{n+1}_j-u^n_j}{\Delta t}+c\frac{u^n_{j+1}-u^n_{j-1}}{2\Delta x}=0\)

A solution to a linear equation can be expressed as a Fourier series where each Fourier component is a solution. Thus, we can test the stability with one single Fourier component of the form:

\(.\\u^n_j=u_0*e^{ik(j\Delta x-C_Dn\Delta t)}\)

for the time derivatives: 

\(u^{n+m}{j}=u^n_j*\lambda ^m\)

----------------------------

\(u^{n+1}_j=\lambda u^n_j\;\;\;where\;\;\;\lambda=e^{-ikC_D\Delta t}\)

for stability it's necessary that:

\(|\lambda|=|e^{-ikC_D\Delta t}|\leq1\)

-----------------------

we get:

\(\lambda =1-i\frac{c\Delta t}{\Delta x}*sin(k\Delta x)\)

The absolute value of \(\lambda\) is:

\(|\lambda|^2=1+[\frac{c\Delta t}{\Delta x}*sin(k\Delta x)]^2 \)

\(|\lambda|>1\)

The scheme is unconditionally unstable!

\(\frac{u^{n+1}_j-u^{n-1}_j}{2\Delta t}+c\frac{u^n_{j+1}-u^n_{j-1}}{2\Delta x}=0\)

\(|\lambda|^2= \left(\frac{c\Delta t}{\Delta x}*sin(k\Delta x) \right)^2+1-\frac{c\Delta t}{\Delta x}*sin(k\Delta x) \)

The scheme is conditionally stable with the Courant-Fredrichs-Lewy or CFL criterion! 

A spurious solution to a finite-difference approximation to a differential equation that is not related to the physical solutions of the differential equation.

For instance, the leapfrog differencing scheme can introduce a computational mode.

--------

three level schemes like the leap-frog-scheme requires more than one initial condition (the physical initial condition \(u^{n=0}\) and the computational initial condition \(u^{n=1}\) .

taken the oscillation equation: \(\frac{du}{dt}=i\omega u\;\;where\;\;u=u(t)\)

with the solution:\(u=u_0e^{i\omega t}\)

the leap-frog scheme can be written as: \(u^{n+1}=u^{n-1}+2i\omega \Delta tu^n\)

Then there exists two amplification roots: \(\lambda_{1,2}=i\omega \Delta t\pm \sqrt{1-(\omega \Delta t)^2} \)

Thus there are two solutions of the form: \(u^{n+1}=\lambda u^n \)

Since we are solving a linear equation its solution will be a linear combination of the two solutions:

\(u^n=a\lambda_1^nu^0_1+b\lambda_2^nu^0_2\)

The physical mode is the one where \(\lambda \rightarrow 1\) when \(\Delta t \rightarrow 0\) . The other is the computational mode, which is introduced by the numerical scheme.

The upstream or upwind scheme: 

\(\frac{u_j^{n+1}-u_j^n}{\Delta t}+c\frac{u^n_j-u^n_{j-1}}{\Delta x}\;\;\;if\;c>0\\ \frac{u_j^{n+1}-u_j^n}{\Delta t}+c\frac{u^n_{j+1}-u^n_{j}}{\Delta x}\;\;\;if\;c<0\)

The von Neumann stability analysis gives us:

\(\lambda=1-\frac{c\Delta t}{\Delta x}[1-cos(k\Delta x)+i*sin(k\Delta x)]\)

it follows:

\(0\leq\frac{c\Delta t}{\Delta x}\leq 1\)

which implies that c has to be positive to have stability. The scheme is hence conditionally stable. 

advantage: A disturbance cannot propagate in the direction opposite to the physical advection. Thus, no parasitic waves will contaminate the numerical solution. 

disadvantage: The upstream scheme is highly diffusive!

\(u(x,t)=u_0*exp(i(kx-\omega t))\)

with:

• angular frequency... \(\omega=2\pi f\)  mit  \(f=1/T\)
• wave number... \(k=2\pi/\lambda\)
• Dispersion-relation: \(\omega =kc\)

\(u(x,t)=u_0*exp(ik(x-ct))\)

group velocity: \(c_g=\frac{\partial \omega }{\partial k}\)

Overview: 

1. The leap-frog scheme is not a good idea since the stability analysis of Neumann gives unconditionally instability
2. Instead one should use the Crank- Nicholson scheme
   1. This is unconditional stable but gives meaningful solutions only for \(\gamma \Delta t<1\)
   2. The scheme is strictly seen implicit, but since the discretization doesn't vary in space, one can rearrange the solution to an "explicit scheme".

\(\frac{\partial u}{\partial t}=-\gamma u\) where \(\gamma >0\)

the analytical solution is: \(u(t)=u_0*e^{-\gamma t}\)

leap-frog discretization brings:

\(\frac{u^{n+1}_j-u^{n-1}_j}{2\Delta t}=-\gamma u_j\)

The stability analysis by Neumann shows:

\(\lambda_{1,2}=-\Delta t*\gamma \pm \sqrt{(\gamma \Delta t)^2+1}\)

\(|\lambda |>1 \), which shows unconditional instability!

Typically used for the Rayleigh friction is the Crank-Nickolson scheme: here the right hand side is taken as an average of n-1 and n+1 such that:

\(\frac{u_j^{n+1}-u_j^{n-1}}{2\Delta t}=-\frac{\gamma }{2}(u_j^{n+1}+u_j^{n-1})\)

The stability analysis by Neumann gives:

\(\lambda^2=\frac{1-\gamma\Delta t}{1+\gamma \Delta t}<1\)

Hence, the scheme is absolutely stable. Although, one needs to choose \(\gamma \Delta t<1\) to get a realisitc solution. 

Furthermore, the scheme seems to be implicit, but since the discretization doesn't vary in space one can just rearrange the equation to an explicit scheme! 

Arakawa:

• Euler scheme: conditionally stable
• backward scheme: stable if \(\gamma \Delta t>0\) (unconditional stable)
• trapezoidal scheme: unconditionally stable

Fazit: I would use the Crank-Nicolson scheme, since this is unconditional stable and has an order of accuracy of 2. The popular leap-frog scheme is unconditional unstable and can't be used. 

\(\frac{\partial u}{\partial t}=A\frac{\partial^2u}{\partial x^2}\;\;where\;\;A>0\)

has the analytical solution: \(u(x,t)=u_0e^{\pm ikx-Ak^2t}\)

Euler forward: conditionally stable

Leap-frog: unconditionally unstable

Crank-Nicholson scheme/trapezoidal scheme: unconditionally stable

heat equation with leap frog scheme:

\(\frac{u_j^{n+1}-u_j^{n-1}}{2\Delta t}=A\frac{u_{j+1}-2u_j+u_{j-1}}{(\Delta x)^2}\)

for Neumann method set: \(u^n_j=u_0\lambda^ne^{ikj\Delta x}\)

depends on the time step (n-1,n,n+1) for the spatial derivative on the right hand side: 

• for n: (leap frog)
  • \(\lambda^2 +\frac{8A\Delta t}{(\Delta x)^2}sin^2(\frac{k\Delta x}{2})\lambda-1=0\)
  • the scheme is unconditionally unstable
  • (identity of use: \(sin^2(x)=\frac{1}{2}(1-cos(2x))\))
• for n-1: (Euler forward)
  • the scheme is conditionally stable with \(\frac{A\Delta t}{(\Delta x)^2}<1/4\) (\(-1\leq\lambda\leq 1\))
  • to avoid imaginary roots and oscillating of the solution we recommend: \(\frac{A\Delta t}{(\Delta x)^2}<1/8\) (\(\lambda^2>0\))
  • \(\lambda^2=1-\frac{8A\Delta t}{(\Delta x)^2}sin^2(\frac{k\Delta x}{2})\)
• for an average of n-1 and n+1 (Crank-Nicholson Scheme)
  • absolute stable but imaginary roots which are leading to oscillating solutions can occur
  • \(\lambda^2=\frac{1-\frac{4A\Delta t}{(\Delta x)^2}sin^2(\frac{k\Delta x}{2})}{1+\frac{4A\Delta t}{(\Delta x)^2}sin^2(\frac{k\Delta x}{2})}\)
  • to avoid oscillating imaginary roots on should use: \(\frac{A\Delta t}{(\Delta x)^2}<1/4\)

There leap-frog can be replaced by an Euler forward in time (replace all n-1 by n). The stability analysis doesn't change. 

Fazit: The leap-frog scheme is unconditionally unstable. Although the Crank-Nickolson scheme is unconditionally stable but one should use \(\frac{A\Delta t}{(\Delta x)^2}<1/4\) to avoid oscillating solutions. 

Poisson's equation in two dimensions: 

\(\nabla^2 u=(\frac{\partial^2}{\partial x^2}+\frac{\partial ^2}{\partial y^2})u=f(x,y)\)

If \(f(x,y)=0\) then it is the Laplace equation. 

The equation can be discretised as 

\(\frac{u_{i-1,j}-2u_{i,j}+u_{i+1,j}}{(\Delta x)^2}+\frac{u_{i,j-1}-2u_{i,j}+u_{i,j+1}}{(\Delta y)^2 }=f_{i,j}\)

If we consider a square grid \((\Delta x=\Delta y)\) we get:

\(u_{i,j}=\frac{1}{4}[u_{i-1,j}+u_{i+1,j}+u_{i,j-1}+u_{i,j+1}-(\Delta x)^2f_{i,j}]\)

To solve this we need boundary values for the domain and initial values at the iteration level \(u_{i,j}^m\;\;(m=0)\)

most simple form to solve this by iteration is the Jacobi iteration. Quite inefficient and not used for practical purposes. 

\(u_{i,j}^{m+1}=\frac{1}{4}[u_{i-1,j}^m+u_{i+1,j}^m+u_{i,j-1}^m+u_{i,j+1}^m-(\Delta x)^2f_{i,j}]\)

An improvement in efficiency is the Gauss-Seidel iteration where newly computed values are partly used in the same iteration level: iteration level m+1 values are available for nodes (i-1,j) and (i,j-1).

\(u_{i,j}^{m+1}=\frac{1}{4}[u_{i-1,j}^{m+1}+u_{i+1,j}^m+u_{i,j-1}^{m+1}+u_{i,j+1}^m-(\Delta x)^2f_{i,j}]\)

Gauss-Seidel iteration can be further improved by increasing the convergence rate using the method of SOR (Successive Over Relaxation). The change between two successive Gauss-Seidel iterations is called the residual c, which is defined as: 

\(c=u_{i,j}^{m+1}-u_{i,j}^m\)

The Gauss-Seidel residual is multiplied by a relaxation factor \(\omega\) and new iteration value is obtained from: 

\(u_{i,j}^{m+1}=u_{i,j}^m+\omega c=(1-\omega )u^m_{i,j}+\omega\widehat{u}_{i,j}^{m+1}\)

where  \(\widehat{u}_{i,j}^{m+1}\) denotes the new iteration value obtained from Gauss-Seidel method. 

\(u_{i,j}^{m+1}=(1-\omega)u_{i,j}^m+\)