1. Toy模板网首页
  2. > Toy博客

正演(1): 二维声波正演模拟程序（中心差分）Python实现

8月前 作者：苗妮 分类：Toy博客 阅读(30) 违法举报

这篇具有很好参考价值的文章主要介绍了正演(1): 二维声波正演模拟程序（中心差分）Python实现。希望对大家有所帮助。如果存在错误或未考虑完全的地方，请大家不吝赐教，您也可以点击"举报违法"按钮提交疑问。

目录

1、原理： 

1）二维声波波动方程:

​编辑

2）收敛条件（不是很明白）

3）雷克子波

4）二维空间衰减函数

 5）边界吸收条件 (不是很明白。。)

 2、编程实现

1）参数设置：

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

3）边界吸收条件

4）波动方程，迭代公式：

5）全部代码如下：

3、基于matlab的二维波动方程实现 

波动方程数值解是波动方程正演、逆时偏移和全波形反演的核心技术之一。本文采用二阶有限差分对波动方程进行了离散，进而实现了对波动方程的数值求解，模拟出其在介质中的传播过程。

NumPy 通常与 SciPy（Scientific Python）和 Matplotlib（绘图库）一起使用， 这种组合广泛用于替代 MatLab，是一个强大的科学计算环境，有助于我们通过 Python 学习数据科学或者机器学习。

SciPy 包含的模块有最优化、线性代数、积分、插值、特殊函数、快速傅里叶变换、信号处理和图像处理、常微分方程求解和其他科学与工程中常用的计算。

Matplotlib 是 Python 编程语言及其数值数学扩展包 NumPy 的可视化操作界面。它为利用通用的图形用户界面工具包，如 Tkinter, wxPython, Qt 或 GTK+ 向应用程序嵌入式绘图提供了应用程序接口（API）。 

本文代码部分才用了 NumPy和Matplotlib包。

1、原理： 

1）二维声波波动方程:

中心差分法 python代码,地震数据处理,算法,python,机器学习

其中u 为声压，f 为震源中心声压，x/z 为x/z 方向的采样点，t 为时间，v 为速度。

利用泰勒公式对进行展开得到：

中心差分法 python代码,地震数据处理,算法,python,机器学习


两式相加得：

中心差分法 python代码,地震数据处理,算法,python,机器学习


则有：

中心差分法 python代码,地震数据处理,算法,python,机器学习

近似得二阶差分算子：

中心差分法 python代码,地震数据处理,算法,python,机器学习

利用二阶中心差分算子对二阶导数进行离散：

中心差分法 python代码,地震数据处理,算法,python,机器学习

将上式代入声波方程得到二阶中心差分格式：

中心差分法 python代码,地震数据处理,算法,python,机器学习

中心差分法 python代码,地震数据处理,算法,python,机器学习

 其空间和时间差分格式示意图如下图所示：

中心差分法 python代码,地震数据处理,算法,python,机器学习                                          中心差分法 python代码,地震数据处理,算法,python,机器学习

2）收敛条件（不是很明白）

中心差分法 python代码,地震数据处理,算法,python,机器学习

3）雷克子波

雷克子波是地震子波中的一种。由震源激发、经地下传播并被人们在地面或井中接收到的地震波通常是一个短的脉冲振动，称该振动为振动子波， 如下图所示。

中心差分法 python代码,地震数据处理,算法,python,机器学习

公式为： f(t)=(1-2*(pi*f*t)^2)*exp(-(pi*f*t)^2) ，其中t 为时间，f 为主频。

4）二维空间衰减函数

震源中心为1，不衰减，距离中心越远，衰减程度越大。

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))  # 二维空间衰减函数

下图显示的是h_x_z[150] 的曲线，整个的二维空间衰减系数h_x_z，以[150,150] 为中心（震源），向四周衰减。

中心差分法 python代码,地震数据处理,算法,python,机器学习

 5）边界吸收条件 (不是很明白。。)

作用：声波传播时是没有边界的，因此也不存在边界反射问题。但由于模拟正演时，观测范围有限，因此必然是有边界的，边界吸收条件就是尽可能的将能量吸收，将边界反射降到最低。（我的理解哈，欢迎讨论）以下是两种边界吸收条件。现在用的比较多的是PML条件。

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

中心差分法 python代码,地震数据处理,算法,python,机器学习

其中：n 为边界的外法线方向； s 边界的切线方向。
对上式进行离散得到上、下、左、右边界差分格式如下：
中心差分法 python代码,地震数据处理,算法,python,机器学习
中心差分法 python代码,地震数据处理,算法,python,机器学习
其中： 中心差分法 python代码,地震数据处理,算法,python,机器学习   N 、 M 为边界的网格数。

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

中心差分法 python代码,地震数据处理,算法,python,机器学习
对上式进行离散可得上下左右边界计算公式：

中心差分法 python代码,地震数据处理,算法,python,机器学习

中心差分法 python代码,地震数据处理,算法,python,机器学习

 2、编程实现

1）参数设置：

  1. x/z方向长度1500m，x/z方向空间步长5m，每个方向的采样点数为301；
  2. 模拟时长1s，时间步长0.001s，时间采样数1000；
  3. 震源频率25Hz；
  4. 空间衰减因子0.5；
  5. 波速固定，任何位置都为3000m/s
  6. 震源位置在中心；初始声压为0。
# 区域大小
Nx = 301
Nz = 301
# 空间间隔
dx = h
dy = h

# 时间采样数
Nt = 1000
# 时间步
dt = 1 / Nt

# 速度模型
v = np.ones((Nx+1,Nz+1)) * 3000
u = np.zeros((Nt+1,Nx+1,Nz+1))
h = 5

# 子波主频
fm = 25 
# 空间衰减因子
alpha = 0.5

# 迭代公式中的r
A = v **2 * dt ** 2 / h ** 2 
C = v * dt / h

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

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)

3）边界吸收条件

# 边界条件
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 边界条件
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]

4）波动方程，迭代公式：

# 迭代公式
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]

5）全部代码如下：

import numpy as np
import 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]*v[II, JJ]**2*dt**2+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)

参考：

波动方程数值求解（一）_声波方程的解_MaT--2018的博客-CSDN博客文章来源地址https://www.toymoban.com/news/detail-785303.html

3、基于matlab的二维波动方程实现

close all
clear
clc

% 此程序是有限差分法实现声波方程数值模拟


%% 参数设置
delta_t = 0.001; % s
delta_s = 10; % 空间差分：delta_s = delta_x = delta_y (m)
nx = 800;
ny = 800;
nt = 1000;
fmain = 12.5;
%loop:按照10000s为一次大循环；slice代表每隔1000s做一个切片
loop_num = 3;
slice = 1000;
slice_index = 1;
%% 初始化
%震源
t = 1:nt;
t0 = 50;
s_t = (1-2*(pi*fmain*delta_t*(t-t0)).^2).*exp(-(pi*fmain*delta_t*(t-t0)).^2);%源
%一个loop代表向后计算10000s,目的是减少内存消耗
%间隔1000s保存一张切片
num_slice = nt*loop_num/slice;
U_loop(1:ny,1:nx,1:num_slice) = 0;
U_next_loop(1:ny,1:nx,1:2) = 0;
%初始化数组变量
w(1:ny,1:nx) = 0;
U(1:ny,1:nx,1:nt) = 0;
w(400,400) = 1;
V(1:ny,1:nx) = 2000;
A = V.^2*delta_t^2/delta_s^2;
B = V*delta_t/delta_s;

%% 开始计算
JJ = 2:ny-1;
II = 2:nx-1;
start_time = clock;
for loop = 1:loop_num
    fprintf('Loop=%d\n',loop)
    for i_t = 2:nt-1
        if(loop>1)
            s_t(i_t) = 0;
        end
        %上边界
        U(1,II,i_t+1) = (2-2*B(1,II)-A(1,II)).*U(1,II,i_t)+2*B(1,II).*(1+B(1,II)).*U(2,II,i_t)...
            -A(1,II).*U(3,II,i_t)+(2*B(1,II)-1).*U(1,II,i_t-1)-2*B(1,II).*U(2,II,i_t-1);
        %下边界
        U(ny,II,i_t+1) = (2-2*B(ny,II)-A(ny,II)).*U(ny,II,i_t)+2*B(ny,II).*(1+B(ny,II)).*U(ny-1,II,i_t)...
            -A(ny,II).*U(ny-2,II,i_t)+(2*B(ny,II)-1).*U(ny,II,i_t-1)-2*B(1,II).*U(ny-1,II,i_t-1);
        %左边界
        U(JJ,1,i_t+1) = (2-2*B(JJ,1)-A(JJ,1)).*U(JJ,1,i_t)+2*B(JJ,1).*(1+B(JJ,1)).*U(JJ,1+1,i_t)...
                            -A(JJ,1).*U(JJ,1+2,i_t)+(2*B(JJ,1)-1).*U(JJ,1,i_t-1)-2*B(JJ,1).*U(JJ,1+1,i_t-1);
        %右边界
        U(JJ,nx,i_t+1) = (2-2*B(JJ,nx)-A(JJ,nx)).*U(JJ,nx,i_t)+2*B(JJ,nx).*(1+B(JJ,nx)).*U(JJ,nx-1,i_t)...
                            -A(JJ,nx).*U(JJ,nx-2,i_t)+(2*B(JJ,nx)-1).*U(JJ,nx,i_t-1)-2*B(JJ,nx).*U(JJ,nx-1,i_t-1);
        %递推公式
        U(JJ,II,i_t+1) = s_t(i_t).*w(JJ,II)+A(JJ,II).*(U(JJ,II+1,i_t)+U(JJ,II-1,i_t)+U(JJ+1,II,i_t)+U(JJ-1,II,i_t))+...
                            (2-4*A(JJ,II)).*U(JJ,II,i_t)-U(JJ,II,i_t-1);
        if(mod(i_t,100)==0)
            run_time = etime(clock,start_time);
            fprintf('step=%d,total=%d,累计耗时%.2fs\n',i_t+(loop-1)*nt,nt*loop_num,run_time)
            U_loop(:,:,slice_index) = U(:,:,i_t);
            slice_index = slice_index +1;
        end
    end
    %处理四个角点
    KK = 1:nt;
    U(1,1,KK) = 1/2*(U(1,2,KK)+U(2,1,KK));
    U(1,nx,KK) = 1/2*(U(1,nx-1,KK)+U(2,nx,KK));
    U(ny,1,KK) = 1/2*(U(ny-1,1,KK)+U(ny,2,KK));
    U(ny,nx,KK) = 1/2*(U(ny-1,nx,KK)+U(ny,nx-1,KK));
    %% 为下一次loop做准备
    fprintf('step=%d,total=%d,累计耗时%.2fs\n',i_t+1+(loop-1)*nt,nt*loop_num,run_time)
    U_loop(:,:,slice_index) = U(:,:,nt);
    slice_index = slice_index +1;
    U_next_loop(:,:,1)=U(:,:,nt-1);
    U_next_loop(:,:,2)=U(:,:,nt);
    U(:,:,:) = 0;
    U(:,:,1) = U_next_loop(:,:,1);
    U(:,:,2) = U_next_loop(:,:,2);
end

%% 制作动图
fmat=moviein(num_slice);
filename = 'FDM_4_homogenerous.gif';
for II = 1:num_slice
    pcolor(U_loop(:,:,II));
    shading interp;
    axis tight;
    set(gca,'yDir','reverse');
    str_title = ['FDM-4-homogenerous  t=',num2str(delta_t*II*100),'s'];
    title(str_title)
    drawnow; %刷新屏幕
    F = getframe(gcf);%捕获图窗作为影片帧
    I = frame2im(F); %返回图像数据
    [I, map] = rgb2ind(I, 256); %将rgb转换成索引图像
    if II == 1
        imwrite(I,map, filename,'gif', 'Loopcount',inf,'DelayTime',0.1);
    else
        imwrite(I,map, filename,'gif','WriteMode','append','DelayTime',0.1);
    end
    fmat(:,II)=getframe;
end
movie(fmat,10,5);
%% 绘图为gif
pcolor(U_loop(:,:,num_slice))
shading interp;
axis tight;
set(gca,'yDir','reverse');
str_title = ['FDM-4-homogenerous  t=',num2str(delta_t*num_slice*100),'s'];
title(str_title)
colormap('Gray')
filename = [str_title,'.jpg'];
saveas(gcf,filename)

%% 耗时
toc

到了这里，关于正演(1): 二维声波正演模拟程序（中心差分）Python实现的文章就介绍完了。如果您还想了解更多内容，请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章，希望大家以后多多支持TOY模板网！

本文来自互联网用户投稿，该文观点仅代表作者本人，不代表本站立场。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如若转载，请注明出处： 如若内容造成侵权/违法违规/事实不符，请点击违法举报进行投诉反馈，一经查实，立即删除！

分享到：
领支付宝红包 赞助服务器费用

相关文章

  • 二维差分详解

    二维差分详解

    前言 上一期我们分享了一维差分的使用方法，这一期我们将接着上期的内容带大家了解二位差分的使用方法，话不多说，LET’S GO!（上一期链接） 二维差分 二维差分我们可以用于对矩阵区间进行多次操作的题。 二维差分我们还可以采用一维差分的思想，如图假如我们要对区

    2024年02月03日
    浏览(36)

  • 二维差分算法详解

    二维差分模板 给定一个n行m列的矩阵，下标从1开始。接下来有q次操作，每次操作输入5个参数x1, y1, x2, y2, k，表示把以(x1, y1)为左上角,(x2,y2)为右下角的子矩阵的每个元素都加上k。请输出操作后的矩阵。 第一行包含三个整数n,m,q。接下来n行，每行m个整数，代表矩阵的元素。接

    2024年01月20日
    浏览(27)
  • 洛谷 P3397 地毯 刷题笔记 二维差分矩阵

    洛谷 P3397 地毯 刷题笔记 二维差分矩阵

    P3397 地毯 - 洛谷 | 计算机科学教育新生态 (luogu.com.cn) 方法1 欺负数据小  暴力水过 #includeiostream using namespace std; const int N=1010; int a[N][N]; int main(){     int n,m;     cinnm;     for(int i=0;im;i++){         int x1,y1,x2,y2;         cinx1y1x2y2;         for(int q=x1;q=x2;q++){           

    2024年02月03日
    浏览(59)
  • 微信小程序在模拟器可以显示，手机扫编译二维码显示不了 解决办法

    微信小程序在模拟器可以显示，手机扫编译二维码显示不了 解决办法

    1，本地服务器即是使用localhost，如果希望在手机上可以演示，请把微信开发者工具代码中wx.request请求数据的所有带有localhost的路径全部替换本机的ip地址，如http://localhost:8080/test/hhh, 改成http://172.00.00.000/test/hhh, 这里只是举个例子，查看ip地址方法： 打开cmd，输入ipconfig，红色

    2024年02月12日
    浏览(28)

  • 数学模型：Python实现差分方程

    上篇文章：微分方程 文章摘要：差分方程的Python实现。 参考书籍：数学建模算法与应用(第3版)司守奎 孙玺菁。 PS1：只涉及了具体实现并不涉及底层理论。没有给出底层理论参考书籍的原因是不想做这个方向吧。所以对我只要掌握基本模型有个概念那就好了。 PS2：这里跳过

    2024年02月16日
    浏览(25)
  • 利用中心差分格式求解一阶波动方程(附Matlab代码)

    利用中心差分格式求解一阶波动方程(附Matlab代码)

    ∂ 2 u ∂ t 2 − ∂ 2 u ∂ x 2 = 0 , 0 x 1 , t 0 , frac{partial^2u}{partial t^2}-frac{partial^2u}{partial x^2}=0,0x1,t0, ∂ t 2 ∂ 2 u ​ − ∂ x 2 ∂ 2 u ​ = 0 , 0 x 1 , t 0 , 初始边值条件为: u ( 0 , x ) = 2 s i n ( π x ) , uleft(0,xright)=2sinleft(pi xright), u ( 0 , x ) = 2 s in ( π x ) , u t ( 0 , x ) = 0 , 0 x 1 , u_tle

    2024年02月06日
    浏览(65)
  • 高校数据中心网络规划设计及思科模拟器CISCO模拟实现（网络安全、数据冗余）

    高校数据中心网络规划设计及思科模拟器CISCO模拟实现（网络安全、数据冗余）

            数据中心在现代社会中的地位愈加重要，这得益于信息技术的迅速发展。信息处理的能力、安全性等方面的要求也在不断攀升。因此，在服务器的计算能力、稳定性、可靠性、安全性、未来扩展性以及方便管理等多个方面，都应对其要求更高水平。         高校没有

    2024年02月01日
    浏览(45)

  • 有限差分法-一维热传导方程及其Matlab程序实现

    2.2.2 一维热传导方程 热传导方程是描述热量在介质中传导的数学模型。在许多实际应用中，我们需要预测物体的温度随时间和空间的演化情况，这就需要用到热传导方程。 热传导方程的背景可以追溯到18世纪，当时科学家们对热的本质和热量如何传递产生了浓厚的兴趣。傅里

    2024年02月10日
    浏览(30)
  • 【嵌入式开发】基于树莓派实现超声波避障小车（Python）

    【嵌入式开发】基于树莓派实现超声波避障小车（Python）

    1.1 所需硬件 （1）烧制好的树莓派4B （2）小车车架（可在网上购买）。 （3）直流电机*4：用于驱动小车行驶。 （4）L298N电机驱动模块：用于实现对电机的控制。 （5）超声波测距模块：用于实时测距，以实现自主避障。 （6）其余辅助器件：包括充电宝（树莓派供电）、干电

    2023年04月14日
    浏览(57)

  • 周赛354（模拟、差分、排序+双指针、枚举）

    难度简单1 给你一个下标从 1 开始、长度为 n 的整数数组 nums 。 对 nums 中的元素 nums[i] 而言，如果 n 能够被 i 整除，即 n % i == 0 ，则认为 num[i] 是一个 特殊元素 。 返回 nums 中所有 特殊元素 的 平方和 。 示例 1： 示例 2： 提示： 1 = nums.length == n = 50 1 = nums[i] = 50 模拟 python 难

    2024年02月16日
    浏览(30)

觉得文章有用就打赏一下文章作者

支付宝扫一扫打赏

博客赞助

微信扫一扫打赏

请作者喝杯咖啡吧~博客赞助

支付宝扫一扫领取红包，优惠每天领

二维码1

领取红包

二维码2

领红包