二维地震波有限差分走时计算matlab

以下是一个简单的二维地震波有限差分走时计算的matlab程序：

clear all
clc

%定义常数
vp=2000;    %波速
dt=0.002;   %时间步长
dx=10;      %空间步长
nx=100;     %x方向网格数
ny=100;     %y方向网格数
nt=1000;    %时间步数

%初始化地震波场和速度场
u=zeros(nx,ny);
v=zeros(nx,ny);
v(1:nx/2,1:ny/2)=vp;
v(nx/2+1:nx,ny/2+1:ny)=vp;

%计算稳定性限制条件
c=vp*dt/dx;
if c>=1
    fprintf('Stability condition is not satisfied!');
    return;
end

%计算有限差分系数
A=(dt^2)/(dx^2);

%进行时间迭代
for n=1:nt
    %计算波场更新
    for i=2:nx-1
        for j=2:ny-1
            u(i,j)=2*v(i,j)*(1-c^2)*u(i,j)-u(i,j)*(1-c^2)*v(i,j)...
                +A*(v(i+1,j)*u(i+1,j)+v(i,j+1)*u(i,j+1)-4*v(i,j)*u(i,j)+v(i-1,j)*u(i-1,j)+v(i,j-1)*u(i,j-1));
        end
    end
    
    %绘制波场
    imagesc(u);
    colormap(gray);
    drawnow;
end

这段程序中，我们首先定义了常数，包括波速、时间步长、空间步长、网格数等。然后，我们初始化了地震波场和速度场，并计算了稳定性限制条件和有限差分系数。接着，我们进行时间迭代，通过有限差分计算波场更新，并绘制出波场。最终，我们得到了地震波在二维空间中的传播情况。