numerics for geofluid dynamics

finite difference approximation or finite difference scheme. 

The difference between the numerical and the true solution. 

\(u_j^n-u(j\Delta x,n\Delta t)\)

Most of the time we cannot find the error of the solution. However, we can always find a meassure of accuracy. 

The truncation error of a consitent scheme can be made abitrarily small by a suffiecient reduction of the increments \(\Delta x\;and\;\Delta t\).

Although, this cannot be expected for the error of the numerical solution: \(u^n_j-u(j\Delta x, n \Delta t)\)

Definition: If the error of the numerical solution approaches zero as the grid is refined, the solution is called convergent. If this holds for any initial conditions, than the scheme is called convergent too. 

Consistency od a scheme does not guarantee convergence!

For convergence of a scheme, the initial conditions needs to lie within the domain of dependency.

For the advection equations  (Euler forward and upstream) this is the case if the CFL criterion holds. That is:

\(c\Delta t \leq\Delta x \)

(If a scheme is convergent it is also stable!)

consistency plus stability gives convergence!

consistency plus stability gives convergence!

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. 

         

         

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+\)

The linearized shallow water equations or the gravity-inertia wave equations are:

\(\frac{\partial u}{\partial t}-fv=-g\frac{\partial h}{\partial x}\\ \frac{\partial v}{\partial t}+fu=-g\frac{\partial h}{\partial y}\\ \frac{\partial h}{\partial t}+H(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y})=0\)

If we consider solutions in wave form than we obtain:

\(.\\\omega^2=f^2+gH(k^2+l^2)\)

For each of the three grids we use the simplest centred approximations for the space derivative and the Coriolis terms.

If we substitute the wave solution into the numerical schemes we get the frequency equations. This are basically dependent of two paramters: \(k\Delta x\) and \(gH/f^2\) (Rossby Radius to the square)

The equations are not shown but discussed with the following graphic:

Grid A:

1. The frequency reaches a maximum at  \(k\Delta x=\pi/2\) , i.e. a wave length of 4 grids. The group velocity for this wave length is zero, which means the wave energy stays near that point. 
2. For \(\pi/2 < k \Delta x < \pi\)the frequency decreases as the wave number increases. Thus, the group velocity vector has the wrong sign. 
3. For the two-grid-interval wave (\(k\Delta x=\pi\)) the group velocity vector is again zero. 

Grid B:

1. The frequency increases monotonically. 
2. Reaches a maximum for \(k\Delta x=\pi\) . There the group velocity is again zero. 

Grid C: 

1. The frequency increases monotonically in a similiar way to the one for the B-grid, when the Rossby Radius is longer than a half grid (\((\sqrt{(gH)/f}> \frac{\Delta x}{2})\) .
2. If the Rossby-Radius resolution fall beneath that, then the frequency will decrease for increasing wavenumber throughout \(0< k\Delta x <\pi \)
3. The big advantage of the C-grid: Velocities are perpendicular on the walls of the grid-box, which makes differentiating straightforward. 

Taylor series of f(x) at the point a:

\(Tf(x;a)=\sum^{\infty}_{n=0}\frac{f^{(n)}(a)}{n!}(x-a)^n\)

The grid point method

• A numerical scheme needs to be above all consistent.
• This means that the approximation should approach the derivative when the grid interval approaches zero. 
• For a consistent scheme the truncation error needs to be at least of order one.