有限差分模拟二维波动方程

有限差分法是一种数值计算方法，用于近似连续函数的导数，常用于解决偏微分方程，如波动方程。在二维波动方程中，通常描述的是在二维空间中的波动现象，例如声波、电磁波或水波的传播。

二维波动方程的标准形式是：

\[ \frac{\partial^2 u}{\partial t^2} = c^2(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}) \]

其中，\(u(x, y, t)\) 是位移或扰动，\(c\) 是波速，\(t\) 是时间。

有限差分方法是通过将时间和空间上的偏导数替换为它们在网格点处的离散近似来求解这个方程。常用的方法有向前差分（用当前值估计未来的值）、向后差分（用过去的值估计未来的值）以及中心差分（使用前后两点的值进行平均），这取决于阶数的选择。

以下是一个简单的二维有限差分示例，用四阶中心差分方法（即二阶对时间，二阶对空间）来逼近方程：

% 定义网格参数
dx = 0.1; % 空间步长
dt = 0.01; % 时间步长
Lx = 10; % 横向长度
Ly = 10; % 纵向长度
N = Lx/dx; % 空间网格点数
T = 5; % 总时间

% 初始化二维数组存储波场
u = zeros(N, N, T+1);

% 边界条件（例如固定边界或零初始条件）
u(:, 1, :) = 0;
u(:, N, :) = 0;
u(1, :, :) = 0;
u(N, :, :) = 0;

% 波速
c = 1;

% 二维中心差分
for t = 1:T
    u(:,:,t+1) = u(:,:,t) + dt/c^2 * (4*u(:,:,t) - u(:,:,t-1) + ...
        central_diff2(u, 2, dx, 'y') + central_diff2(u, 2, dx, 'x'));
end

% 二维中心差分函数（对于x和y方向）
central_diff2(u, k, dx, direction) =
    if strcmp(direction, 'x')
        (u(2:end-1, :, t) - 2*u(1:end-2, :, t))/dx^2
    else
        (u(:, 2:end-1, t) - 2*u(:, 1:end-2, t))/dx^2
    end;

在这个例子中，`central_diff2`函数用于计算空间偏导数，`direction`参数决定了是沿x轴还是y轴。