二维声波方程的有限差分法数值模拟

文章目录

二维声波方程的有限差分法数值模拟
- 一、实现效果
二、matlab代码分享
三、python代码分享

一、实现效果

二、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 = 100;
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_numfprintf('Loop=%d\n',loop)for i_t = 2:nt-1if(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;endend%处理四个角点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_slicepcolor(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 == 1imwrite(I,map, filename,'gif', 'Loopcount',inf,'DelayTime',0.1);elseimwrite(I,map, filename,'gif','WriteMode','append','DelayTime',0.1);endfmat(:,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

三、python代码分享

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from PIL import Imageclass WaveEquationFD:def __init__(self, N, D, Mx, My):self.N = Nself.D = Dself.Mx = Mxself.My = Myself.tend = 6self.xmin = 0self.xmax = 2self.ymin = 0self.ymax = 2self.initialization()self.eqnApprox()def initialization(self):self.dx = (self.xmax - self.xmin)/self.Mxself.dy = (self.ymax - self.ymin)/self.Myself.x = np.arange(self.xmin, self.xmax+self.dx, self.dx)self.y = np.arange(self.ymin, self.ymax+self.dy, self.dy)#----- Initial condition -----#self.u0 = lambda r, s: 0.2*np.exp(-((r-(self.xmax + self.xmin)/2)**2+(s-(self.ymax + self.ymin)/2)**2)/0.01)#----- Initial velocity -----#self.v0 = lambda a, b: 0#----- Boundary conditions -----#self.bxyt = lambda left, right, time: 0self.dt = (self.tend - 0)/self.Nself.t = np.arange(0, self.tend+self.dt/2, self.dt)# Assertion for the condition of r < 1, for stabilityr = 4*self.D*self.dt**2/(self.dx**2+self.dy**2);assert r < 1, "r is bigger than 1!"def eqnApprox(self):#----- Approximation equation properties -----#self.rx = self.D*self.dt**2/self.dx**2self.ry = self.D*self.dt**2/self.dy**2self.rxy1 = 1 - self.rx - self.ry self.rxy2 = self.rxy1*2#----- Initialization matrix u for solution -----#self.u = np.zeros((self.Mx+1, self.My+1))self.ut = np.zeros((self.Mx+1, self.My+1))self.u_1 = self.u.copy()#----- Fills initial condition and initial velocity -----#for j in range(1, self.Mx):for i in range(1, self.My):self.u[i,j] = self.u0(self.x[i], self.y[j])self.ut[i,j] = self.v0(self.x[i], self.y[j])def solve_and_animate(self):u_2 = np.zeros((self.Mx+1, self.My+1))xx, yy = np.meshgrid(self.x, self.y)fig = plt.figure()        ax = fig.add_subplot(111, projection='3d')wframe = Nonecount=1k = 0nsteps = self.Nwhile k < nsteps:if wframe:ax.collections.remove(wframe)self.t = k*self.dt#----- Fills in boundary condition along y-axis (vertical, columns 0 and Mx) -----#for i in range(self.My+1):self.u[i, 0] = self.bxyt(self.x[0], self.y[i], self.t)self.u[i, self.Mx] = self.bxyt(self.x[self.Mx], self.y[i], self.t)for j in range(self.Mx+1):self.u[0, j] = self.bxyt(self.x[j], self.y[0], self.t)self.u[self.My, j] = self.bxyt(self.x[j], self.y[self.My], self.t)if k == 0:for j in range(1, self.My):for i in range(1, self.Mx):self.u[i,j] = 0.5*(self.rx*(self.u_1[i-1,j] + self.u_1[i+1,j])) \+ 0.5*(self.ry*(self.u_1[i,j-1] + self.u_1[i,j+1])) \+ self.rxy1*self.u[i,j] + self.dt*self.ut[i,j]else:for j in range(1, self.My):for i in range(1, self.Mx):self.u[i,j] = self.rx*(self.u_1[i-1,j] + self.u_1[i+1,j]) \+ self.ry*(self.u_1[i,j-1] + self.u_1[i,j+1]) \+ self.rxy2*self.u[i,j] - u_2[i,j]u_2 = self.u_1.copy()self.u_1 = self.u.copy()wframe = ax.plot_surface(xx, yy, self.u, cmap='rainbow', linewidth=2, antialiased=False)ax.set_xlim3d(0, 2.0)ax.set_ylim3d(0, 2.0)ax.set_zlim3d(-1.5, 1.5)ax.set_xticks([0, 0.5, 1.0, 1.5, 2.0])ax.set_yticks([0, 0.5, 1.0, 1.5, 2.0])ax.set_xlabel("x")ax.set_ylabel("y")ax.set_zlabel("U")plt.savefig(str(count))plt.pause(0.01)k += 0.5count+=1def main():simulator = WaveEquationFD(200, 0.1, 100, 100)simulator.solve_and_animate()plt.show()if __name__ == "__main__":main()#N = 200
#D = 0.25
#Mx = 50
#My = 50
imgs = []

[^1]本文仅做学习交流，详细原理请参考@link@link

二维声波方程的有限差分法数值模拟相关推荐

二维有限元方程matlab,有限元法求解二维Poisson方程的MATLAB实现
有限元法求解二维 Poisson 方程的 MATLAB 实现陈莲a ,郭元辉b ,邹叶童a ( 西华师范大学 a．数学与信息学院; b．教育信息技术中心,四川南充 6437009) 摘要: ...
二维burgers方程_二维Burgers方程的RKDG有限元解法
二维 Burgers 方程的 RKDG 有限元解法 ∗ 马艳春 1, 张寅虎 2, 冯新龙 1 [摘要] 摘要 : 本文应用 RKDG 有限元方法求解具有周期边界条件的二维非粘性 Burgers ...
二维Poisson方程五点差分格式及简单求解方法Python实现
二维Poisson方程简介给出二维 PoissonPoissonPoisson 方程 DirichletDirichletDirichlet 边值问题: {−Δu=f(x,y)(x,y)∈Ωu=φ ...
galerkin有限元法matlab实现,有限元法求解二维Poisson方程的MATLAB实现
有限元法求解二维Poisson方程的MATLAB实现陈莲a,郭元辉b,邹叶童a [摘要]文章讨论了圆形区域上的三角形单元剖分.有限元空间,通过变分形式离散得到有限元方程. 用MATLAB编程求得数值 ...
使用matlab求解二维浅水方程的数值解（一）—浅水波
最近在读<ocean modelling for beginners>这本书,对于做海洋数值模拟工作的小白来说,这绝对是一本好书.强烈推荐给理论基础较弱的学习者,这本书循序渐进,由简入繁的 ...
二维burgers方程_用格子Boltzmann方法研究二维Burgers方程
用格子 Boltzmann 方法研究二维 Burgers 方程张伟 ; 李文杰 [期刊名称] <天津城市建设学院学报> [年 ( 卷 ), 期] 2012(018)001 [ 摘要 ] ...
有限元方法基础-以二维拉普拉斯方程为例(附程序)
文章目录前言问题描述问题区域偏微分方程变分问题(弱形式) 有限元离散二维线性有限元 : 参考基函数 2D linear finite element : affine mapping 有限 ...
二维各向同性介质弹性波数值模拟(交错网格有限差分法)
一.一阶速度-应力弹性波方程在二维各向同介质xoz平面内,假定体力为0. 从上面方程当中,我们为了得到各点的应力和速度值,就需要得到关于对时间t和空间x,z的偏导. 二.时间上的2M阶差分由Tay ...
二维非稳态导热微分方程_二维非稳态传热的温度场数值模拟
背景:这是本学期凝固实验课的实验之一.这节课有两个数值模拟实验,第一个是二维常物性的,只有一种介质.而第二个实验是模拟凝固过程,稍微复杂一些.这篇文章是针对第一个实验写的,实验书上是按照显示差分进行的 ...

二维声波方程的有限差分法数值模拟

二维声波方程的有限差分法数值模拟

文章目录

一、实现效果

二、matlab代码分享

三、python代码分享

二维声波方程的有限差分法数值模拟相关推荐

最新文章

热门文章