Explicit 4th order space wave equation not stable implementation?

Question

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

# additional not needed 
Simulation  = np.zeros([numberiter, Nz, Nx])

# source activation, center of grid
Uprevious[10][10] = Source[0]

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

    # tringular source position center grid
    if(i < np.size(Source)):
        Ucurrent[10][10] = 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

Answer 1

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:

http://www.youtube.com/watch?v=MA7Rwcq3z9E