二維標量波方程的顯式四階離散化由下式給出：

\ {開始} eqnarray的 U_ {jk} ^ {n + 1} = \ left（\ frac {\ Delta t V_ {jk}} {\ Delta s} \ right）^ 2 \ left（\ sum_ {a = -N} ^ N w_a U_ {j + ak} ^ 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 + ak} ^ n + U_ {j k + a} ^ n \ right）+ 2 U_ {jk} ^ {n} - U_ {jk} ^ {n-1} \ label {1} \ {端} eqnarray的

其中$ \ Delta s = \ Delta x = \ Delta z $

根據Jing-Bo Cheen的說法 1 穩定性限制由下式給出，其中$ V $是$ \ max（V_ {jk}）$ $ N = 2 $空間衍生物的第四個訂單：

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

其中$ w_a $是居中的差異權重，在下面的代碼中定義為：

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

對於下面的代碼，我使用numpy編寫，我有以下內容：

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

網格尺寸：

$$ Nz = Nx = 20 $$

資源：

三角小波放置在（10,10）處，具有23個樣本或（23ms）。 峰值幅度為1.0

叠代次數（足以看到一些傳播）：

$$ numberiter = 200 $$

然後我會期望：

$$ \ Delta t \ leq 0.0030618621784789723 $$

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

$$ \ Delta t = 0.001 $$

更新：。

閱讀了一些帖子我猜是與沒有光滑的巫術或使用空間衍生物高階相關的事情。我仍然沒有足夠的技能來使這裏工作。

我所看到的是很多分散：在6次叠代中，我網格上的最大值（非零）是E-14順序。

同時查看我的網格中的特定點，靠近源（13,13）隨著時間的推移（前50次叠代）我會得到以下圖片

我做錯了什麽？使其穩定的建議非常適用

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