# Explicit 4th order space wave equation not stable implementation?

The explicit 4th order discretization for the 2D scalar wave equation is given by:

\begin{eqnarray} U_{jk}^{n+1} = \left( \frac{\Delta t V_{jk} }{\Delta s} \right) ^2 \left( \sum_{a=-N}^N w_a U_{j+a k}^n + \sum_{a=-N}^N w_a U_{j k+a}^n \right) + 2 U_{jk}^{n} - U_{jk}^{n-1} \nonumber \\ U_{jk}^{n+1} = \left( \frac{\Delta t V_{jk}}{\Delta s} \right) ^2 \sum_{a=-N}^N w_a \left( U_{j+a k}^n + U_{j k+a}^n \right) + 2 U_{jk}^{n} - U_{jk}^{n-1} \label{1} \end{eqnarray}

Where $\Delta s = \Delta x = \Delta z$

According to Jing-Bo Cheen1 the stability limit is given by, where $V$ is $\max(V_{jk})$ for $N = 2$ forth order for space derivatives:

$$\Delta t \leq \frac{2 \Delta s}{ V \sqrt{\sum_{a=-N}^{N} (|w_a^1| + |w_a^2|)}}$$

Where $w_a$ are the centered differences weights, In the code bellow it's defined be:

$$w_a = [-1.0/12, 4.0/3, -5.0/2, 4.0/3, -1.0/12]$$

For the code bellow, I wrote using numpy, I had the following:

$$\Delta s = 10$$ $$V = 2000.0$$

Grid dimensions:

$$Nz = Nx = 20$$

Source:

Triangular Wavelet placed at (10, 10) with 23 samples or (23 ms). Amplitude of peak 1.0

Number of iterations (enough to see some propagation):

$$numberiter = 200$$

Then I would expect that :

$$\Delta t \leq 0.0030618621784789723$$

This case, should be ok no? since it's :

$$\Delta t = 0.001$$

Updated:.

After reading some posts I guess is something related to not having an smooth sorce or using spatial derivatives high order. Still I don't have all the skills to make this here work.

What I am seing is a lot of dispersion : just in 6 iterations my maximum values on my grid (non zero) are on E-14 order.

Also looking on a specific point of my grid, close to the source (13,13) I get the following picture as time evolves (first 50 iterations) What am I doing wrong? Suggetions to make it stable will be very appreaciated

import numpy as np
import sys

def Triangle(fc, dt, n=None):
r"""
Triangle Wave one Period.
Defined by frequency and sample rate or by size

* fc        : maximum desired frequency
* dt        : sample rate
* n         : half length of triangle
"""
if(n==None):
n=int(1/float(fc*dt))

t = np.arange(0+1.0/n, 1, 1.0/n)
y = 1-t
y = np.append(y, 0.0)
y_ = 1-t[::-1]
y_ = np.insert(y_, 0, 0.0)

return np.append(y_, np.append(1, y))

Nz = Nx = 20
Dt = 0.001
Ds = 10
numberiter = 200

Source = Triangle(90, 0.001)
Uprevious = np.zeros([Nz, Nx]) # previous time
Ucurrent = np.zeros([Nz, Nx]) # current time
Ufuture = np.zeros([Nz, Nx]) # future time
V = np.zeros([Nz, Nx])
V[:][:] = 2000.0

Simulation  = np.zeros([numberiter, Nz, Nx])

# source activation, center of grid
Uprevious = Source

for i in range(1, numberiter+1):

# tringular source position center grid
if(i < np.size(Source)):
Ucurrent = Source[i]

for k in range(Nz):
for j in range(Nx):
# u0k u1k*ujk*u3k u4k
# Boundary fixed 0 outside
u0k=u1k=u3k=u4k=0.0
uj0=uj1=uj3=uj4=0.0
ujk = Ucurrent[k][j]

if(j-2 > -1):
u0k = Ucurrent[k][j-2]
if(j-1 > -1):
u1k = Ucurrent[k][j-1]
if(j+1 < Nx):
u3k = Ucurrent[k][j+1]
if(j+2 < Nx):
u4k = Ucurrent[k][j+2]
if(k-2 > -1):
uj0 = Ucurrent[k-2][j]
if(k-1 > -1):
uj1 = Ucurrent[k-1][j]
if(k+1 < Nz):
uj3 = Ucurrent[k+1][j]
if(k+2 < Nz):
uj4 = Ucurrent[k+2][j]

d2u_dx2 = (-u0k+16*u1k-30*ujk+16*u3k-u4k)/12.0
d2u_dz2 = (-uj0+16*uj1-30*ujk+16*uj3-uj4)/12.0
Ufuture[k][j] = (d2u_dx2+d2u_dz2)*(Dt*V[k][j]/Ds)**2
Ufuture[k][j] += 2*Ucurrent[k][j]-Uprevious[k][j]

# make the update in the time stack
Uprevious = Ucurrent
Ucurrent = Ufuture
Simulation[i-1] = Ucurrent

sys.stdout.write("\r %d" %(i) )


1 A stability formula for Lax-Wendroff methods with fourth-order in time and general-order in space for the scalar wave equation - Geophysics Vol. 76 No. 2 2011

After more than 2 months and no one found the answer to my problem. I found myself the errors and the solution. To help the community I'am sharing it.

• The source term is wrong

The correct discrete equation using a source term $S(t)$ is given by :

\begin{multline} U_{jk}^{n+1} = \left( \frac{\Delta t V_{jk}}{\Delta s} \right) ^2 \left[ \sum_{a=-N}^N w_a \left( U_{j+a k}^n + U_{j k+a}^n \right) \right] + S_{jk}^n \Delta t^2 V_{jk}^2 + 2 U_{jk}^{n} - U_{jk}^{n-1} \label{xb1} \end{multline}

Not by the equation given in the beginning that doesn't take in account the source term. The code doesn't implement (above) it just set's $U_{10,10}^n = S(t)$. That's not right.

• Python class/pointer error (biggest mistake)

When you do the time stack update (bellow), you are not copying the numpy arrays content. You are just changing the class pointers.

# make the update in the time stack
Uprevious = Ucurrent #  Uprevious is pointing to Ucurrent (previous ok).
Ucurrent = Ufuture # Ucurrent is pointing to Ufuture.
Simulation[i-1] = Ucurrent


So what the Laplacian loop is doing:

Ufuture[k][j] += 2*Ucurrent[k][j]-Uprevious[k][j]


is in fact:

Ufuture[k][j] += 2*Ufuture[k][j]-Uprevious[k][j]


Content copy would be:

Ucurrent[:][:] = Ufuture[:][:]


Running the code above we get, simple animation: