波动方程 - 在三维图中动态显示二维波动方程的解就像水面波澜起伏

flyfish

波动方程的求解结果通常不是一个单一的数值，而是一个函数或一组函数，这些函数描述了波随时间和空间的传播情况。具体来说，波动方程的解可以是关于时间和空间变量的函数，表示在这些变量下波的振幅或位置。

二维波动方程

二维波动方程的形式为：
$$\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)$$其中，$u(x,y,t)$ 是时间 $t$ 和空间 $(x,y)$ 上的波函数，$c$ 是波速。

差分公式

为了数值求解这个偏微分方程，我们可以将其离散化。离散化的过程中使用以下有限差分近似：

1. 时间方向的二阶偏导数：
   $$\frac{{\partial }^{2}u}{\partial t^2}\approx \frac{u_{i,j}^{n+1}-2u_{i,j}^n+u_{i,j}^{n-1}}{\Delta t^2}$$

2. 空间 $x$ 方向的二阶偏导数：
   $$\frac{{\partial }^{2}u}{\partial x^2}\approx \frac{u_{i+1,j}^n-2u_{i,j}^n+u_{i-1,j}^n}{\Delta x^2}$$

3. 空间 $y$ 方向的二阶偏导数：
   $$\frac{{\partial }^{2}u}{\partial y^2}\approx \frac{u_{i,j+1}^n-2u_{i,j}^n+u_{i,j-1}^n}{\Delta y^2}$$

将这些差分公式代入波动方程中，得到离散化的波动方程：
$$\frac{u_{i,j}^{n+1}-2u_{i,j}^n+u_{i,j}^{n-1}}{\Delta t^2}=c^2\left(\frac{u_{i+1,j}^n-2u_{i,j}^n+u_{i-1,j}^n}{\Delta x^2}+\frac{u_{i,j+1}^n-2u_{i,j}^n+u_{i,j-1}^n}{\Delta y^2}\right)$$
整理得到显式更新公式：
$$u_{i,j}^{n+1}=2u_{i,j}^n-u_{i,j}^{n-1}+c^2\Delta t^2\left(\frac{u_{i+1,j}^n-2u_{i,j}^n+u_{i-1,j}^n}{\Delta x^2}+\frac{u_{i,j+1}^n-2u_{i,j}^n+u_{i,j-1}^n}{\Delta y^2}\right)$$

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()