为什么我的粘性波动方程在python模拟中会爆炸？

Question

我正在尝试从this paper实现公式16：

上面的方程式16是粘性波动方程式，它应该从大开始并随着时间的流逝而消亡，最终逼近0。但是，当我运行模拟时，它似乎爆炸了。如果您看下面的图片（迭代0-4），它似乎工作正常（即，正弦波看上去越来越小，这很好）。但是，在5.6和更高的位置开始爆炸（您可以看到y轴表示压力以数量级增加）。

这里是输出：

这里是使用Finite Difference Method的代码：

import numpy as np
import math
import matplotlib
import matplotlib.pyplot as plt
import time

x_mesh_param = 1000
# function that increments time steps (solves for the next state of p)
# returns 'p2'
# assuming no boundary conditions because we only care about what a hydrophone
# receives.  Don't care about what happens to wave beyond hydrophone.
def iterate(p0,p1,x_mesh,dt,v,c0 ):

    # step 3 on pg. 9

    dx=x_mesh[1]-x_mesh[0]
    p2 = np.zeros(p1.shape) # this is where we store next state of p

    for i in range(1,len(x_mesh)-1):
        p_txx = (p1[i+1]-p0[i+1]-2*(p1[i]-p0[i])+p1[i-1]-p0[i-1])/ (dt* dx**2)
        p_xx  = (p1[i+1]-2*p1[i]+p1[i-1]) / (dx**2)
        p2[i] = (dt)**2 * (4.0/3 * v * p_txx + c0**2*p_xx)+2*p1[i]-p0[i]

    # what happens at p2[0] (begin) and p2[-1] (end)

    # forward difference (NO BOUNDARY conditions)
    p_txx = (p1[2]-p0[2]-2*(p1[1]-p0[1])+p1[0]-p0[0])/(dt*dx**2)
    p_xx = (p1[2]-2*p1[1]+p1[0]) / (dx**2)
    #p2[0] = (dt)**2 * (4.0/3 * v * p_txx + c0**2*p_xx)+2*p1[0]-p0[0]

    # taylor series
    #p2[0]=1/(1-1/dt+1/(2*(dt)**2))* ((1-1/dt+1/((dt)**2))*p2[1]-p2[2]/(2*(dt**2)))

    #p2[0] = p1[1] # NEUMANN 

    # if you comment out lines 34,37,39 you get Dirichlet Conditions

    # backward difference
    p_txx = (p1[-1]-p0[-1]-2*(p1[-2]-p0[-2])+p1[-3]-p0[-3])/(dt*dx**2)
    p_xx = (p1[-1]-2*p1[-2]+p1[-3]) / (dx**2)
    #p2[-1] = (dt)**2 * (4.0/3 * v * p_txx + c0**2*p_xx)+2*p1[-1]-p0[-1]





    # Dirichlet if line 46 commented out
    return p2

def firstIteration(p1,x_mesh,dt,v,c0 ):


    # step 3 on pg. 9

    dx=x_mesh[1]-x_mesh[0]
    p2 = np.zeros(p1.shape) # this is where we store next state of p

    for i in range(1,len(x_mesh)-1):
        p_txx = 0
        p_xx  = (p1[i+1]-2*p1[i]+p1[i-1]) / (dx**2)
        p2[i] = (dt)**2 * (4.0/3 * v * p_txx + c0**2*p_xx)+p1[i] # assuming p1[i]-p0[i]=0 (coming from assumption that d/dt (p(x,0)) =0 (initial cond. of motion of fluid)

    # what happens at p2[0] (begin) and p2[-1] (end)


    # forward difference (NO BOUNDARY conditions)
    p_txx = 0
    p_xx = (p1[2]-2*p1[1]+p1[0]) / (dx**2)
    #p2[0] = (dt)**2 * (4.0/3 * v * p_txx + c0**2*p_xx)+p1[0]

    # taylor series
    #p2[0]=1/(1-1/dt+1/(2*(dt)**2))* ((1-1/dt+1/((dt)**2))*p2[1]-p2[2]/(2*(dt**2)))

    #p2[0] = p1[1] # NEUMANN 

    # if you comment out lines 34,37,39 you get Dirichlet Conditions

    # backward difference
    p_txx = 0
    p_xx = (p1[-1]-2*p1[-2]+p1[-3]) / (dx**2)
    #p2[-1] = (dt)**2 * (4.0/3 * v * p_txx + c0**2*p_xx)+p1[-1]

    return p2

def simulate():
    states = []
    L=100 
    #dt=.02
    dt=0.001
    v = .001/998
    c0 = 1480

    x_mesh = np.linspace(0,L,x_mesh_param)


    state_initial = np.array([math.sin(x/20*math.pi) for x in x_mesh])
    states.append(state_initial)

    states.append(firstIteration(states[0], x_mesh,dt,v,c0 ))

    for i in range(50):
        new_state = iterate(states[-2], states[-1], x_mesh,dt,v,c0)
        states.append(new_state)

    return states



if __name__=="__main__":
    states = simulate()
    L=100
    x_mesh = np.linspace(0,L,x_mesh_param)

    counter=0
    for s in states:
        fig, ax = plt.subplots()
        ax.plot(x_mesh, s)

        ax.set(xlabel='distance (m)', ylabel='pressure (atm)',
               title='Pressure vs. distance, iteration = '+str(counter))
        ax.grid()

        plt.show()
        print counter
        counter = counter + 1

您将看到有一个主程序调用simulate（）。然后Simulation（）调用firstIteration（）尝试处理边界条件。然后其余的由iterate（）处理。我基本上是使用中央，向前和向后的差异来计算波动方程的导数。

主程序只是在不同的时间步长上打印出整个图形。

Answer 1

您为线性（无p的乘积项）偏微分方程实现了一个数值求解器，该方程似乎在多个时间步中都起作用，然后以快速振荡爆炸。

从数学上讲，您的离散方程组具有随着时间呈指数增长的解。浮点数的最低有效数字的数字噪声将呈指数增长，直到它接管了您要寻找的解决方案。

如果P是由当前和前两个时间步长的N压力值组成的向量（3N值的数组），则代码中所有值的加法和减法都可以用矩阵形式重写（表示为一步）：

P := D @ P

其中D是形状（3N，3N）的矩阵。如果P0是系统的初始状态，则在n时间步长之后的系统状态为[]
Pn = np.linalg.matrix_power(D, n) @ P0
这里，matrix_power是矩阵求幂，即matrix_power(D, 3) == D @ D @ D。如果D的特征值的实部> 1，则解将呈指数增长。您可以使用np.linalg.eigvals检查特征值，最好使用较小的矩阵，例如N=100。如果幸运的话，减小步长dx和dt可能足以使特征值进入非正值范围。
如果不是，则必须从D中消除较大的特征值：
import numpy as np
evals, evecs = np.linalg.eig(D)
evals_1 = np.where(evals.real > 1, 0, evals)
D_1 = evecs @ np.diag(evals_1) @ np.linalg.inv(evecs)
例如，使用小的矩阵（特征值1.01和0.99）：
D = np.array([[ 1.  , -0.01],
              [-0.01,  1.  ]])

# result from eigenvalue elimination as above
D_1 = np.array([[0.495, 0.495],
                [0.495, 0.495]])
尝试使用与特征值0.99的特征向量相差很小的初始向量：
>>> matpow = np.linalg.matrix_power

>>> matpow(D, 200) @ [1, 1.001]
array([0.13038866, 0.13770467])

>>> matpow(D_1, 200) @ [1, 1.001]
array([0.13404666, 0.13404666])

>>> matpow(D, 2000) @ [1, 1.001]
array([-219643.102525,  219643.102525])

>>> matpow(D_1, 2000) @ [1, 1.001]
array([1.86468848e-09, 1.86468848e-09])
在经过200个时间步长之后，数字是相似的，但是在经过2000个时间步长之后，原始矩阵的输出就发散了。
[如果遇到查找特征值和特征向量的麻烦（对于大型矩阵来说可能很慢，则可以通过用matpow(D_1, n)替换为]来大大加快后续计算的速度）>
Pn = evecs @ np.diag(evals_1**n) @ np.linalg.inv(evecs) @ P0