一个简单的二维波动方程有限差分正演示例（均匀介质，无吸收边界条件及震源）

从FWI.yml文件中导入环境配置

from IPython.display import Latex
from IPython.display import HTML

二维波动方程格式

$u_{tt}(x,y,t)=u_{xx}(x,y,t)+u_{yy}(x,y,t),(x,y)\in \Omega =(0,1)\times (0,1),t\in (0,1.4)(1)$

$u(x,y,0)=\sin (\pi x)\sin (\pi y),u_t(x,y,0)=0,(x,y)\in (0,1)\times (0,1)(2)$

$u(x,y,t)=0,(x,y)\in \partial \Omega t\in [0,1.4](3)$

解析解为$u(x,y,t)=\cos (\sqrt{2}\pi t)\sin (\pi x)\sin (\pi y)$

算法原理

网格划分取$h=△x=△y=0.01\tau =△t=0.1h$ 由此，空间采样点$N=\frac{1}{h}=100$, 时间采样点$M=\frac{1.4}{0.1h}=1400$,令$x_i=ih,y_j=jh,t_k=k\tau ,i\in [0,N],j\in [0,N],k\in [0,M]$

1中（1）、（2）、（3）可离散为

$u_{i,j}^{k+1}=r^2(u_{i+1,j}^k+u_{i-1,j}^k+u_{i,j+1}^k+u_{i,j-1}^k)+(2-4r^2)u_{i,j}^k-u_{i,j}^{k-1}i\in [1,N-1],j\in [1,N-1],k\in [1,M-1]$,$r=\frac{\tau }{h}=0.1$

$u_{i,j}^0=\sin (\pi x_i)\sin (\pi y_j)=\sin (\pi ih)\sin (\pi jh),i\in [0,N],j\in [0,N]$

$u_{i,j}^1=\frac{r^2}{2}(u_{i+1,j}^0+u_{i-1,j}^0+u_{i,j+1}^0+u_{i,j-1}^0)+(1-2r^2)u_{i,j}^{0}i\in （0,N）,j\in （0,N）$

$u_{0,0}^k=u_{0,N}^k=u_{N,0}^k=u_{N,N}^k=0k\in [0,M]$

具体推导如下

离散化波动方程

对波动方程(1)建立二阶有限差分格式可以得到：

$\frac{u_{i,j}^{k+1}-2u_{i,j}^k+u_{i,j}^{k-1}}{{\tau }^2}=\frac{u_{i+1,j}^k-2u_{i,j}^k+u_{i-1,j}^k}{h^2}+\frac{u_{i,j+1}^k-2u_{i,j}^k+u_{i,j-1}^k}{h^2}(4)$

整理可得：

$u_{i,j}^{k+1}=r^2(u_{i+1,j}^k+u_{i-1,j}^k+u_{i,j+1}^k+u_{i,j-1}^k)+(2-4r^2)u_{i,j}^k-u_{i,j}^{k-1}(5)$

其中，$i\in [1,N-1],j\in [1,N-1],k\in [1,M-1]$,$r=\frac{\tau }{h}=0.1$,局部阶段误差为$o(h^2)$

离散化初值条件

考虑初直条件$u(x,y,0)=\sin (\pi x)\sin (\pi y),(x,y)\in (0,1)\times (0,1)$，差分格式为

$u_{i,j}^0=\sin (\pi x_i)\sin (\pi y_j)=\sin (\pi ih)\sin (\pi jh),i\in [0,N],j\in [0,N](6)$

考虑到初直条件$u_t(x,y,0)=0,(x,y)\in (0,1)\times (0,1)$,利用二阶差商近似：

$\frac{u_{i,j}^1-u_{i,j}^{-1}}{2\tau }=0,(x,y)\in (0,1)\times (0,1)(7)$

将k=0代入（5）,则有：

$u_{i,j}^1=r^2(u_{i+1,j}^0+u_{i-1,j}^0+u_{i,j+1}^0+u_{i,j-1}^0)+(2-4r^2)u_{i,j}^0-u_{i,j}^{-1}(8)$

带入（7）中$u_{i,j}^1=u_{i,j}^{-1}$整理得：

$u_{i,j}^1=\frac{r^2}{2}(u_{i+1,j}^0+u_{i-1,j}^0+u_{i,j+1}^0+u_{i,j-1}^0)+(1-2r^2)u_{i,j}^0(9)$

离散化边界条件

根据边界条件（3）可建立以下差分格式;

$u_{0,0}^k=u_{0,N}^k=u_{N,0}^k=u_{N,N}^k=0(10)$

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import animation

#定义计算区域
x0 = 0.0
x1 = 1.0
y0 = 0.0
y1 = 1.0
N = 100
ds = (x1 - x0)/N

t0 = 0.0
t1 = 1.4
M = 1400
dt = (t1 - t0)/M
x = np.linspace(x0, x1, N+1)
y = np.linspace(y0, y1, N+1)
t = np.linspace(t0, t1, M+1)

r = dt/ds;
#初始化计算区域
(X, Y, T) = np.meshgrid(x,y,t)
#解析解
uu = np.multiply(np.multiply(np.cos(np.sqrt(2) * np.pi * T),np.sin(np.pi * X)),np.sin(np.pi * Y))
#有限差分求解
#创建差分数组
u = np.zeros((N+1,N+1,M+1))
#加入初始条件1
for i in range(0,N+1):
    for j in range(0,N+1):
        u[i,j,0] = np.sin(np.pi*i*ds) * np.sin(np.pi*j*ds)
#加入初始条件2
for i in range(1,N):
    for j in range(1,N):
        u[i,j,1] = r**2 / 2 * (u[i+1,j,0] + u[i-1,j,0] + u[i,j+1,0] + u[i,j-1,0]) + (1 - 2*r**2) * u[i,j,0] 
#加入边界条件
for k in range(0,M+1):
    u[0,0,k] = u[0,N,k] = u[N,0,k] = u[N,N,k] = 0
#递推求[1,N-1]时刻波场
for k in range(1,M):
    for i in range(1,N):
        for j in range(1,N):
            u[i,j,k+1] = r**2*(u[i+1,j,k]+u[i-1,j,k]+u[i,j+1,k]+u[i,j-1,k]) + (2-4*r**2)*u[i,j,k] -u[i,j,k-1]

(X, Y) = np.meshgrid(x,y)
t = 100

fig = plt.figure(figsize=(20,16))
ax1 = fig.add_subplot(131,projection='3d')
surf1 = ax1.plot_surface(X, Y, uu[:,:,t], cmap='viridis')
fig.colorbar(surf1)
ax1.set_xlabel('X Label')
ax1.set_ylabel('Y Label')
ax1.set_zlabel('Z Label')
ax1.set_title(f"analysis solution at time {t*dt}")

ax2 = fig.add_subplot(132,projection='3d')
surf2 = ax2.plot_surface(X, Y, u[:,:,t], cmap='viridis')
fig.colorbar(surf2)
ax2.set_xlabel('X Label')
ax2.set_ylabel('Y Label')
ax2.set_zlabel('Z Label')
ax2.set_title(f"FD solution at time {t*dt}")

ax3 = fig.add_subplot(133,projection='3d')
surf3 = ax3.plot_surface(X, Y, uu[:,:,t]-u[:,:,t], cmap='viridis')
fig.colorbar(surf3)
ax3.set_xlabel('X Label')
ax3.set_ylabel('Y Label')
ax3.set_zlabel('Z Label')
ax3.set_title(f"error at time {t*dt}")

plt.show()

生成动图：

# 创建初始数据
data = uu[:, :, 0]

# 创建图形对象
fig, ax = plt.subplots()

# 创建pcolor对象
pc = ax.pcolor(X, Y, data, cmap='rainbow')

# 定义更新函数
def update(i):
    # 生成新的数据
    new_data = u[:, :, i]
    # 更新pcolor对象
    pc.set_array(new_data.ravel())
    ax.set_title(f"{i*dt} s")
    # 返回更新后的pcolor对象
    return pc,
plt.close()
# 创建动画对象
ani = animation.FuncAnimation(fig, update, frames=range(1, M+1), interval=1, blit=True)

# 显示动画
HTML(ani.to_html5_video())

下一目标：实现密度为常数的二维声波方程的有限差分正演
差分格式为：

$\frac{1}{v(x,z{)}^2}\frac{{\partial }^{2}u}{\partial t^2}-(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2})=f(t)\delta （x-x_s）\delta （z-z_s）$

其中$u(x,z,t)$表示压力，$v(x,z)$为介质速度，$f(t)$是震源函数，$(x_s,z_s)$为震源位置。
方程满足初始条件

$u(x,z,t=0)=0,\frac{\partial u}{\partial t}(x,z,t=0)=0,(x,z)\in \partial \Omega$