写一个二维声波方程正演四阶精度的MATLAB程序

以下是一个二维声波方程正演四阶精度的MATLAB程序：

% 定义常量
nx = 101;               % x轴方向网格数
ny = 101;               % y轴方向网格数
dx = 10;                % x轴方向网格间距
dy = 10;                % y轴方向网格间距
dt = 0.001;             % 时间步长
nt = 1000;              % 时间步数
vp = 1500;              % 声速

% 初始化波场
u = zeros(nx, ny, nt);
un = zeros(nx, ny);
unn = zeros(nx, ny);

% 初始化震源
sx = 50;
sy = 50;
f0 = 100;
Q = 1/(2*pi*f0);
t = (0:dt:(nt-1)*dt);
src = exp(-Q*(t-1/f0).^2);
src = diff(src);
src = [0, src];
src = repmat(src, [nx, ny, 1]);

% 进行正演
for n = 1:nt
    for i = 4:nx-3
        for j = 4:ny-3
            % 计算四阶差分
            laplacex = (-1/12*u(i+2,j,n) + 4/3*u(i+1,j,n) - 5/2*u(i,j,n) + 4/3*u(i-1,j,n) - 1/12*u(i-2,j,n))/dx^2;
            laplacey = (-1/12*u(i,j+2,n) + 4/3*u(i,j+1,n) - 5/2*u(i,j,n) + 4/3*u(i,j-1,n) - 1/12*u(i,j-2,n))/dy^2;
            
            % 更新波场
            u(i,j,n+1) = 2*u(i,j,n) - un(i,j) + vp^2*dt^2*(laplacex + laplacey) + src(i,j,n)*dt^2;
            
            % 更新波场历史值
            unn(i,j) = un(i,j);
            un(i,j) = u(i,j,n);
        end
    end
end

% 绘制波场图像
figure;
for i = 1:nt
    imagesc(u(:,:,i));
    colorbar;
    caxis([-0.01 0.01]);
    title(sprintf('Time Step: %d', i));
    drawnow;
end

该程序使用了四阶差分的方式来计算二维声波方程的正演，同时使用了震源函数来模拟震源的效果。程序中还包含了波场的初始化和历史值的更新等处理过程。你可以根据需要修改程序中的参数来调整计算的精度和效果。