各向同性介质二维声波波动方程有限差分法数值模拟（波场模拟动图）

这是一个补充代码，是声波二阶位移方程的有限差分法波场模拟程序，没有交错网格。

弹性波一阶速度应力-应力方程及其交错网格波场模拟见我前面的帖子。C/C++，MATLAB代码的帖子都可以找到。

clc;
clear;
Nx=401; Nz=401; Nt=650;%设置采样点数，采样时间点数
h=5;    %x方向和z方向的步长
dt=0.001;   %时间步长
c=zeros(Nx,Nz);
A=zeros(Nx,Nz);
for i=1:Nx
    for j=1:Nz
        c(i,j)=2000;    %波传播速度为3km/s
        if(j>240)
            c(i,j)=3500;
        end
        f=30;           %震源主频
        t0=1.2/f;
        A(i,j)=(dt*c(i,j))^2/h^2;
%         A(i,j)=1;
    end
end

u=zeros(Nx,Nz,Nt);
for k=2:Nt-1
    for i=3:Nx-2
        for j=3:Nz-2
            if i==200 & j==200 & k<100
                u(i,j,k+1)=exp(-(pi*f*(k*dt-t0)).^2).*(1-2.*(pi*f*(k*dt-t0)).^2);
                %在(100km，100km)处设置一个振动源
            else u(i,j,k+1)=A(i,j)*(u(i+1,j,k)+u(i-1,j,k)+u(i,j+1,k)+u(i,j-1,k)-4*u(i,j,k))-u(i,j,k-1)+2*u(i,j,k);
            end
            u(3,j,k+1)=u(4,j,k+1);
        end
    end
end


filename='二维波场快照.gif';
figure('units','normalized','position',[0,0,0.5,0.75])%设置成图框的大小
for k=1:4:Nt
    aa=u(:,:,k);
    aa=aa';
    pcolor(aa);
    shading interp;
    colormap('gray');
    caxis([-.01 .01]);
    axis equal;
    axis([1,Nx,1,Nz]);
    set(gca,'Ydir','reverse');
    xlabel('x'); ylabel('z');
    title('二维声波传播快照');
    % if(k==201) keyboard; end
    f=getframe(gcf); % 捕获画面
    imind=frame2im(f);
    [imind,cm] = rgb2ind(imind,256);
    if k==1
        imwrite(imind,cm,filename,'gif', 'Loopcount',inf,'DelayTime',0.05); %采用延迟时间为0.05秒写入给定的文件
    else
        imwrite(imind,cm,filename,'gif','WriteMode','append','DelayTime',0.1); %采用延迟时间为0.1秒写入给定的文件
    end
end

来个双层速度的波场快照：

声明：

该代码非本人原创，本人学习之余略有改动。原文请参考二维声波传播方程的有限差分模拟_交错网格法二维声波方程有限差分正演模拟代码-CSDN博客