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 $$


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 $$


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)

oscilation 50th 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):
    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    

    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        
            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


1 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:



Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.