二维波动方程数值模拟(非均匀介质)

2024-02-12 23:38

本文主要是介绍二维波动方程数值模拟(非均匀介质),希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

二维波动方程数值模拟(非均匀介质)

0. 前言

最近学习了波动方程在弹性介质中的传播,笔者对其中的面波、反射波、折射波的产生原理产生了兴趣,本文想要就波运动的本质方程 – 波动方程, 直接通过数值模拟的方法, 看看是否在介质边界是否能产生反射波、折射波,并观察其传播特性。

两层波速不同的介质实验结果如下:

在这里插入图片描述

1. 二维波动方程

∂ 2 u ∂ t 2 = v 2 ( ∂ 2 u ∂ x 2 + ∂ 2 u ∂ y 2 ) (1) \frac{\partial^2 u}{\partial t^2} = v^2(\frac{\partial^2u}{\partial x^2} + \frac{\partial^2u}{\partial y^2}) \tag{1} t22u=v2(x22u+y22u)(1)

2. 二维波动方程差分形式建立

2.1 一阶偏导差分形式

对于任意的函数 f f f, 若其在 x x x 邻域内可导,则
∂ f ∂ x = lim ⁡ Δ x → 0 f ( x + Δ x ) − f ( x ) Δ x \dfrac{\partial f}{\partial x} = \lim_{\Delta x \to 0} \frac{f(x+ \Delta x) - f(x)}{\Delta x} xf=Δx0limΔxf(x+Δx)f(x)
x = x k , x + Δ x = x k + 1 x = x_k, \ \ x + \Delta x = x_{k + 1} x=xk,  x+Δx=xk+1, 那么有:
∂ f ∂ x = lim ⁡ Δ x → 0 f ( x k + 1 ) − f ( x k ) Δ x \frac{\partial f}{\partial x} = \lim_{\Delta x \to0}\frac{f(x_{k+1}) - f(x_k)}{\Delta x} xf=Δx0limΔxf(xk+1)f(xk)
Δ x \Delta x Δx 充分小时, 可以近似的认为
∂ f ∂ x = f ( x k + 1 ) − f ( x k ) Δ x \frac{\partial f}{\partial x} = \frac{f(x_{k+1}) - f(x_k)}{\Delta x} xf=Δxf(xk+1)f(xk)

2.2 二阶偏导的差分形式

同理, 对于任意函数 f f f, 若其在 x x x 邻域内二阶可导, 则可近似认为:
∂ 2 f ∂ x 2 = f ′ ( x k + 1 ) − f ′ ( x k ) Δ x = f ( x k + 1 ) − 2 f ( x k ) + f ( x k − 1 ) Δ x 2 (2) \frac{\partial^2 f}{\partial x^2} = \frac{f'(x_{k+1}) - f'(x_k)}{\Delta x}= \frac{f(x_{k+1}) - 2f(x_{k}) + f(x_{k-1})}{\Delta x^2} \tag{2} x22f=Δxf(xk+1)f(xk)=Δx2f(xk+1)2f(xk)+f(xk1)(2)
式中 f ′ f' f f f f x x x 的偏导。

2.3 二维均匀介质波动方程差分形式

对于均匀介质

​ 我们将二维空间 x , y x, y x,y 分别分割为 n , m , l n, m, l n,m,l个等分段,且满足这些等分段张成的小正方体区域足够小, 每一个区域所在的坐标为 x i , y j , t k x_i, y_j, t_k xi,yj,tk, 每小段长度分别为 Δ x , Δ y , Δ t \Delta x, \Delta y, \Delta t Δx,Δy,Δt, 那么可以将波场 u u u 按上式作差分形式展开:

记号 u i , j k u_{i, j}^k ui,jk 表达经数值差分得到的波场 u u u 在网格点 ( x i , y j , t k ) (x_i, y_j, t_k) (xi,yj,tk) 处的值

u i , j k − 1 − 2 u i , j k + u i , j k + 1 Δ t 2 = v 2 [ u i − 1 , j k − 2 u i , j k + u i + 1 , j k Δ x 2 + u i , j − 1 k − 2 u i , j k + u i , j + 1 k Δ y 2 ] \frac{u_{i, j}^{k - 1} - 2{u_{i, j}^k} + u_{i, j}^{k + 1}}{\Delta t^2} = v^2\left[\frac{u_{i - 1, j}^{k } - 2{u_{i, j}^k} + u_{i + 1, j}^{k}}{\Delta x^2} + \frac{u_{i, j - 1}^{k} - 2{u_{i, j}^k} + u_{i, j + 1}^{k}}{\Delta y^2}\right] Δt2ui,jk12ui,jk+ui,jk+1=v2[Δx2ui1,jk2ui,jk+ui+1,jk+Δy2ui,j1k2ui,jk+ui,j+1k]

经化简得:
u i , j k + 1 = r x ( u i − 1 , j k + u i + 1 , j k ) + r y ( u i , j − 1 k , u i , j + 1 k ) + 2 ( 1 − r x − r y ) u i , j k − u i , j k − 1 (3) u_{i, j}^{k + 1} = r_x{(u_{i -1, j}^k + u_{i+1, j}^k)} + r_y(u_{i, j-1}^k, u_{i, j +1}^k) + 2(1 - r_x - r_y) u_{i, j}^{k} - u_{i, j}^{k - 1} \tag{3} ui,jk+1=rx(ui1,jk+ui+1,jk)+ry(ui,j1k,ui,j+1k)+2(1rxry)ui,jkui,jk1(3)
式中
r x = Δ t 2 Δ x 2 v 2 , r y = Δ t 2 Δ y 2 v 2 r_x = \frac{\Delta t^2}{\Delta x^2}v^2, \ r_y = \frac{\Delta t^2}{\Delta y^2}v^2 rx=Δx2Δt2v2, ry=Δy2Δt2v2
根据 ( 3 ) (3) (3) 式, 如果我们有第 k k k 步 和 第 k − 1 k - 1 k1 步的波, 就可以推算出 k + 1 k+1 k+1 步的波场。

2.4 二维非均匀介质波动方程差分形式

对于非均匀介质, 其本质在于各点的波速不同, 观察 ( 3 ) (3) (3) 式, 实际上对于不同的点只有 r x , r y r_x, r_y rx,ry 不同。

我们令 v = ρ ( x , y ) v = \rho(x, y) v=ρ(x,y)

记记号 v i , j v_{i, j} vi,j 表示 ρ ( x i , y j ) \rho(x_i, y_j) ρ(xi,yj)

记号 r x i , j r_{x}^{i, j} rxi,j 表示 Δ t 2 Δ x 2 v i , j 2 \dfrac{\Delta t^2}{\Delta x^2}v_{i, j}^2 Δx2Δt2vi,j2,

记号 r y i , j r_{y}^{i, j} ryi,j 表示 Δ t 2 Δ y 2 v i , j 2 \dfrac{\Delta t^2}{\Delta y^2}v_{i, j}^2 Δy2Δt2vi,j2.

对于 ( 3 ) (3) (3) 式中波长的每一个点 u i , j u_{i, j} ui,j, 将 r x , r y r_x, r_y rx,ry 改写为 r x i , j , r y i , j r_x^{i, j}, r_y^{i, j} rxi,j,ryi,j

得到:
u i , j k + 1 = r x i − 1 , j ⋅ u i − 1 , j k + r x i + 1 , j ⋅ u i + 1 , j k + r y i , j − 1 ⋅ u i , j − 1 k + r y i , j + 1 ⋅ u i , j + 1 k + 2 ( 1 − r x i , j − r y i , j ) u i , j k − u i , j k − 1 (4) u_{i, j}^{k + 1} = r_x^{i-1, j}\cdot u_{i -1, j}^k + r_x^{i+1, j}\cdot u_{i+1, j}^k \\ \ + r_y^{i, j-1}\cdot u_{i, j-1}^k + r_y^{i, j+1}\cdot u_{i, j +1}^k\\ \ + 2(1 - r_x^{i, j} - r_y^{i, j}) u_{i, j}^{k} - u_{i, j}^{k - 1} \tag{4} ui,jk+1=rxi1,jui1,jk+rxi+1,jui+1,jk +ryi,j1ui,j1k+ryi,j+1ui,j+1k +2(1rxi,jryi,j)ui,jkui,jk1(4)
得到非均匀介质的波动方程差分形式。

3. 初值条件及边界条件

3.1 初值条件

选用震源:
u ( x , y , 0 ) = a e − b [ ( x − x 0 ) 2 + ( y − y 0 ) 2 ] u(x, y, 0) = ae^{-b\left[(x-x_0)^2 + (y-y_0)^2\right]} u(x,y,0)=aeb[(xx0)2+(yy0)2]
x 0 , y 0 x_0, y_0 x0,y0 为震源位置, a , b a, b a,b为系数。

特别的, 从 u i , j 0 u_{i, j}^0 ui,j0 递推到 u i , j 1 u_{i, j}^1 ui,j1 时, 需要用到 u i , j − 1 u_{i, j}^{-1} ui,j1 时的值, 但我们并没有当前时步的波场情况。这里我们使用外插法来获得 − 1 -1 1 时步的波长, 即
u i , j − 1 = u i , j 0 − 2 Δ t v 0 ( x i , y j ) u_{i, j}^{-1} = u_{i,j}^0 - 2\Delta t v_0(x_i, y_j) ui,j1=ui,j0tv0(xi,yj)
v 0 v_0 v0 为初始状态下每一点的波动速度。

另外, 三点差分格式必须满足收敛条件:
4 v 2 Δ t 2 Δ x 2 + Δ y 2 ≤ 1 \frac{4 v^2 \Delta t^2}{\Delta x^2 + \Delta y^2} \le 1 Δx2+Δy24v2Δt21

3.2 边界条件

波在传播到实验区域边界后, 如果将边界的波场值设定恒定 0 0 0, 那么会产生反射波,影响波长结果, 为此我们需要一种吸收边界来抵消边界产生的反射波。

其方程满足如下形式:
∂ u ∂ n − 1 v ∂ u ∂ t = 0 \frac{\partial u}{\partial n} - \frac{1}{v}\frac{\partial u}{\partial t} = 0 nuv1tu=0
n n n 为边界的法向量。

以顶部的边界条件为例, 其差分形式满足:
u 1 , j k − u 0 , j k Δ y = 1 v u 0 , j k + 1 − u 0 , j k Δ t \frac{u_{1, j}^k - u_{0, j}^k}{\Delta y} = \frac{1}{v}\frac{u^{k +1}_{0, j} - u^k_{0, j}}{\Delta t} Δyu1,jku0,jk=v1Δtu0,jk+1u0,jk

u 0 , j k + 1 = ( 1 − A y ) u 0 , j k + A y ⋅ u 1 , j k (5) u_{0, j}^{k+1} = (1 - A_y) u_{0, j}^k + A_y \cdot u_{1, j}^k \tag{5} u0,jk+1=(1Ay)u0,jk+Ayu1,jk(5)

A y = Δ t Δ y v A_y = \frac{\Delta t}{\Delta y}v Ay=ΔyΔtv

我们用与 ( 4 ) (4) (4)同样的方法可以将 A y A_y Ay 代换为 A y i , j A_y^{i, j} Ayi,j 得到非均匀介质下的吸收条件。

其余三个边界的吸收条件与 ( 5 ) (5) (5) 式类似,不再赘述。

4. 代码实现

4.1 约定

mask蒙层

在传递参数的过程中, 我们传递一个数组 v v v 表示一组波速情况, 程序会构建一个蒙层 m a s k mask mask 来记录每一个点的波速索引。其中 m a s k mask mask 与波场具有同样的大小, m a s k [ i , j ] mask[i, j] mask[i,j] 表示 ( x i , y j ) (x_i, y_j) (xi,yj) 处的波的波速所在 v v v数组的下标。

比如:

我们假设
v = [ 100 , 200 , 300 ] v = [100, 200, 300]\\ v=[100,200,300]
如果 m a s k [ 1 , 2 ] = 1 mask[1, 2] = 1 mask[1,2]=1, 表示 ( x 1 , y 2 ) (x_1, y_2) (x1,y2) 处的波速为 v [ 1 ] v[1] v[1], 即 200 200 200,

如果 m a s k [ 5 , 3 ] = 0 mask[5, 3] = 0 mask[5,3]=0, 表示 ( x 5 , y 3 ) (x_5, y_3) (x5,y3) 处的波速为 v [ 0 ] v[0] v[0], 即 100 100 100

sep函数

在传递参数的过程中需要传递一个 s e p sep sep参数来指定生成蒙层的方法。 s e p sep sep 是一个函数, 该函数传入两个参数 i , j i, j i,j 返回 m a s k [ i , j ] mask[i, j] mask[i,j] 的值。

比如:

sep = lambda x, y: 0 if x > 2 else 1

表示当 x i > 2 时 x_i > 2 时 xi>2 , m a s k [ i , j ] mask[i, j] mask[i,j] 的值为 0 0 0, 否则为 1 1 1

4.2 代码

import osimport numpy as np
import matplotlib.pyplot as plt
import seaborn as snsclass WaveEquation2D:# initify horizion wave equation 2d simulator.def __init__(self, mt, mx, my, nt, nx, ny, px, py, v, sep):self.mt = mtself.mx = mxself.my = myself.nt = ntself.nx = nxself.ny = nyself.v = np.array(v)self.sep = sepself.dx = self.mx / self.nxself.dy = self.my / self.nyself.dt = self.mt / self.ntself.x = np.arange(0, self.mx + self.dx, self.dx)self.y = np.arange(0, self.my + self.dy, self.dy)self.t = np.arange(0, self.mt + self.dt, self.dt)self.u0 = lambda r, c: 0.2 * \np.exp(-((r - px)**2 + (c - py)**2) / 0.01)self.v0 = lambda r, c: 0r = 4 * max(self.v) ** 2 * self.dt ** 2 / (self.dx ** 2 + self.dy ** 2)assert r < 1, f"r should be less then 1, r is {r}."self.Ax = self.v * self.dt / self.dxself.Ay = self.v * self.dt / self.dyself.rx = self.Ax ** 2self.ry = self.Ay ** 2self.rxy = 1 - self.rx - self.ryself.u = np.zeros((self.nt + 1, self.nx + 1, self.ny + 1))for i in range(1, self.nx):for j in range(1, self.ny):self.u[0, i, j] = self.u0(self.x[i], self.y[j])self.mask = np.array([[sep(self.x[i], self.y[j]) for j in range(0, ny + 1)]for i in range(0, nx + 1)])def info(self):print("wave arguments:")print(f"dx: {self.dx}, dy: {self.dy}, dt:{self.dt}")print(f"ax: {self.Ax}, ay: {self.Ay}")print(f"rx: {self.rx}, ry: {self.ry}, rxy: {self.rxy}")def apply_boundary(self, t):# absort 1 dimfor i in range(self.nx + 1):self.u[t, i, 0] =\self.Ay[self.mask[i, 1]] * self.u[t-1, i, 1] + \(1 - self.Ay[self.mask[i, 0]]) * self.u[t-1, i, 0]self.u[t, i, self.ny] =\self.Ay[self.mask[i, self.ny - 1]] * \self.u[t-1, i, self.ny-1] + \(1-self.Ay[self.mask[i, self.ny]]) * self.u[t-1, i, self.ny]for j in range(self.ny + 1):self.u[t, 0, j] =\self.Ax[self.mask[1, j]] * self.u[t-1, 1, j] + \(1 - self.Ax[self.mask[0, j]]) * self.u[t-1, 0, j]self.u[t, self.nx, j] =\self.Ax[self.mask[self.nx-1, j]] * \self.u[t-1, self.nx-1, j] + \(1-self.Ax[self.mask[self.nx, j]]) * self.u[t-1, self.nx, j]def solve(self):# solve t == 1for i in range(1, self.nx):for j in range(1, self.ny):self.u[1, i, j] =\0.5 * self.rx[self.mask[i-1, j]] * self.u[0, i-1, j] + \0.5 * self.rx[self.mask[i+1, j]] * self.u[0, i+1, j] + \0.5 * self.ry[self.mask[i, j-1]] * self.u[0, i, j-1] + \0.5 * self.ry[self.mask[i, j+1]] * self.u[0, i, j+1] + \self.rxy[self.mask[i, j]] * self.u[0, i, j] + \self.dt * self.v0(i, j)for t in range(2, self.nt):print('\rsolving step {} / {}'.format(t + 1, self.nt), end="")# fills boundary conditiosself.apply_boundary(t)for i in range(1, self.nx):for j in range(1, self.ny):self.u[t, i, j] = \self.rx[self.mask[i-1, j]] * self.u[t-1, i-1, j] + \self.rx[self.mask[i+1, j]] * self.u[t-1, i+1, j] + \self.ry[self.mask[i, j-1]] * self.u[t-1, i, j-1] + \self.ry[self.mask[i, j+1]] * self.u[t-1, i, j+1] + \2 * self.rxy[self.mask[i, j]] * self.u[t-1, i, j] -\self.u[t-2, i, j]def plot(self, img_path="./image", show=True, save=True):if not os.path.exists(img_path):os.mkdir(img_path)for t in range(0, self.nt):sns.heatmap(self.u[t, :, :],vmin=-np.max(self.u)/10,vmax=np.max(self.u)/10,cmap='seismic',xticklabels=False,yticklabels=False)if save:plt.savefig(os.path.join(img_path, str(t)))if show:plt.pause(0.01)plt.cla()plt.clf()if __name__ == '__main__':we = WaveEquation2D(mt=2,mx=4,my=4,nt=300,nx=80,ny=80,px=0,py=2,v=[3.0, 5.0],sep=lambda r, c: 1 if r > 2 else 0)we.solve()we.plot(show=True, save=True)

附录

A. 收敛性条件

为了方便后面描述,定义波场精确值
u i , j , k = u ( x i , y j , t k ) u_{i, j, k} = u(x_i, y_j, t_k) ui,j,k=u(xi,yj,tk)
与波场的差分值 u i , j k u_{i, j}^{k} ui,jk区分。

和二阶差分算子

δ t 2 U = U i , j , k + 1 − 2 U i , j , k + U i , j , k − 1 δ x 2 U = U i + 1 , j , k − 2 U i , j , k + U i − 1 , j , k δ y 2 U = U i , j + 1 , k − 2 U i , j , k + U i , j − 1 , k \delta_{t}^{2} U = U_{i, j, k + 1} - 2U_{i, j, k} + U_{i, j, k-1}\\ \delta_{x}^{2} U = U_{i+1, j, k} - 2U_{i, j, k} + U_{i-1, j, k}\\ \delta_{y}^{2} U = U_{i, j+1, k} - 2U_{i, j, k} + U_{i, j-1, k} δt2U=Ui,j,k+12Ui,j,k+Ui,j,k1δx2U=Ui+1,j,k2Ui,j,k+Ui1,j,kδy2U=Ui,j+1,k2Ui,j,k+Ui,j1,k

算子的右下部分表示所要二阶差分的坐标。

精确值的差分格式可以写为

δ t 2 u Δ t 2 = v 2 [ δ x 2 u Δ x 2 + δ y 2 u Δ y 2 ] \dfrac{\delta^2_t u}{\Delta t^2}=v^2\left[ \dfrac{\delta^2_x u}{\Delta x^2} + \dfrac{\delta^2_y u}{\Delta y^2} \right] Δt2δt2u=v2[Δx2δx2u+Δy2δy2u]

在差分格式下,我们取截断误差
T i , j , k = δ t 2 u Δ t 2 − v 2 [ δ x 2 u Δ x 2 + δ y 2 u Δ y 2 ] T_{i, j, k} = \dfrac{\delta^2_t u}{\Delta t^2}-v^2\left[ \dfrac{\delta^2_x u}{\Delta x^2} + \dfrac{\delta^2_y u}{\Delta y^2} \right] Ti,j,k=Δt2δt2uv2[Δx2δx2u+Δy2δy2u]

此时截断误差是 u u u 为精确值时,由于差分格式所引入的误差。

上式经过化简可得
u i , j , k + 1 = r x ( u i − 1 , j , k + u i + 1 , j , k ) + r y ( u i , j − 1 , k , u i , j + 1 , k ) + 2 ( 1 − r x − r y ) u i , j , k − u i , j , k − 1 − ( Δ t ) 2 T i , j , k (6) \begin{aligned} u_{i, j, k+1} &= r_x{(u_{i -1, j, k}+ u_{i+1, j, k})} \\ & + r_y(u_{i, j-1, k}, u_{i, j +1, k}) \\ &+ 2(1 - r_x - r_y) u_{i, j, k}\\ &- u_{i, j, k-1}\\ & - (\Delta t)^2 T_{i, j, k} \end{aligned} \tag{6} ui,j,k+1=rx(ui1,j,k+ui+1,j,k)+ry(ui,j1,k,ui,j+1,k)+2(1rxry)ui,j,kui,j,k1(Δt)2Ti,j,k(6)

定义真实误差
e i , j , k = u i , j k − u i , j , k e_{i, j, k} = u_{i, j}^k - u_{i, j, k} ei,j,k=ui,jkui,j,k
也就是在网格点上函数真实值与差分估计值的差。

( 6 ) (6) (6) 式 减去 ( 3 ) (3) (3) 式,可得

e i , j , k + 1 = r x ( e i − 1 , j , k + e i + 1 , j , k ) + r y ( e i , j − 1 , k , e i , j + 1 , k ) + 2 ( 1 − r x − r y ) e i , j , k − e i , j , k − 1 − ( Δ t ) 2 T i , j , k (7) \begin{aligned} e_{i, j, k+1} &= r_x{(e_{i -1, j, k}+ e_{i+1, j, k})} \\ & + r_y(e_{i, j-1, k}, e_{i, j +1, k}) \\ &+ 2(1 - r_x - r_y) e_{i, j, k}\\ &- e_{i, j, k-1}\\ & - (\Delta t)^2 T_{i, j, k} \end{aligned} \tag{7} ei,j,k+1=rx(ei1,j,k+ei+1,j,k)+ry(ei,j1,k,ei,j+1,k)+2(1rxry)ei,j,kei,j,k1(Δt)2Ti,j,k(7)

如果满足 1 − r x − r y ≥ 0 1-r_x-r_y \ge 0 1rxry0, 则上式所有与 e e e有关的系数均为整数,打绝对值以后可以提到绝对值号的外面。
由三角不等式有
∣ e i , j , k + 1 ∣ ≤ r x ( ∣ e i − 1 , j , k ∣ + ∣ e i + 1 , j , k ∣ ) + r y ( ∣ e i , j − 1 , k ∣ + ∣ e i , j + 1 , k ∣ ) + 2 ( 1 − r x − r y ) ∣ e i , j , k ∣ + ∣ e i , j , k − 1 ∣ + ∣ ( Δ t ) 2 T i , j , k ∣ (8) \begin{aligned} |e_{i, j, k+1}| & \le r_x{(|e_{i -1, j, k}|+ |e_{i+1, j, k}|)} \\ & + r_y(|e_{i, j-1, k}|+ |e_{i, j +1, k}|) \\ &+ 2(1 - r_x - r_y) |e_{i, j, k}|\\ & + |e_{i, j, k-1}|\\ & + |(\Delta t)^2 T_{i, j, k}| \end{aligned} \tag{8} ei,j,k+1rx(ei1,j,k+ei+1,j,k)+ry(ei,j1,k+ei,j+1,k)+2(1rxry)ei,j,k+ei,j,k1+(Δt)2Ti,j,k(8)

将 截断误差 T T T 的定义式进行泰勒展开,可知 T = O ( Δ t 2 + Δ x 2 + Δ y 2 ) T = O(\Delta t^2 + \Delta x^2 + \Delta y^2) T=O(Δt2+Δx2+Δy2), 那么必存在 M > 0 M > 0 M>0, 使得 任意的 T T T 都有
∣ T ∣ ≤ M ( Δ t 2 + Δ x 2 + Δ y 2 ) |T| \le M(\Delta t^2 + \Delta x^2 + \Delta y^2) TM(Δt2+Δx2+Δy2).
再令
E k = sup ⁡ i , j ∣ e i , j , k ∣ E_k = \sup_{i, j}|e_{i, j, k}| Ek=i,jsupei,j,k
代入 (8) 得
∣ e i , j , k + 1 ∣ ≤ r x ( E k + E k ) + r y ( E k + E k ) + 2 ( 1 − r x − r y ) E m + E k − 1 + ( Δ t ) 2 M ( Δ t 2 + Δ x 2 + Δ y 2 ) \begin{aligned} |e_{i, j, k+1}| & \le r_x{(E_k+ E_k)} \\ & + r_y(E_k+ E_k) \\ &+ 2(1 - r_x - r_y) E_m\\ & + E_{k-1}\\ & + (\Delta t)^2 M(\Delta t^2 + \Delta x^2 + \Delta y^2) \end{aligned} ei,j,k+1rx(Ek+Ek)+ry(Ek+Ek)+2(1rxry)Em+Ek1+(Δt)2M(Δt2+Δx2+Δy2)
于是有
E k + 1 ≤ E k + E k − 1 + ( Δ t ) 2 M ( Δ t 2 + Δ x 2 + Δ y 2 ) E_{k+1} \le E_k + E_{k-1} + (\Delta t)^2 M(\Delta t^2 + \Delta x^2 + \Delta y^2) Ek+1Ek+Ek1+(Δt)2M(Δt2+Δx2+Δy2)
递推可得
E k ≤ 2 E 0 + k ( Δ t ) 2 M ( Δ t 2 + Δ x 2 + Δ y 2 ) E_{k} \le 2 E_0 + k (\Delta t)^2 M(\Delta t^2 + \Delta x^2 + \Delta y^2) Ek2E0+k(Δt)2M(Δt2+Δx2+Δy2)
由于在初值给定,必然有 E 0 = 0 E_0 = 0 E0=0,
于是此时有
E k ≤ k ( Δ t ) 2 M ( Δ t 2 + Δ x 2 + Δ y 2 ) E_{k} \le k (\Delta t)^2 M(\Delta t^2 + \Delta x^2 + \Delta y^2) Ekk(Δt)2M(Δt2+Δx2+Δy2)

Δ t , Δ x , Δ y → 0 \Delta t, \Delta x, \Delta y \to 0 Δt,Δx,Δy0时, E k → 0 E_k \to 0 Ek0, 可以说明收敛性。

综上,在条件 1 − r x − r y ≥ 0 1 - r_x - r_y \ge 0 1rxry0 时, 即
r x + r y ≤ 1 Δ t 2 v 2 Δ x 2 + Δ t 2 v 2 Δ y 2 ≤ 1 Δ t 2 v 2 ( Δ x 2 + Δ y 2 ) Δ x 2 Δ y 2 ≤ 1 r_x + r_y \le 1 \\ \dfrac{\Delta t^2 v^2}{\Delta x^2} + \dfrac{\Delta t^2 v^2}{\Delta y^2} \le 1\\ \dfrac{\Delta t^2 v^2 (\Delta x^2 + \Delta y^2)}{\Delta x^2 \Delta y^2} \le 1 rx+ry1Δx2Δt2v2+Δy2Δt2v21Δx2Δy2Δt2v2(Δx2+Δy2)1

由于
a ⋅ b ≤ ( a + b ) 2 4 a^\cdot b \le \frac{(a+b)^2}{4} ab4(a+b)2


4 v 2 Δ t 2 Δ x 2 + Δ y 2 ≤ 1 \frac{4 v^2 \Delta t^2}{\Delta x^2 + \Delta y^2} \le 1 Δx2+Δy24v2Δt21
该条件是一个十分严苛的条件。

这篇关于二维波动方程数值模拟(非均匀介质)的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/703920

相关文章

poj2576(二维背包)

题意:n个人分成两组,两组人数只差小于1 , 并且体重只差最小 对于人数要求恰好装满,对于体重要求尽量多,一开始没做出来,看了下解题,按照自己的感觉写,然后a了 状态转移方程:dp[i][j] = max(dp[i][j],dp[i-1][j-c[k]]+c[k]);其中i表示人数,j表示背包容量,k表示输入的体重的 代码如下: #include<iostream>#include<

hdu2159(二维背包)

这是我的第一道二维背包题,没想到自己一下子就A了,但是代码写的比较乱,下面的代码是我有重新修改的 状态转移:dp[i][j] = max(dp[i][j], dp[i-1][j-c[z]]+v[z]); 其中dp[i][j]表示,打了i个怪物,消耗j的耐力值,所得到的最大经验值 代码如下: #include<iostream>#include<algorithm>#include<

【C++】_list常用方法解析及模拟实现

相信自己的力量,只要对自己始终保持信心,尽自己最大努力去完成任何事,就算事情最终结果是失败了,努力了也不留遗憾。💓💓💓 目录   ✨说在前面 🍋知识点一:什么是list? •🌰1.list的定义 •🌰2.list的基本特性 •🌰3.常用接口介绍 🍋知识点二:list常用接口 •🌰1.默认成员函数 🔥构造函数(⭐) 🔥析构函数 •🌰2.list对象

usaco 1.2 Transformations(模拟)

我的做法就是一个一个情况枚举出来 注意计算公式: ( 变换后的矩阵记为C) 顺时针旋转90°:C[i] [j]=A[n-j-1] [i] (旋转180°和270° 可以多转几个九十度来推) 对称:C[i] [n-j-1]=A[i] [j] 代码有点长 。。。 /*ID: who jayLANG: C++TASK: transform*/#include<

【生成模型系列(初级)】嵌入(Embedding)方程——自然语言处理的数学灵魂【通俗理解】

【通俗理解】嵌入(Embedding)方程——自然语言处理的数学灵魂 关键词提炼 #嵌入方程 #自然语言处理 #词向量 #机器学习 #神经网络 #向量空间模型 #Siri #Google翻译 #AlexNet 第一节:嵌入方程的类比与核心概念【尽可能通俗】 嵌入方程可以被看作是自然语言处理中的“翻译机”,它将文本中的单词或短语转换成计算机能够理解的数学形式,即向量。 正如翻译机将一种语言

HDU 2159 二维完全背包

FATE 最近xhd正在玩一款叫做FATE的游戏,为了得到极品装备,xhd在不停的杀怪做任务。久而久之xhd开始对杀怪产生的厌恶感,但又不得不通过杀怪来升完这最后一级。现在的问题是,xhd升掉最后一级还需n的经验值,xhd还留有m的忍耐度,每杀一个怪xhd会得到相应的经验,并减掉相应的忍耐度。当忍耐度降到0或者0以下时,xhd就不会玩这游戏。xhd还说了他最多只杀s只怪。请问他能

hdu4431麻将模拟

给13张牌。问增加哪些牌可以胡牌。 胡牌有以下几种情况: 1、一个对子 + 4组 3个相同的牌或者顺子。 2、7个不同的对子。 3、13幺 贪心的思想: 对于某张牌>=3个,先减去3个相同,再组合顺子。 import java.io.BufferedInputStream;import java.io.BufferedReader;import java.io.IOExcepti

【每日一题】LeetCode 2181.合并零之间的节点(链表、模拟)

【每日一题】LeetCode 2181.合并零之间的节点(链表、模拟) 题目描述 给定一个链表,链表中的每个节点代表一个整数。链表中的整数由 0 分隔开,表示不同的区间。链表的开始和结束节点的值都为 0。任务是将每两个相邻的 0 之间的所有节点合并成一个节点,新节点的值为原区间内所有节点值的和。合并后,需要移除所有的 0,并返回修改后的链表头节点。 思路分析 初始化:创建一个虚拟头节点

每日一题|牛客竞赛|四舍五入|字符串+贪心+模拟

每日一题|四舍五入 四舍五入 心有猛虎,细嗅蔷薇。你好朋友,这里是锅巴的C\C++学习笔记,常言道,不积跬步无以至千里,希望有朝一日我们积累的滴水可以击穿顽石。 四舍五入 题目: 牛牛发明了一种新的四舍五入应用于整数,对个位四舍五入,规则如下 12345->12350 12399->12400 输入描述: 输入一个整数n(0<=n<=109 ) 输出描述: 输出一个整数

【算法专场】模拟(下)

目录 前言 38. 外观数列 算法分析 算法思路 算法代码 1419. 数青蛙 算法分析 算法思路 算法代码  2671. 频率跟踪器 算法分析 算法思路 算法代码 前言 在前面我们已经讲解了什么是模拟算法,这篇主要是讲解在leetcode上遇到的一些模拟题目~ 38. 外观数列 算法分析 这道题其实就是要将连续且相同的字符替换成字符重复的次数+