为什么我的粘性波动方程在 python 模拟中会爆炸？
Why does my viscous wave equation blow up in the python simulation?
我正在尝试根据 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() 的主程序。然后 simulate() 调用 firstIteration() 尝试处理边界条件。然后剩下的由 iterate() 处理。我基本上使用中心差分、前差分和后差分来计算波动方程的导数。
主程序只是在不同的时间步长打印出整个图形。
您为线性（无 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
并执行一个时间步长。如果这导致 P 中的值 >1，您的特征值可能太大。 （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)
我正在尝试根据 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() 的主程序。然后 simulate() 调用 firstIteration() 尝试处理边界条件。然后剩下的由 iterate() 处理。我基本上使用中心差分、前差分和后差分来计算波动方程的导数。
主程序只是在不同的时间步长打印出整个图形。
您为线性（无 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
并执行一个时间步长。如果这导致 P 中的值 >1，您的特征值可能太大。 （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)