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

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

一、实现效果

二、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_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

三、python代码分享

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from PIL import Image

class WaveEquationFD:
    
    def __init__(self, N, D, Mx, My):
        self.N = N
        self.D = D
        self.Mx = Mx
        self.My = My
        self.tend = 6
        self.xmin = 0
        self.xmax = 2
        self.ymin = 0
        self.ymax = 2
        self.initialization()
        self.eqnApprox()
        
        
    def initialization(self):
        self.dx = (self.xmax - self.xmin)/self.Mx
        self.dy = (self.ymax - self.ymin)/self.My
        
        self.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: 0
        
        self.dt = (self.tend - 0)/self.N
        self.t = np.arange(0, self.tend+self.dt/2, self.dt)
        
        # Assertion for the condition of r < 1, for stability
        r = 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**2
        self.ry = self.D*self.dt**2/self.dy**2
        self.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 = None
        count=1
        k = 0
        nsteps = self.N
        
        while 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.5
            count+=1
    
    
def 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