计算科学 - 使用后向与中心有限差分近似 - 吾爱随笔录

我正在求解简单的二阶波动方程：

\frac{\partial^{2} E}{\partial t^{2}} = c^{2} \frac{\partial^{2} E}{\partial z^{2}}

$\frac {\partial ^2 E}{\partial t^2} = c^2 \frac {\partial ^2 E}{\partial z^2}$

在（以 SI 单位）的域上：

$z = [0,L=5]$ m, s, m/s $t = [0,t_{max} = 5]$ $c = 1$

和边界/初始条件：

E (z = 0) = E (z = L) = 0

$E(z=0) = E(z=L) = 0$

E (t = 0) = - E (t = t_{m a x}) = s i n (\frac{π z}{L})

$E(t=0) = -E(t=t_{max}) = sin(\frac{\pi z}{L})$

我知道解析解，但我试图用数值求解。（数值和解析解与测试精度进行了比较。）我有几个问题：

1）当我通过对时间的二阶导数应用 2 阶中心差分近似来求解波动方程时，代码工作得非常好。尽管当我及时对二阶导数应用更高阶的反向有限差分近似时，我的解决方案会发散。为什么使用反向差分近似比同阶的中心差分近似更差，有什么特别的原因吗？我的印象是使用更高阶的反向差分方法会给我更高的准确性，但它似乎根本不起作用。如果需要，我可以提供代码，尽管它是一个相当基本的实现。

m/s 时，我能够求解上述方程，但如果我使用 m/s，则解会发散。这是有道理的，因为 Courant 数 ( ) 远大于一。，是否有任何可能的有限差分方案可用于数值求解方程 $c = 1$ $c = 3 \times 10^8$ $C_o = c \frac {\Delta t}{\Delta x}$ $C_o >> 1$ ? 我已经阅读了一些关于 BDF 方法来帮助解决刚性方程的内容，但在这种情况下，唯一的问题是需要更精确的网格，对吗？我也听说过光谱方法，但我仍在尝试学习如何实现它们，以及如果波动方程中包含其他（非线性）项，它们是否足够适用/准确。

编辑：

经过进一步测试，即使没有后向差异，输出实际上也略有偏差。代码：

import matplotlib
matplotlib.use('pdf')
import os 
import matplotlib.pyplot as plt
import numpy as np
from tqdm import tqdm
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import matplotlib.pyplot as plt
from matplotlib.colors import Normalize

c = 1.#3.*(10.**8.)

#z space
Ngridz = 400 #number of intervals in z-axis
zmax = L = 5.
dz = zmax/Ngridz 
z = np.linspace(0.0,zmax,Ngridz+1) 

#time
Ngridt = 400 #number of intervals in t-axis
tmax = L/c
dt = tmax/Ngridt
t = np.linspace(0.0,tmax,Ngridt+1)  

#def dt2(X,k,i,z,t,kX,TYPE): 
#    """Approximation for second time derivative""" - BACKWARD DIFFERENCE
#    ht = t[1]-t[0]
#    if TYPE == 'CONJ':
#        if i == 0:
#            kTT = np.conj(X[k,i]-2.*X[k,i-1]+X[k,i-2])/(ht**2.)
#        elif i == 1:
#            kTT = np.conj(X[k,i-1]-2.*X[k,i]-X[k,i+1])/(ht**2.)
#        elif i == 2:
#            kTT = np.conj(-(X[k,i])+2.*(X[k,i-1])-(X[k,i-2]))/(ht**2.)
#        elif i == 3:
#            kTT = np.conj(-2.*(X[k,i])+5.*(X[k,i-1])-
#            4.*(X[k,i-2])+(X[k,i-3]))/(ht**2.)
#        elif i == 4:
#            kTT = np.conj((-35/12.)*(X[k,i])+(26./3.)*(X[k,i-1])-
#            (19./2.)*(X[k,i-2])+(14./3.)*(X[k,i-3])-
#            (11./12.)*(X[k,i-4]))/(ht**2.)
#        elif i == 5:
#            kTT = np.conj((-15./4.)*(X[k,i])+(77./6.)*(X[k,i-1])-
#            (107./6.)*(X[k,i-2])+13.*(X[k,i-3])-
#            (61./12.)*(X[k,i-4])+(5./6.)*(X[k,i-5]))/(ht**2.)
#        else:
#            kTT = np.conj((-203./45.)*(X[k,i])+(87./5.)*(X[k,i-1])-
#            (117./4.)*(X[k,i-2])+(254./9.)*(X[k,i-3])-
#            (33./2.)*(X[k,i-4])+(27./5.)*(X[k,i-5])-
#            (137./180.)*(X[k,i-6]))/(ht**2.)       
#    else: 
#        if i == 0:
#            kTT = (X[k,i]-2.*X[k,i-1]+X[k,i-2])/(ht**2.)
#        elif i == 1:
#            kTT = (X[k,i-1]-2.*X[k,i]-X[k,i+1])/(ht**2.)
#        elif i == 2:
#            kTT = (-(X[k,i])+2.*(X[k,i-1])-(X[k,i-2]))/(ht**2.)
#        elif i == 3:
#            kTT = (-2.*(X[k,i])+5.*(X[k,i-1])-
#            4.*(X[k,i-2])+(X[k,i-3]))/(ht**2.)
#        elif i == 4:
#            kTT = ((-35/12.)*(X[k,i])+(26./3.)*(X[k,i-1])-
#            (19./2.)*(X[k,i-2])+(14./3.)*(X[k,i-3])-
#            (11./12.)*(X[k,i-4]))/(ht**2.)
#        elif i == 5:
#            kTT = ((-15./4.)*(X[k,i])+(77./6.)*(X[k,i-1])-
#            (107./6.)*(X[k,i-2])+13.*(X[k,i-3])-
#            (61./12.)*(X[k,i-4])+(5./6.)*(X[k,i-5]))/(ht**2.)
#        else:
#            kTT = ((-203./45.)*(X[k,i])+(87./5.)*(X[k,i-1])-
#            (117./4.)*(X[k,i-2])+(254./9.)*(X[k,i-3])-
#            (33./2.)*(X[k,i-4])+(27./5.)*(X[k,i-5])-
#            (137./180.)*(X[k,i-6]))/(ht**2.)      
#    return kTT

def dt2(X,k,i,z,t,kX,TYPE): 
    """Approximation for second time derivative""" - CENTRAL DIFFERENCE
    ht = t[1]-t[0]
    if TYPE == 'CONJ':
        if i == 0:
            kTT = np.conj(X[k,i]-2.*X[k,i+1]+X[k,i+2])/(ht**2.)
        else:
            kTT = np.conj(X[k,i-1]-2.*X[k,i]-X[k,i+1])/(ht**2.) 
    else: 
        if i == 0:
            kTT = (X[k,i]-2.*X[k,i+1]+X[k,i+2])/(ht**2.) 
        else: 
            kTT = (X[k,i-1]-2.*X[k,i]+X[k,i+1])/(ht**2.) 
    return kTT

Ep = np.zeros((len(z),len(t)),dtype=np.complex_)  
TEST = np.zeros((len(z),len(t)),dtype=np.complex_)  

progress = tqdm(total=100.) #this provides a progress bar
total = 100.
progressz = (total)/(len(z))

k = 0
while k < (len(z) - 1):
    progress.update(progressz)

    EpM = Ep

    hz = z[k+1]-z[k] #hz is positive

    i = 0
    while i < (len(t) - 1):
        ht = t[i+1] - t[i]

        EpM[0,i] = 0. #at z = 0 
#        EpM[Ngridz-1,:] = 0. 
        EpM[k,0] = np.sin(np.pi*z[k]) 
#        EpM[k+1,0] = np.sin(np.pi*z[k+1])
        EpM[k,Ngridt-1] = -np.sin(np.pi*z[k])    

        EpM[k+1,i] = (-EpM[k,i-1] + 2.*EpM[k,i] + ((hz/(c))**2.)*dt2(EpM,k,i,z,t,0.,'x'))
#((hz/(c*ht))**2.)*(EpM[k,i+1] - 2.*EpM[k,i] + EpM[k,i-1]))


#        print 'k='+str(k)+',i='+str(i)+',dE='+str(EpM[k+1,i]-EpM[k,i])

        EpM[0,i] = 0. #at z = 0 
#        EpM[Ngridz-1,:] = 0. 
        EpM[k,0] = np.sin(np.pi*z[k]) 
#        EpM[k+1,0] = np.sin(np.pi*z[k+1])
        EpM[k,Ngridt-1] = -np.sin(np.pi*z[k])   

        TEST[k,i] = np.sin(np.pi*z[k])*np.cos(np.pi*t[i])

        i = i + 1 

        Ep[k,:] = EpM[k,:]

    k = k + 1

    if k == (len(z)-1):  
        progress.update(progressz) 

Ereal = (Ep).real

newpath = r'test_wave_equation'
if not os.path.exists(newpath):
    os.makedirs(newpath)

plt.figure(1)
fig, ax = plt.subplots(figsize=(20, 20))
plt.subplot(221)
plt.plot(t,Ereal[0,:],'k:',linewidth = 1.5)
plt.ylabel('Normalized E.real')
plt.subplot(222)
plt.plot(t,Ereal[int(Ngridz*0.33),:],'k:',linewidth = 1.5)
plt.subplot(223)
plt.plot(t,Ereal[int(Ngridz*0.66),:],'k:',linewidth = 1.5)
plt.subplot(224)
plt.plot(t,Ereal[int(Ngridz),:],'k:',linewidth = 1.5)
plt.xlabel(r" t (s)")
plt.savefig(str(newpath)+'/E.real_vs_t.pdf')
#plt.show()

plt.figure(2)
fig, ax = plt.subplots(figsize=(20, 20))
plt.subplot(221)
plt.plot(z,Ereal[:,0],'k:',linewidth = 1.5)
plt.ylabel('Normalized E.real')
plt.subplot(222)
plt.plot(z,Ereal[:,int(Ngridt*0.33)],'k:',linewidth = 1.5)
plt.subplot(223)
plt.plot(z,Ereal[:,int(Ngridt*0.66)],'k:',linewidth = 1.5)
plt.subplot(224)
plt.plot(z,Ereal[:,Ngridt],'k:',linewidth = 1.5)
plt.xlabel(r" z (m)")
plt.savefig(str(newpath)+'/E.real_vs_z.pdf') 

plt.figure(3)
fig, ax = plt.subplots(figsize=(20, 20))
plt.subplot(221)
plt.plot(t,TEST[0,:],'k:',linewidth = 1.5)
plt.ylabel('True E')
plt.subplot(222)
plt.plot(t,TEST[int(Ngridz*0.33),:],'k:',linewidth = 1.5)
plt.subplot(223)
plt.plot(t,TEST[int(Ngridz*0.66),:],'k:',linewidth = 1.5)
plt.subplot(224)
plt.plot(t,TEST[int(Ngridz),:],'k:',linewidth = 1.5)
plt.xlabel(r" t (s)")
plt.savefig(str(newpath)+'/E.real_vs_t.pdf')
#plt.show()

plt.figure(4)
fig, ax = plt.subplots(figsize=(20, 20))
plt.subplot(221)
plt.plot(z,TEST[:,0],'k:',linewidth = 1.5)
plt.ylabel('Normalized E.real')
plt.subplot(222)
plt.plot(z,TEST[:,int(Ngridt*0.33)],'k:',linewidth = 1.5)
plt.subplot(223)
plt.plot(z,TEST[:,int(Ngridt*0.66)],'k:',linewidth = 1.5)
plt.subplot(224)
plt.plot(z,TEST[:,Ngridt],'k:',linewidth = 1.5)
plt.xlabel(r" z (m)")
plt.savefig(str(newpath)+'/E.real_vs_z.pdf')

输出：http: //imgur.com/a/GwPCt

后向差分近似的输出：http: //imgur.com/a/66jem

如您所见，我的输出似乎与解析解不匹配，并且这种差异仅在使用中心差分近似时才出现。随着时间和 z 分量的推进，我的解决方案似乎确实收敛到了正确的答案。我已经尝试增加网格点的数量（即减小间隔大小）和调试，但我似乎没有看到任何明显的错误。关于可能导致问题的任何想法？我已经注释掉了用于二阶时间导数的后向差分近似的函数，但是当我取消注释并使用它而不是 dt2 的中心差分函数时，解决方案会发散并且原始问题仍然存在。