从FWI.yml文件中导入环境配置
from IPython.display import Latex
from IPython.display import HTML
二维波动方程格式
u t t ( x , y , t ) = u x x ( x , y , t ) + u y y ( x , y , t ) , ( x , y ) ∈ Ω = ( 0 , 1 ) × ( 0 , 1 ) , t ∈ ( 0 , 1.4 ) ( 1 ) u_{tt}(x,y,t)=u_{xx}(x,y,t)+u_{yy}(x,y,t),\quad (x,y)\in\Omega = (0,1)\times(0,1),t\in(0,1.4)\quad(1) utt(x,y,t)=uxx(x,y,t)+uyy(x,y,t),(x,y)∈Ω=(0,1)×(0,1),t∈(0,1.4)(1)
u ( x , y , 0 ) = sin ( π x ) sin ( π y ) , u t ( x , y , 0 ) = 0 , ( x , y ) ∈ ( 0 , 1 ) × ( 0 , 1 ) ( 2 ) u(x,y,0) = \sin(\pi x)\sin(\pi y),\quad u_t(x,y,0) = 0, \quad(x,y)\in(0,1)\times(0,1)\quad(2) u(x,y,0)=sin(πx)sin(πy),ut(x,y,0)=0,(x,y)∈(0,1)×(0,1)(2)
u ( x , y , t ) = 0 , ( x , y ) ∈ ∂ Ω t ∈ [ 0 , 1.4 ] ( 3 ) u(x,y,t) = 0,\quad(x,y)\in\partial \Omega\quad t\in[0,1.4]\quad(3) u(x,y,t)=0,(x,y)∈∂Ωt∈[0,1.4](3)
解析解为 u ( x , y , t ) = cos ( 2 π t ) sin ( π x ) sin ( π y ) u(x,y,t)=\cos(\sqrt2\pi t)\sin(\pi x)\sin(\pi y) u(x,y,t)=cos(2πt)sin(πx)sin(πy)
算法原理
网格划分取 h = △ x = △ y = 0.01 τ = △ t = 0.1 h h =\triangle x = \triangle y = 0.01\quad \tau = \triangle t = 0.1h h=△x=△y=0.01τ=△t=0.1h 由此,空间采样点 N = 1 h = 100 N=\frac{1}{h}=100 N=h1=100, 时间采样点 M = 1.4 0.1 h = 1400 M = \frac{1.4}{0.1h}= 1400 M=0.1h1.4=1400,令 x i = i h , y j = j h , t k = k τ , i ∈ [ 0 , N ] , j ∈ [ 0 , N ] , k ∈ [ 0 , M ] x_i = ih,y_j =jh,t_k=k\tau,\quad i\in[0,N],j\in[0,N],k\in[0,M] xi=ih,yj=jh,tk=kτ,i∈[0,N],j∈[0,N],k∈[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 − 4 r 2 ) u i , j k − u i , j k − 1 i ∈ [ 1 , N − 1 ] , j ∈ [ 1 , N − 1 ] , k ∈ [ 1 , M − 1 ] 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}\quad i\in[1,N-1],j\in[1,N-1],k\in[1,M-1] ui,jk+1=r2(ui+1,jk+ui−1,jk+ui,j+1k+ui,j−1k)+(2−4r2)ui,jk−ui,jk−1i∈[1,N−1],j∈[1,N−1],k∈[1,M−1], r = τ h = 0.1 r=\frac{\tau}{h}=0.1 r=hτ=0.1
u i , j 0 = sin ( π x i ) sin ( π y j ) = sin ( π i h ) sin ( π j h ) , i ∈ [ 0 , N ] , j ∈ [ 0 , N ] u_{i,j}^0 = \sin(\pi x_i)\sin(\pi y_j)=\sin(\pi ih)\sin(\pi jh), \quad i\in[0,N],j\in[0,N] ui,j0=sin(πxi)sin(πyj)=sin(πih)sin(πjh),i∈[0,N],j∈[0,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 i ∈ ( 0 , N ) , j ∈ ( 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\quad i\in(0,N),j\in(0,N) ui,j1=2r2(ui+1,j0+ui−1,j0+ui,j+10+ui,j−10)+(1−2r2)ui,j0i∈(0,N),j∈(0,N)
u 0 , 0 k = u 0 , N k = u N , 0 k = u N , N k = 0 k ∈ [ 0 , M ] u_{0,0}^k=u_{0,N}^k=u_{N,0}^k=u_{N,N}^k=0 \quad k\in[0,M] u0,0k=u0,Nk=uN,0k=uN,Nk=0k∈[0,M]
具体推导如下
离散化波动方程
对波动方程(1)建立二阶有限差分格式可以得到:
u i , j k + 1 − 2 u i , j k + u i , j k − 1 τ 2 = u i + 1 , j k − 2 u i , j k + u i − 1 , j k h 2 + u i , j + 1 k − 2 u i , j k + u i , j − 1 k h 2 ( 4 ) \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}\quad(4) τ2ui,jk+1−2ui,jk+ui,jk−1=h2ui+1,jk−2ui,jk+ui−1,jk+h2ui,j+1k−2ui,jk+ui,j−1k(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 − 4 r 2 ) u i , j k − u i , j k − 1 ( 5 ) 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}\quad(5) ui,jk+1=r2(ui+1,jk+ui−1,jk+ui,j+1k+ui,j−1k)+(2−4r2)ui,jk−ui,jk−1(5)
其中, i ∈ [ 1 , N − 1 ] , j ∈ [ 1 , N − 1 ] , k ∈ [ 1 , M − 1 ] i\in[1,N-1],j\in[1,N-1],k\in[1,M-1] i∈[1,N−1],j∈[1,N−1],k∈[1,M−1], r = τ h = 0.1 r=\frac{\tau}{h}=0.1 r=hτ=0.1,局部阶段误差为 o ( h 2 ) o(h^2) o(h2)
离散化初值条件
考虑初直条件 u ( x , y , 0 ) = sin ( π x ) sin ( π y ) , ( x , y ) ∈ ( 0 , 1 ) × ( 0 , 1 ) u(x,y,0) = \sin(\pi x)\sin(\pi y),\quad(x,y)\in(0,1)\times(0,1) u(x,y,0)=sin(πx)sin(πy),(x,y)∈(0,1)×(0,1),差分格式为
u i , j 0 = sin ( π x i ) sin ( π y j ) = sin ( π i h ) sin ( π j h ) , i ∈ [ 0 , N ] , j ∈ [ 0 , N ] ( 6 ) u_{i,j}^0 = \sin(\pi x_i)\sin(\pi y_j)=\sin(\pi ih)\sin(\pi jh), \quad i\in[0,N],j\in[0,N]\quad(6) ui,j0=sin(πxi)sin(πyj)=sin(πih)sin(πjh),i∈[0,N],j∈[0,N](6)
考虑到初直条件 u t ( x , y , 0 ) = 0 , ( x , y ) ∈ ( 0 , 1 ) × ( 0 , 1 ) \quad u_t(x,y,0) = 0, \quad(x,y)\in(0,1)\times(0,1) ut(x,y,0)=0,(x,y)∈(0,1)×(0,1),利用二阶差商近似:
u i , j 1 − u i , j − 1 2 τ = 0 , ( x , y ) ∈ ( 0 , 1 ) × ( 0 , 1 ) ( 7 ) \frac{u_{i,j}^1-u_{i,j}^{-1}}{2\tau} =0, \quad(x,y)\in(0,1)\times(0,1)\quad(7) 2τui,j1−ui,j−1=0,(x,y)∈(0,1)×(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 − 4 r 2 ) u i , j 0 − u i , j − 1 ( 8 ) 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}\quad(8) ui,j1=r2(ui+1,j0+ui−1,j0+ui,j+10+ui,j−10)+(2−4r2)ui,j0−ui,j−1(8)
带入(7)中 u i , j 1 = u i , j − 1 u_{i,j}^1=u_{i,j}^{-1} ui,j1=ui,j−1整理得:
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 ( 9 ) 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\quad(9) ui,j1=2r2(ui+1,j0+ui−1,j0+ui,j+10+ui,j−10)+(1−2r2)ui,j0(9)
离散化边界条件
根据边界条件(3)可建立以下差分格式;
u 0 , 0 k = u 0 , N k = u N , 0 k = u N , N k = 0 ( 10 ) u_{0,0}^k=u_{0,N}^k=u_{N,0}^k=u_{N,N}^k=0 \quad(10) u0,0k=u0,Nk=uN,0k=uN,Nk=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())
下一目标:实现密度为常数的二维声波方程的有限差分正演
差分格式为:
1 v ( x , z ) 2 ∂ 2 u ∂ t 2 − ( ∂ 2 u ∂ x 2 + ∂ 2 u ∂ y 2 ) = f ( t ) δ ( x − x s ) δ ( z − z s ) \frac{1}{v(x,z)^2}\frac{\partial^2u}{\partial t^2}-(\frac{\partial^2u}{\partial x^2}+\frac{\partial^2u}{\partial y^2})=f(t)\delta(x-x_s)\delta(z-z_s) v(x,z)21∂t2∂2u−(∂x2∂2u+∂y2∂2u)=f(t)δ(x−xs)δ(z−zs)
其中 u ( x , z , t ) u(x,z,t) u(x,z,t)表示压力, v ( x , z ) v(x,z) v(x,z)为介质速度, f ( t ) f(t) f(t)是震源函数, ( x s , z s ) (x_s,z_s) (xs,zs)为震源位置。
方程满足初始条件
u ( x , z , t = 0 ) = 0 , ∂ u ∂ t ( x , z , t = 0 ) = 0 , ( x , z ) ∈ ∂ Ω u(x,z,t=0)=0,\frac{\partial u}{\partial t}(x,z,t=0)=0,(x,z)\in \partial \Omega u(x,z,t=0)=0,∂t∂u(x,z,t=0)=0,(x,z)∈∂Ω