波动方程数值求解（一）

已于 2022-07-11 16:36:59 修改
阅读量5.3k 收藏 52
点赞数 9
分类专栏： 波动方程 文章标签： python numpy
于 2022-07-11 14:28:19 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/qq_45935215/article/details/125719067
版权

波动方程数值解是波动方程正演、逆时偏移和全波形反演的核心技术之一。本文采用二阶有限差分对波动方程进行了离散，进而实现了对波动方程的数值求解，模拟出其在介质中的传播过程。
1、二维声波波动方程离散
在这里插入图片描述
利用泰勒公式进行展开得到：
在这里插入图片描述
两式相减得：
在这里插入图片描述
则有：
在这里插入图片描述
近似得二阶差分算子：
在这里插入图片描述
利用二阶中心差分算子对二阶导数进行离散：
在这里插入图片描述
将上式代入声波方程得到二阶中心差分格式：
在这里插入图片描述
其中：在这里插入图片描述
收敛满足：
在这里插入图片描述

其空间和时间差分格式示意图如下图所示：
在这里插入图片描述
在这里插入图片描述
2、均匀介质二维声波方程求解
初始参数：

# 区域大小
Nx = 301
Nz = 301
# 时间采样数
Nt = 1000
# 速度模型
v = np.ones((Nx+1,Nz+1)) * 3000
u = np.zeros((Nt+1,Nx+1,Nz+1))
h = 5
# 子波主频
fm = 25 
# 空间衰减因子
alpha = 0.5
# 时间步
dt = 1 / Nt
# 空间间隔
dx = h
dy = h
# 迭代公式中的r
A = v **2 * dt ** 2 / h ** 2 
C = v * dt / h

收敛条件判断：

r = np.max(v)*dt/h
assert r < 0.707,f'r should < 0.707, but is {r}'

雷克子波及二维空间衰减函数：

t = np.arange(0, Nt+1)
t0 = 0  # 延迟时间，相当于在t=t0时激发 ，震幅在t0时最大，相位也在此
s_t = (1 - 2 * (np.pi * fm * dt * (t - t0)) ** 2) * np.exp( - (np.pi * fm * dt * (t - t0)) ** 2)

x = np.arange(0,Nx+1)
z = np.arange(0,Nz+1)
x,z = np.meshgrid(x,z)

h_x_z =  np.zeros((Nx+1,Nz+1))  #np.exp(-0.25 * ((x - Nx//2)**2 + (z - Nz//2)**2))  二维空间衰减函数
h_x_z[Nx//2,Nz//2] = 1   # 在Nx//2,Nz//2处激发
h_x_z =  np.exp(-alpha ** 2 * ((x - Nx//2)**2 + (z - Nz//2)**2))  # 二维空间衰减函数

x0 = Nx // 2
z0 = Nz // 2
u0 = lambda r, s: 0.25*np.exp(-((r-x0)**2+(s-z0)**2))
JJ = np.arange(1,Nz)
II = np.arange(1,Nx)
II,JJ = np.meshgrid(II,JJ)

吸收边界条件：
Clayton-Engquist 单成波吸收边界条件:最早是由 Clayton 等人发现并推广的， 其微分表达式为：

在这里插入图片描述
其中：n 为边界的外法线方向； s 边界的切线方向。
对上式进行离散得到上、下、左、右边界差分格式如下：
在这里插入图片描述
在这里插入图片描述
其中： 在这里插入图片描述N 、 M 为边界的网格数。
则有：

# 边界条件
ii = np.arange(Nx+1)
jj = np.arange(Nz+1)
# Clayton-Engquist-majda 二阶吸收边界条件
u[t+1,  0, jj] = (2 - 2 * C[ 0, jj] - C[ 0, jj] ** 2) * u[t,  0, jj] \
                        + 2 * C[ 1, jj] * (1 + C[ 1, jj]) * u[t,  1, jj] \
                        - C[ 2, jj] ** 2 * u[t,  2, jj] \
                        + (2 * C[ 0, jj] - 1) * u[t - 1,  0, jj] \
                        - 2 * C[ 1, jj] * u[t - 1,  1, jj]


# 下部
u[t+1, -1, jj] = (2 - 2 * C[ -1, jj] - C[ -1, jj] ** 2) * u[t,  -1, jj] \
                        + 2 * C[ -2, jj] * (1 + C[ -2, jj]) * u[t,  -2, jj] \
                        - C[ -3, jj] ** 2 * u[t,  -3, jj] \
                        + (2 * C[ -1, jj] - 1) * u[t - 1,  -1, jj] \
                        - 2 * C[ -2, jj] * u[t - 1,  -2, jj]

# 左部
u[t+1, ii,  0] = (2 - 2 * C[ii,  0] - C[ii,  0] ** 2) * u[t,  ii, 0] \
                        + 2 * C[ii,  1] * (1 + C[ii,  1]) * u[t,  ii, 1] \
                        - C[ii,  2] ** 2 * u[t,  ii, 2] \
                        + (2 * C[ii,  0] - 1) * u[t - 1,  ii, 0] \
                        - 2 * C[ii,  1] * u[t - 1,  ii, 1]
# 右部
u[t+1, ii, -1] = (2 - 2 * C[ii,  -1] - C[ii,  -1] ** 2) * u[t,  ii, -1] \
                        + 2 * C[ii,  -2] * (1 + C[ii,  -2]) * u[t,  ii, -2] \
                        - C[ii,  -3] ** 2 * u[t,  ii, -3] \
                        + (2 * C[ii,  -1] - 1) * u[t - 1,  ii, -1] \
                        - 2 * C[ii,  -2] * u[t - 1,  ii, -2]

Reynolds 边界条件：对于二维声波波动方程，应用二维声波方的微分算子可以 得到：

在这里插入图片描述
对上式进行离散可得上下左右边界计算公式：
在这里插入图片描述
在这里插入图片描述
则有：

#Reynolds 边界条件
u[t+1,ii, 0] = u[t,ii, 0] + u[t,ii, 1] - u[t-1,ii, 1] + C[ii, 1]*u[t,ii, 1] - C[ii, 0]*u[t,ii, 0] -C[ii, 2]*u[t-1,ii, 2] +C[ii, 1]*u[t-1,ii, 1]

u[t+1,ii,-1] = u[t,ii,-1] + u[t,ii,-2] - u[t-1,ii,-2] + C[ii,-2]*u[t,ii,-2] - C[ii,-1]*u[t,ii,-1] -C[ii,-3]*u[t-1,ii,-3] +C[ii,-2]*u[t-1,ii,-2]

u[t+1, 0,jj] = u[t, 0,jj] + u[t, 1,jj] - u[t-1, 1,jj] + C[ 1,jj]*u[t, 1,jj] - C[ 0,jj]*u[t, 0,jj] -C[ 2,jj]*u[t-1, 2,jj] +C[ 1,jj]*u[t-1, 1,jj]

u[t+1,-1,jj] = u[t,-1,jj] + u[t,-2,jj] - u[t-1,-2,jj] + C[-2,jj]*u[t,-2,jj] - C[-1,jj]*u[t,-1,jj] -C[-3,jj]*u[t-1,-3,jj] +C[-1,jj]*u[t-1,-2,jj]

迭代公式：

# 迭代公式
u[t+1,II,JJ] = s_t[t]*h_x_z[II,JJ]+A[II,JJ]*(u[t,II,JJ+1]+u[t,II,JJ-1]+u[t,II+1,JJ]+u[t,II-1,JJ])+(2-4*A[II,JJ])*u[t,II,JJ]-u[t-1,II,JJ]

所有代码如下：

import numpy as np
import imageio.v2 as imageio
import os
import pandas as pd
from matplotlib import pyplot as plt
# 解决中文问题
plt.rcParams['font.sans-serif'] = ['SimHei']
# 解决负号显示问题
plt.rcParams['axes.unicode_minus'] = False


Nx = 301
Nz = 301
Nt = 1000
v = np.ones((Nx+1,Nz+1)) * 3000
h = 5
fm = 25
alpha = 0.5
dt = 1 / Nt
dx = h
dy = h
A = v **2 * dt ** 2 / h ** 2 
C = v * dt / h

r = np.max(v)*dt/h
assert r < 0.707,f'r should < 0.707, but is {r}'
u = np.zeros((Nt+1,Nx+1,Nz+1))


t = np.arange(0, Nt+1)
t0 = 0  # 延迟时间，相当于在t=t0时激发 ，震幅在t0时最大，相位也在此
s_t = (1 - 2 * (np.pi * fm * dt * (t - t0)) ** 2) * np.exp( - (np.pi * fm * dt * (t - t0)) ** 2)

plt.plot(s_t)
plt.show()

x = np.arange(0,Nx+1)
z = np.arange(0,Nz+1)
x,z = np.meshgrid(x,z)
# print(((z - Nz//2)**2).shape)


h_x_z =  np.zeros((Nx+1,Nz+1))  #np.exp(-0.25 * ((x - Nx//2)**2 + (z - Nz//2)**2))  二维空间衰减函数
h_x_z[Nx//2,Nz//2] = 1   # 在Nx//2,Nz//2处激发
h_x_z =  np.exp(-alpha ** 2 * ((x - Nx//2)**2 + (z - Nz//2)**2))  # 二维空间衰减函数


JJ = np.arange(1,Nz)
II = np.arange(1,Nx)
II,JJ = np.meshgrid(II,JJ)

mode = 'c_e'
img_path = './2_order'
if not os.path.exists(img_path):
    os.makedirs(img_path)

for t in range(2,Nt):
    print('\rstep {} / {}'.format(t ,Nt), end="")
    
    # 边界条件
    ii = np.arange(Nx+1)
    jj = np.arange(Nz+1)
    if mode == 'c_e':
        # Clayton-Engquist-majda 二阶吸收边界条件
        u[t+1,  0, jj] = (2 - 2 * C[ 0, jj] - C[ 0, jj] ** 2) * u[t,  0, jj] \
                                + 2 * C[ 1, jj] * (1 + C[ 1, jj]) * u[t,  1, jj] \
                                - C[ 2, jj] ** 2 * u[t,  2, jj] \
                                + (2 * C[ 0, jj] - 1) * u[t - 1,  0, jj] \
                                - 2 * C[ 1, jj] * u[t - 1,  1, jj]


        # 下部
        u[t+1, -1, jj] = (2 - 2 * C[ -1, jj] - C[ -1, jj] ** 2) * u[t,  -1, jj] \
                                + 2 * C[ -2, jj] * (1 + C[ -2, jj]) * u[t,  -2, jj] \
                                - C[ -3, jj] ** 2 * u[t,  -3, jj] \
                                + (2 * C[ -1, jj] - 1) * u[t - 1,  -1, jj] \
                                - 2 * C[ -2, jj] * u[t - 1,  -2, jj]

        # 左部
        u[t+1, ii,  0] = (2 - 2 * C[ii,  0] - C[ii,  0] ** 2) * u[t,  ii, 0] \
                                + 2 * C[ii,  1] * (1 + C[ii,  1]) * u[t,  ii, 1] \
                                - C[ii,  2] ** 2 * u[t,  ii, 2] \
                                + (2 * C[ii,  0] - 1) * u[t - 1,  ii, 0] \
                                - 2 * C[ii,  1] * u[t - 1,  ii, 1]
        # 右部
        u[t+1, ii, -1] = (2 - 2 * C[ii,  -1] - C[ii,  -1] ** 2) * u[t,  ii, -1] \
                                + 2 * C[ii,  -2] * (1 + C[ii,  -2]) * u[t,  ii, -2] \
                                - C[ii,  -3] ** 2 * u[t,  ii, -3] \
                                + (2 * C[ii,  -1] - 1) * u[t - 1,  ii, -1] \
                                - 2 * C[ii,  -2] * u[t - 1,  ii, -2]
    if mode == 're':
        #Reynolds 边界条件
        u[t+1,ii, 0] = u[t,ii, 0] + u[t,ii, 1] - u[t-1,ii, 1] + C[ii, 1]*u[t,ii, 1] - C[ii, 0]*u[t,ii, 0] -C[ii, 2]*u[t-1,ii, 2] +C[ii, 1]*u[t-1,ii, 1]

        u[t+1,ii,-1] = u[t,ii,-1] + u[t,ii,-2] - u[t-1,ii,-2] + C[ii,-2]*u[t,ii,-2] - C[ii,-1]*u[t,ii,-1] -C[ii,-3]*u[t-1,ii,-3] +C[ii,-2]*u[t-1,ii,-2]

        u[t+1, 0,jj] = u[t, 0,jj] + u[t, 1,jj] - u[t-1, 1,jj] + C[ 1,jj]*u[t, 1,jj] - C[ 0,jj]*u[t, 0,jj] -C[ 2,jj]*u[t-1, 2,jj] +C[ 1,jj]*u[t-1, 1,jj]

        u[t+1,-1,jj] = u[t,-1,jj] + u[t,-2,jj] - u[t-1,-2,jj] + C[-2,jj]*u[t,-2,jj] - C[-1,jj]*u[t,-1,jj] -C[-3,jj]*u[t-1,-3,jj] +C[-1,jj]*u[t-1,-2,jj]


    # 迭代公式
    u[t+1,II,JJ] = s_t[t]*h_x_z[II,JJ]+A[II,JJ]*(u[t,II,JJ+1]+u[t,II,JJ-1]+u[t,II+1,JJ]+u[t,II-1,JJ])+(2-4*A[II,JJ])*u[t,II,JJ]-u[t-1,II,JJ]


    
    plt.imshow(u[t], cmap='gray_r') # 'seismic' gray_r
    # plt.axis('off')  # 隐藏坐标轴
    plt.colorbar()
    if t % 10 == 0:
        plt.savefig(os.path.join(img_path, str(t) + '.png'),
                    bbox_inches="tight", # 去除坐标轴占用空间
                    pad_inches=0.0) # 去除所有白边
        plt.pause(0.05)
        
    plt.cla()
    plt.clf()  # 清除所有轴，但是窗口打开，这样它可以被重复使用
plt.close()

# 保存gif
filenames=[]
for files in os.listdir(img_path ):
    if files.endswith('jpg') or files.endswith('jpeg') or files.endswith('png'):
        file=os.path.join(img_path ,files)
        filenames.append(file)

images = []
for filename in filenames:
    images.append(imageio.imread(filename))

imageio.mimsave(os.path.join(img_path,  'wave.gif'), images,duration=0.0001)

最终得到有边界为0、吸收条件结果如下：
在这里插入图片描述

在这里插入图片描述
生成gif的代码如上，但有闪烁的问题如上图所示，本人无法解决，若有知道如何解决的，在线求教。

波动方程正演模型
04-27
波动方程正演实验适合地震层析和正演方向的同学，模型建立，源规范，时间步长，Siesmograms，波场快照
二维波动方程数值模拟(非均匀介质)
Sumbrella_的博客
12-09 7091
文章用有限差分法模拟并可视化了非均匀介质波场，可以观测到反射波、折射波等。
二维波动方程求解
06-01
利用差值算法逼近解析解，图形显示，对比数值解，可视化。
波动方程数值求解（二）
qq_45935215的博客
07-11 1438
波动方程数值模拟 波动方程数值解是波动方程正演、逆时偏移和全波形反演的核心技术之一。本文采用高阶有限差分对波动方程进行了离散，进而实现了对波动方程的数值求解，模拟出其在介质中的传播过程。...
Python——一维波动方程求解
Thattear的博客
07-02 3422
1、一维波动方程简介 以线的横振动方程作为一维波动方程的例子，其运动学方程满足： 在上式中为线密度；T为张力；代表强迫力的大小。当没有外力，且线的密度均匀等于，我们可以取，将一维波动方程化简为：2、一维波动方程的差分解法 在给定的初始条件： 以及给定的边界条件： 我们可以使用差分的方法将一维波动方程写为两种差分格式： 第一种差分格式： 第二种差分格式： 两种差分收敛条件： 2、问题 用一种差分格式计算波动方程混合问题： 3、程序
非线性波动方程：数值解法与实际应用
AI天才研究院
01-07 2077
1.背景介绍 非线性波动方程(Nonlinear Wave Equation, NWE)是一类描述波动过程中非线性现象的部分微分方程。在实际应用中，非线性波动方程广泛用于模拟各种自然现象和工程技术，如气动、流动、热传导、电磁波、光学等。由于非线性波动方程的解是复杂的，因此需要采用数值解法来求解。 本文将从以下六个方面进行阐述： 背景介绍 核心概念与联系 核心算法原理和具体操作步骤以及数学模...
二维波动方程模拟：使用有限差分法对二维波动方程进行数值求解。-matlab开发
06-01
此代码通过有限差分法求解方板上的二维波浪方程，并绘制二维运动和绝对误差的动画。 为简单起见，所有单位都进行了标准化。 它使用 Courant-Friedrich-Levy 稳定性条件。
matlab求解一维波动方程.rar
04-15
在MATLAB中求解一维波动方程，通常采用数值方法，如有限差分法。以下是一些可能涉及的关键步骤： 1. **离散化**：将连续的一维空间和时间轴离散化为网格点，通过差分公式近似偏导数。 2. **初值条件**：根据问题...
波动方程：用有限差分法求解波动方程-matlab开发
06-01
在本项目中，我们将探讨如何使用MATLAB来数值求解波动方程，特别是通过有限差分法这一常见的数值方法。有限差分法是一种将微分方程转化为代数方程组的方法，它通过在离散点上近似导数来解决连续问题。 首先，我们...
波动方程正演
05-06
基于泰勒展开的时空域差分算子的一维波动方程正演，matlab。
求解多阶时间分波动方程的数值方法
05-19
研究论文 求解多阶时间分波动方程的数值方法 Fawang Liu1，Mark M. Meerschaert2， Robert J. McGough3，庄平辉4，刘青霞5
波动方程正演程序
08-30
波动方程正演程序，模拟地震波在地下介质中的传播特征。
二维伪谱法波动方程正演模拟
03-17
用伪谱法模拟波动方程，matlab程序，可以作为波动方程模拟入门学习
MATLAB波动方程数值求解方法
在二维空间中，矩形域内的波动方程通常是一个偏微分方程，其在数值分析领域中的求解通常使用有限差分法、有限元法或谱方法等数值计算技术。MATLAB是一个广泛应用于工程计算、数据分析、算法开发的数学软件，它提供了...
计算海洋声学之波动方程求解思路
oceanljb的博客
02-26 3522
摘要 本文按照计算海洋声学中常见计算流程对波动方程的求解过程进行了梳理，给出了求解要点。本文侧重于整体思路，对求解过程细节之处并未详细讨论。 引言 计算海洋声学的研究内容：海洋声学问题转变为特定边界条件下波动方程的求解问题，利用数值计算手段予以解决。 波动方程的求解思路：通过层层分解，降低波动方程的维度，得到一维波动方程后予以求解，最后通过反变换还原解。 需要关注的要点：方程的维度，边...
波动方程数值模拟
qq_35814924的博客
06-04 6551
波动方程数值模拟一、波动方程简介二、有限差分法三、代码及运行结果四、总结 一、波动方程简介 在实际工作中，经常涉及求解常微分方程边值问题、偏微分方程边值问题。通常，方程的精确解难以求解，而只能得到求解区域内某些点处解的近似值，即数值解。在数学上，二阶线性偏微分方程分为双曲型、抛物型、椭圆型三种，其中双曲型方程的典型代表是波动方程。当满足合适的边界条件和初始条件时，波动方程是可以求解的。常用的求解方...
波动方程（将麦克斯韦方程组求出只有E或H的表达式）
最新发布
m0_53271604的博客
10-01 1315
通过E的波动方程求解E，带入到麦克斯韦中，求解H。或者通过H的波动方程求解H，带入到麦克斯韦中求解E。电场E的波动方程可用来求解电场，求解后带入到麦克斯韦方程组可求解H。3、再根据场的复数表示，得到复数表示的波动方程。复数形式波动方程中场是复数表示的，是r的函数。2、根据麦克斯韦方程组，获得只有E的表达式。先写出线性介质的麦克斯韦方程组。其中k有物理意义，波数。
评论 6
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值