用网站做微信公众号,动漫制作和动漫设计的区别,附近网站建设,广州网站制作选哪家波动方程 - 在三维图中动态显示二维波动方程的解就像水面波澜起伏
flyfish
波动方程的求解结果通常不是一个单一的数值#xff0c;而是一个函数或一组函数#xff0c;这些函数描述了波随时间和空间的传播情况。具体来说#xff0c;波动方程的解可以是关于时间和空间变量的…波动方程 - 在三维图中动态显示二维波动方程的解就像水面波澜起伏
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∂2uc2(∂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^{n1}_{i,j} - 2u^{n}_{i,j} u^{n-1}_{i,j}}{\Delta t^2} ∂t2∂2u≈Δt2ui,jn1−2ui,jnui,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}_{i1,j} - 2u^{n}_{i,j} u^{n}_{i-1,j}}{\Delta x^2} ∂x2∂2u≈Δx2ui1,jn−2ui,jnui−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,j1} - 2u^{n}_{i,j} u^{n}_{i,j-1}}{\Delta y^2} ∂y2∂2u≈Δy2ui,j1n−2ui,jnui,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^{n1}_{i,j} - 2u^{n}_{i,j} u^{n-1}_{i,j}}{\Delta t^2} c^2 \left( \frac{u^{n}_{i1,j} - 2u^{n}_{i,j} u^{n}_{i-1,j}}{\Delta x^2} \frac{u^{n}_{i,j1} - 2u^{n}_{i,j} u^{n}_{i,j-1}}{\Delta y^2} \right) Δt2ui,jn1−2ui,jnui,jn−1c2(Δx2ui1,jn−2ui,jnui−1,jnΔy2ui,j1n−2ui,jnui,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^{n1}_{i,j} 2u^{n}_{i,j} - u^{n-1}_{i,j} c^2 \Delta t^2 \left( \frac{u^{n}_{i1,j} - 2u^{n}_{i,j} u^{n}_{i-1,j}}{\Delta x^2} \frac{u^{n}_{i,j1} - 2u^{n}_{i,j} u^{n}_{i,j-1}}{\Delta y^2} \right) ui,jn12ui,jn−ui,jn−1c2Δt2(Δx2ui1,jn−2ui,jnui−1,jnΔy2ui,j1n−2ui,jnui,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, projection3d)
ax.set_zlim(-4, 4)# 初始化图像
surf ax.plot_surface(X, Y, u[0, :, :], cmapviridis)def update(frame):ax.clear()ax.set_zlim(-4, 4)ax.set_title(fTime step: {frame})surf ax.plot_surface(X, Y, u[frame, :, :], cmapviridis)return surf,# 创建动画
ani FuncAnimation(fig, update, framesm, repeatTrue)
ani.save(wave_quation.gif, writerimagemagick)
plt.show()
