波动方程 - 在三维图中动态显示二维波动方程的解就像水面波澜起伏
flyfish
波动方程的求解结果通常不是一个单一的数值，而是一个函数或一组函数，这些函数描述了波随时间和空间的传播情况。具体来说，波动方程的解可以是关于时间和空间变量的函数，表示在这些变量下波的振幅或位置。
二维波动方程
二维波动方程的形式为：
∂
2
u
∂
t
2
=
c
2
(
∂
2
u
∂
x
2
+
∂
2
u
∂
y
2
)
\frac{\partial^2 u}{\partial t^2} = c^2 \left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \right)
∂t2∂2u=c2(∂x2∂2u+∂y2∂2u)其中，
u
(
x
,
y
,
t
)
u(x, y, t)
u(x,y,t) 是时间
t
t
t 和空间
(
x
,
y
)
(x, y)
(x,y) 上的波函数，
c
c
c 是波速。
差分公式
为了数值求解这个偏微分方程，我们可以将其离散化。离散化的过程中使用以下有限差分近似：
-
时间方向的二阶偏导数：
∂ 2 u ∂ t 2 ≈ u i , j n + 1 − 2 u i , j n + u i , j n − 1 Δ t 2 \frac{\partial^2 u}{\partial t^2} \approx \frac{u^{n+1}_{i,j} - 2u^{n}_{i,j} + u^{n-1}_{i,j}}{\Delta t^2} ∂t2∂2u≈Δt2ui,jn+1−2ui,jn+ui,jn−1 -
空间 x x x 方向的二阶偏导数：
∂ 2 u ∂ x 2 ≈ u i + 1 , j n − 2 u i , j n + u i − 1 , j n Δ x 2 \frac{\partial^2 u}{\partial x^2} \approx \frac{u^{n}_{i+1,j} - 2u^{n}_{i,j} + u^{n}_{i-1,j}}{\Delta x^2} ∂x2∂2u≈Δx2ui+1,jn−2ui,jn+ui−1,jn -
空间 y y y 方向的二阶偏导数：
∂ 2 u ∂ y 2 ≈ u i , j + 1 n − 2 u i , j n + u i , j − 1 n Δ y 2 \frac{\partial^2 u}{\partial y^2} \approx \frac{u^{n}_{i,j+1} - 2u^{n}_{i,j} + u^{n}_{i,j-1}}{\Delta y^2} ∂y2∂2u≈Δy2ui,j+1n−2ui,jn+ui,j−1n
将这些差分公式代入波动方程中，得到离散化的波动方程：
u
i
,
j
n
+
1
−
2
u
i
,
j
n
+
u
i
,
j
n
−
1
Δ
t
2
=
c
2
(
u
i
+
1
,
j
n
−
2
u
i
,
j
n
+
u
i
−
1
,
j
n
Δ
x
2
+
u
i
,
j
+
1
n
−
2
u
i
,
j
n
+
u
i
,
j
−
1
n
Δ
y
2
)
\frac{u^{n+1}_{i,j} - 2u^{n}_{i,j} + u^{n-1}_{i,j}}{\Delta t^2} = c^2 \left( \frac{u^{n}_{i+1,j} - 2u^{n}_{i,j} + u^{n}_{i-1,j}}{\Delta x^2} + \frac{u^{n}_{i,j+1} - 2u^{n}_{i,j} + u^{n}_{i,j-1}}{\Delta y^2} \right)
Δt2ui,jn+1−2ui,jn+ui,jn−1=c2(Δx2ui+1,jn−2ui,jn+ui−1,jn+Δy2ui,j+1n−2ui,jn+ui,j−1n)
整理得到显式更新公式：
u
i
,
j
n
+
1
=
2
u
i
,
j
n
−
u
i
,
j
n
−
1
+
c
2
Δ
t
2
(
u
i
+
1
,
j
n
−
2
u
i
,
j
n
+
u
i
−
1
,
j
n
Δ
x
2
+
u
i
,
j
+
1
n
−
2
u
i
,
j
n
+
u
i
,
j
−
1
n
Δ
y
2
)
u^{n+1}_{i,j} = 2u^{n}_{i,j} - u^{n-1}_{i,j} + c^2 \Delta t^2 \left( \frac{u^{n}_{i+1,j} - 2u^{n}_{i,j} + u^{n}_{i-1,j}}{\Delta x^2} + \frac{u^{n}_{i,j+1} - 2u^{n}_{i,j} + u^{n}_{i,j-1}}{\Delta y^2} \right)
ui,jn+1=2ui,jn−ui,jn−1+c2Δt2(Δx2ui+1,jn−2ui,jn+ui−1,jn+Δy2ui,j+1n−2ui,jn+ui,j−1n)
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.animation import FuncAnimation
# 定义参数
t_max = 3 # 时间范围
x_max = 1 # 空间范围
y_max = 1
m = 320 # 时间方向的格子数
n = 32 # 空间 x 方向的格子数
k = 32 # 空间 y 方向的格子数
# 时间和空间步长
dt = t_max / (m - 1)
dx = x_max / (n - 1)
dy = y_max / (k - 1)
# 初始化解数组
u = np.zeros((m, n, k))
# 设置初始条件
x_vals, y_vals = np.linspace(0, x_max, n), np.linspace(0, y_max, k)
X, Y = np.meshgrid(x_vals, y_vals)
u[0, :, :] = np.sin(4 * np.pi * X) + np.cos(4 * np.pi * Y)
u[1, :, :] = np.sin(4 * np.pi * X) + np.cos(4 * np.pi * Y)
# 按照公式进行差分
for ii in range(1, m - 1):
for jj in range(1, n - 1):
for kk in range(1, k - 1):
u[ii + 1, jj, kk] = (dt**2 * (u[ii, jj + 1, kk] + u[ii, jj - 1, kk] - 2 * u[ii, jj, kk]) / dx**2 +
dt**2 * (u[ii, jj, kk + 1] + u[ii, jj, kk - 1] - 2 * u[ii, jj, kk]) / dy**2 +
2 * u[ii, jj, kk] - u[ii - 1, jj, kk])
# 创建图形
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.set_zlim(-4, 4)
# 初始化图像
surf = ax.plot_surface(X, Y, u[0, :, :], cmap='viridis')
def update(frame):
ax.clear()
ax.set_zlim(-4, 4)
ax.set_title(f"Time step: {frame}")
surf = ax.plot_surface(X, Y, u[frame, :, :], cmap='viridis')
return surf,
# 创建动画
ani = FuncAnimation(fig, update, frames=m, repeat=True)
ani.save('wave_quation.gif', writer='imagemagick')
plt.show()
