二维波动方程代码python

以下是一个简单的二维波动方程的 Python 代码：

import numpy as np
import matplotlib.pyplot as plt

# 定义常数
L = 1.0  # 区域长度
T = 1.0  # 模拟时间
c = 1.0  # 波速
dx = 0.01  # 空间步长
dt = 0.0001  # 时间步长

# 计算网格数
nx = int(L / dx) + 1
nt = int(T / dt) + 1

# 初始化网格
u = np.zeros((nx, nx, nt))

# 设置初始条件
u[:, :, 0] = np.exp(-((np.linspace(0, L, nx) - L/2)**2)[:, None] / 0.01)

# 迭代计算
for n in range(1, nt):
    for i in range(1, nx-1):
        for j in range(1, nx-1):
            uxx = (u[i+1,j,n-1] - 2*u[i,j,n-1] + u[i-1,j,n-1]) / dx**2
            uyy = (u[i,j+1,n-1] - 2*u[i,j,n-1] + u[i,j-1,n-1]) / dx**2
            u[i,j,n] = 2*u[i,j,n-1] - u[i,j,n-2] + c**2 * dt**2 * (uxx + uyy)

# 绘制动画
fig = plt.figure()
ax = fig.add_subplot(111)
plt.ion()

for n in range(nt):
    ax.clear()
    ax.imshow(u[:, :, n], cmap='hot', vmin=-0.5, vmax=0.5)
    plt.title('Time: {:.2f}'.format(n*dt))
    plt.pause(0.001)

plt.ioff()
plt.show()

该代码使用了三重循环来迭代计算波动方程的数值解，其中 `uxx` 和 `uyy` 分别表示在 $x$ 和 $y$ 方向上的二阶导数。最后使用 `imshow` 函数绘制了动画，展示了波动的传播过程。