In [1]:
import numpy as np
import matplotlib.pyplot as plt
In [2]:
dx=1
kk = np.arange(0, np.pi, np.pi/100)
dir_fig = '/home2/datahome/qjamet/Figures'

From https://ocw.mit.edu/courses/12-950-atmospheric-and-oceanic-modeling-spring-2004/resources/lec2/

Advection equation

We first consider the advection equation

     $\partial_t \theta + c \partial_x \theta = 0$

which is not subject to staggering issue because only involves 1 variable $\theta$.

Considering waves of the form $\theta=\theta_0 e^{i(kx - \omega t)}$ leads to the following dispersion relation:

     $w=ck$

which are thus non-disperive since both phase ($\frac{\omega}{k}$) and group ($\partial_k \omega$) velocities are equal.

Now, if we discretize this equation with a centered second order difference, we have

     $\partial_t \theta_i + \frac{c}{2\Delta x} \left( \theta_{i+1} - \theta_{i-1} \right)$ = 0

which leads to the dispersion relation

     $\omega_{2^{\text{nd}}} = \frac{c}{\Delta x} sin (k\Delta x)$

The discretized version of the advection equation is then dispersive, even through it will tend to the continuous dispestion relation for $\Delta x \rightarrow 0$

We can then play the same ball game for fourth (and sixth, but does not seem to be correct) order centered difference, leading to (see https://ocw.mit.edu/courses/12-950-atmospheric-and-oceanic-modeling-spring-2004/resources/lec2/ for details of the derivation)

     $\omega_{4^{\text{th}}} = - \frac{c}{6\Delta x} \left( sin(2k\Delta x) - 8 sin(k\Delta x) \right)$

     $\omega_{6^{\text{th}}} = \frac{c}{\Delta x} sin (k\Delta x) \left( 1 + \frac{4}{6}sin(\frac{k \Delta x}{2}) \left( 1+\frac{4}{5}sin(\frac{k\Delta x}{2}) \right) \right)$

In [3]:
c=1
omega_ref = c*kk
omega_2nd = (c/dx)*np.sin(kk*dx)
omega_4th = - (c/(6*dx)) * (np.sin(2*kk*dx)-8*np.sin(kk*dx))
omega_6th =   (c/dx)     * np.sin(kk*dx) * ( 1 + 4/6*np.sin(kk*dx/2) * ( 1 + 4/5*np.sin(kk*dx/2) ) )
In [5]:
p0 = plt.plot(kk*dx/np.pi,omega_ref*dx/c,'k', label='continous')
p1 = plt.plot(kk*dx/np.pi,omega_2nd*dx/c,'r', label=r'2nd order')
p2 = plt.plot(kk*dx/np.pi,omega_4th*dx/c,'b', label=r'4th order')
#p3 = plt.plot(kk*dx/np.pi,omega_6th*dx/c,'g', label=r'6th order')
#
plt.xlabel(r"$\frac{k \Delta x}{\pi}$")
plt.ylabel(r"$\frac{\omega \Delta x}{c}")
plt.grid()
plt.legend()
plt.title(r'Dispersion relations for the advection equation $\partial_t \theta + c \partial_x \theta = 0$')
#
figN = "disp_rel_2nd_4th_centered_scheme_adv_eqn"
plt.savefig(str("%s/%s.png" % (dir_fig, figN)), bbox_inches='tight')

Gravity wave equation

The gravity wave equation

     $\partial_{tt} \eta = gH \partial_{xx} \eta$

with $\eta=e^{i(kx - \omega t)}$ has the dispersion relation

     $w=\pm \sqrt{gH}k$

Sensitivity to discretization

  • 1/ centered second order spatial derivatives on an unstaggered grid.
    The associated dispersion relation is:

     $\omega^2 = \frac{gH}{\Delta x^2}sin^2(k\Delta x)$

  • 2/ centered second order spatial derivatives on an staggered grid.
    The associated dispersion relation is:

     $\omega^2 = \frac{4gH}{\Delta x^2}sin^2(\frac{k\Delta x}{2})$

In [7]:
g=9.81
H=1000
c=np.sqrt(g*H)
omega0 = c*kk
omega1 = np.sqrt((g*H)/(dx**2)*(np.sin(kk*dx))**2)
omega2 = np.sqrt( (4*g*H)/(dx**2) * (np.sin(kk*dx/2)**2) )
#
plt.plot(kk*dx/np.pi,omega0*dx/c,'k', label='continous')
plt.plot(kk*dx/np.pi,omega1*dx/c,'b', label='unstaggered')
plt.plot(kk*dx/np.pi,omega2*dx/c,'r', label='staggered')
#
plt.xlabel(r"$\frac{k \Delta x}{\pi}$")
plt.ylabel(r"$\frac{\omega \Delta x}{c}$")
plt.grid()
plt.legend()
#
figN = "disp_rel_stag_non-stag_grid_gravity_wave_eqn"
plt.savefig(str("%s/%s.png" % (dir_fig, figN)), bbox_inches='tight')
In [ ]: