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

Question

我试图从this paper实现等式16：

上面的方程16是粘性波动方程，它应该开始变大，随着时间的推移逐渐消失，最终接近0。然而，当我运行我的模拟时，它似乎爆炸了。如果您查看下面的图像（迭代0-4），它似乎工作正常（即，正弦波似乎变小了，这很好）。然而，在ierations 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（）的主程序。然后simulate（）调用firstIteration（），该函数尝试处理边界条件。然后迭代（）处理其余部分。我基本上是用中心差分，向前差分和向后差分来计算波动方程的导数

主程序只需在不同的时间步打印出整个图形

Answer 1

您为一个线性（p的无乘积项）偏微分方程实现了一个数值解算器，该方程似乎适用于多个时间步，然后随着快速振荡而爆发

从数学上讲，离散化方程组的解随时间呈指数增长。数字噪声——浮点数的最低有效位将呈指数增长，直到它取代您所寻找的解决方案

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

P := D @ P

其中D是形状为（2N，2N）的矩阵。如果P0是系统的初始状态，那么n时间步之后的系统状态将是

Pn = np.linalg.matrix_power(D, n) @ P0

这里，matrix_power是矩阵求幂，即matrix_power(D, 3) == D @ D @ D。如果D具有绝对值的特征值>；1，那么将出现指数增长的解决方案。您的数据表明，最大特征值约为1000，这意味着其在6个时间步长内增长了因子1e+18，而物理相关特征向量约为0.9。请注意，浮点数的精度约为1e-16

您可以使用np.linalg.eigvals检查特征值，最好使用小矩阵，例如N=100。如果幸运的话，那么减小步长dx和dt就足以使特征值位于1以下的区域

快速脏液

这很快就能实现，但效率不高。您可以尝试从D中消除较大的特征值：

import numpy as np
evals, evecs = np.linalg.eig(D)
evals_1 = np.where(np.abs(evals) > 1, 0, evals)
D_1 = evecs @ np.diag(evals_1) @ np.linalg.inv(evecs)

例如，对于小矩阵（特征值0.9和1000）：

D = np.array([[500.45, -499.55], 
              [-499.55, 500.45]])

# D_1 calculated as above.
D_1 = np.array([[0.45, 0.45],
                [0.45, 0.45]])

P0 = [1, 1+8e-16]
powm = np.linalg.matrix_power

ns = np.arange(7)
Ps_D = np.array([powm(D, n) @ P0 for n in ns])
Ps_D1 = np.array([powm(D_1, n) @ P0 for n in ns])

import matplotlib.pyplot as plt
fig, ax = plt.subplots()
ax.semilogy(ns, Ps_D[:, 1], 'ro-', label='Using D')
ax.semilogy(ns, Ps_D1[:, 1], 'bs-', label='Using D_1')
ax.set_xlabel('Iteration number')
ax.set_ylabel('P[1]')
ax.grid()
ax.legend()
fig.show()

如果您遇到了查找特征值和特征向量的麻烦（对于大型矩阵，这可能会很慢），您可以通过将matpow(D_1, n)替换为

Pn = evecs @ np.diag(evals_1**n) @ np.linalg.inv(evecs) @ P0

快速检查调谐步长是否正常

通过将初始状态设置为，可以判断是否存在大于1的特征值

P = np.zeros(N)
P[N//4] = 1

以及执行一个单一的时间步。如果这导致出现值>；1在P中，您可能具有过大的特征值。（新的P值实际上是矩阵的D[:, N//4]列）

隐式方法

与其像上面描述的那样直接地将大的特征值归零，不如求助于implicit method。这需要写出D矩阵。如果P是当前状态Pnew是下一个时间步的状态，则从方程开始

Pnew = D @ (0.5*(P + Pnew))

求解Pnew的矩阵方程：

# Pnew = inv(1 - 0.5 D) (0.5 D) P
nD = D.shape[0]
D_2 = np.linalg.inv(np.eye(nD) - 0.5*D) @ (0.5*D)

# time step:
Pnew = D_2 @ P

这是Crank-Nicolson method。它要求您构造D矩阵。应用与上述相同的示例给出：

D_2 = np.array([[-0.09191109,  0.91009291],
                [ 0.91009291, -0.09191109]])
evals = np.array[ 0.81818182, -1.00200401]

这将最大特征值从1000减少到-1.002，这更好，但仍然大于1（绝对值），因此从1e-15的舍入误差开始，它将在大约34k个时间步后超过您想要的解。此外，它将0.9特征值减少到0.8。你需要找出哪一个更能描述物理学。不幸的是，数值求解偏微分方程是困难的

编辑：Lutz Lehmann指出，处理大型矩阵的逆矩阵（x中精细网格的结果）不是一个好主意。注意D是一个稀疏矩阵，最好表示为scipy.sparse.csr_matrix、csc_matrix或dia_matrix而不是np.array。然后你就解决了

(1 - 0.5*D) @ P_new = 0.5*D @ P

这是一个矩阵方程，形式为Ax=b，带有x和b向量：

A = scipy.sparse.identity(nD) - 0.5*D

for every time step:
   b = D @ (0.5*P)
   P = scipy.sparse.linalg.spsolve(A, b)

请注意spsolve仍然可能很慢。它可能有助于将P与压力的当前值和先前值交错排列，以便D是带对角的，并且可以存储为这是一个稀疏的dia_matrix

编辑：

因此，最有效的方法是：

import scipy.sparse as sp

n = 200 # Example - number of x points
nD = 2*n # matrix size
D = sp.dok_matrix((nD, nD))

# Initialize pressures
P = np.zeros(nD)

# ... initialize P and nonzero matrix elements here ...

# Prepare matrices
D = sp.csc_matrix(D) # convert to efficient storage format
A = sp.identity(nD) - 0.5*D
A_lu = sp.linalg.splu(A)

for _i in range(num_time_steps):
   b = D @ (0.5*P)
   P = A_lu.solve(b)