一千萬個為什麽

搜索

明確的四階空間波動方程不穩定實現?

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

\ {開始} 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次叠代)我會得到以下圖片

oscilation 50th iterations

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

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

最佳答案

超過2個月後,沒有人找到我的問題的答案。我發現自己的錯誤和解決方案。為了幫助社區我分享它。

  • 源詞是錯誤的

使用源項$ S(t)$的正確離散方程由下式給出:

\ {開始} 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 + ak} ^ 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} \ {端} multline

不是在開頭給出的等式中沒有考慮源項。代碼沒有實現(上面)它只設置$ U_ {10,10} ^ n = S(t)$。那是不對的。

  • Python類/指針錯誤(最大錯誤

當您執行時間堆棧更新(下圖)時,您不會復制numpy數組內容。你只是改變了類指針。

# 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 

那麽拉普拉斯循環正在做什麽:

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

實際上是:

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

內容副本將是:

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

運行上面的代碼我們得到,簡單的動畫:

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

轉載註明原文: 明確的四階空間波動方程不穩定實現?