波动方程 - 在三维图中动态显示二维波动方程的解就像水面波澜起伏
波动方程 - 在三维图中动态显示二维波动方程的解就像水面波澜起伏
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()