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 appreaciatedimport numpy as npimport sysdef 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 = 20Dt = 0.001Ds = 10numberiter = 200Source = Triangle(90, 0.001)Uprevious = np.zeros([Nz, Nx]) # previous timeUcurrent = np.zeros([Nz, Nx]) # current timeUfuture = np.zeros([Nz, Nx]) # future timeV = np.zeros([Nz, Nx]) V[:][:] = 2000.0# additional not needed Simulation  = np.zeros([numberiter, Nz, Nx])# source activation, center of gridUprevious[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
Explicit 4th order space wave equation not stable implementation?
pde;finite difference;python;discretization;numpy
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 wrongThe 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 stackUprevious = 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