有限差分法在matlab中绘制0-100km内地震传播动图，一维波动传播模拟

以下是一维波动传播模拟的 MATLAB 代码，使用有限差分法(Finite Difference Method)来模拟地震波在地下的传播过程。

% 定义参数
nx = 201; % 空间网格数
nt = 1001; % 时间网格数
dt = 0.001; % 时间步长
dx = 500; % 空间步长
vp0 = 2000; % 初始速度
vp = ones(nx,1)*vp0; % 每个网格点上的速度vp
vs0 = 1000; % 初始剪切速度
vs = ones(nx,1)*vs0; % 每个网格点上的速度vs
rho0 = 2000; % 初始密度
rho = ones(nx,1)*rho0; % 每个网格点上的密度rho
f0 = 25; % 激发频率
t0 = 4/f0; % 激发时间
isrc = round(nx/2); % 激发位置
ir = isrc+1; % 接收位置

% 初始化波场
u = zeros(nx,1);
uold = u;
unew = u;
v = u;
vold = v;
vnew = v;

% 计算系数
for i=2:nx
    c1(i) = vp(i)^2 * rho(i);
    c2(i) = (vs(i)^2 - vp(i)^2) * rho(i);
end

% 绘制初始速度模型
figure;
plot(vp,1:nx);
title('Velocity Model');
xlabel('Velocity (m/s)');
ylabel('Depth (m)');
set(gca,'YDir','reverse');
ylim([0 nx*dx]);
xlim([0 3500]);

% 开始模拟
for it=1:nt
    % 计算u, v的时间导数
    for i=2:nx-1
        uxx = u(i+1)-2*u(i)+u(i-1);
        vxx = v(i+1)-2*v(i)+v(i-1);
        unew(i) = 2*u(i)-uold(i) + c1(i)*dt^2/dx^2*uxx + c2(i)*dt^2/dx^2*vxx;
        vnew(i) = 2*v(i)-vold(i) + c1(i)*dt^2/dx^2*vxx + c2(i)*dt^2/dx^2*uxx;
    end
    
    % 边界条件
    unew(1) = 0; unew(nx) = 0;
    vnew(1) = 0; vnew(nx) = 0;
    
    % 添加震源
    if it*dt<t0
        unew(isrc) = unew(isrc) + sin(2*pi*f0*it*dt);
    end
    
    % 更新波场
    uold = u; u = unew;
    vold = v; v = vnew;
    
    % 绘制波场动画
    if mod(it,10)==0
        plot(u,1:nx);
        title(['Seismic Wave Propagation at ',num2str(it*dt),'s']);
        xlabel('Displacement (m)');
        ylabel('Depth (m)');
        set(gca,'YDir','reverse');
        ylim([0 nx*dx]);
        xlim([-0.02 0.02]);
        drawnow;
    end
    
    % 计算接收信号
    if it*dt>t0
        seism(ir,it) = u(ir);
    end
end

% 绘制接收信号
figure;
plot(dt*(1:nt),seism(ir,:));
title('Seismic Record at Receiver');
xlabel('Time (s)');
ylabel('Amplitude');

这段代码模拟了一个长度为 100km，网格数为 201 的介质，初始速度为 2000m/s，初始剪切速度为 1000m/s，初始密度为 2000kg/m3。在距离中心位置为 100km/2=50km 处添加一个频率为 25Hz 的震源，模拟了地震波在介质中的传播过程。

在绘制波场动画的部分，代码每 10 个时间步长绘制一次波场，从而生成了一个 0-100km 内地震传播的动图。接收信号也被记录在了 seism 矩阵中，可以通过绘制 seism(ir,:) 的动态图来观察地震波的传播情况。

需要注意的是，这段代码中使用的是二阶中心差分，因此在计算时间导数时需要使用上一个时间步和下一个时间步的波场。这也是为什么代码中需要定义 uold, unew, vold, vnew 这些变量的原因。