用 matlab 做一个水平范围为500m×500m的正方形湖泊，湖泊的水深为10m，湖水的密度均匀。考虑在湖泊的中心处有一点源初始扰动，即初始时刻湖泊中心处有一1m的水面高程（初始时刻中心处水位为1m），这种由点源扰动产生的波称为圆形波。应用数值模拟方法，对该湖泊中的浅水波传播进行模拟。 由于该水域为浅水，故采用的物理方程为水平二维浅水方程，如方程组1所示，其中， 为水平向速度分量； 为水位变化项； 为海水密度； 为瞬时水深； 为重力加速度。其中，空间步长10m，时间步长0.1s，的程序

以下是一个使用MATLAB进行浅水波传播模拟的示例程序：

% 湖泊参数
L = 500; % 湖泊的水平范围，单位：m
H = 10; % 湖泊的水深，单位：m
dx = 10; % 空间步长，单位：m

% 时间参数
dt = 0.1; % 时间步长，单位：s
t_end = 10; % 模拟的总时间，单位：s

% 网格设置
nx = L/dx + 1; % x方向上的网格数
ny = L/dx + 1; % y方向上的网格数

% 初始条件
eta = zeros(nx, ny); % 水位变化项，初始时刻为0
eta(nx/2, ny/2) = 1; % 中心处水位为1m

% 模拟循环
for t = 0:dt:t_end
    % 计算水平向速度分量和垂直向速度分量（假设为0）
    u = zeros(nx, ny);
    v = zeros(nx, ny);
    
    % 更新水位变化项
    for i = 2:nx-1
        for j = 2:ny-1
            deta_dx = (eta(i+1, j) - eta(i-1, j))/(2*dx);
            deta_dy = (eta(i, j+1) - eta(i, j-1))/(2*dx);
            
            eta(i, j) = eta(i, j) - dt*(u(i, j)*deta_dx + v(i, j)*deta_dy);
        end
    end
    
    % 绘制当前时刻的水位图
    x = 0:dx:L;
    y = 0:dx:L;
    [X, Y] = meshgrid(x, y);
    
    figure;
    surf(X, Y, eta');
    title(['t = ', num2str(t)]);
    xlabel('x (m)');
    ylabel('y (m)');
    zlabel('Water Surface Elevation (m)');
    axis([0 L 0 L -1 1]); % 设置坐标轴范围
    view(2); % 设置视角为俯视图
    colorbar; % 显示颜色条
end

此程序使用有限差分法对浅水方程进行离散化，并使用前向差分法进行时间推进。程序中使用了一个简单的二维数组 `eta` 来表示湖泊中的水位变化项，每个元素对应一个网格点的水位高程。程序通过循环迭代计算每个时间步长内的水位变化，并在每个时间步长结束后绘制当前时刻的水位图。

请注意，这只是一个简单的示例程序，可能需要根据具体情况进行修改和优化，例如边界条件的处理、数值稳定性的考虑等。另外，程序中的绘图部分使用了MATLAB的内置函数 `surf`，你可能需要根据自己的需要进行调整和修改。

希望这能帮助到你开始进行浅水波传播模拟。如有任何问题，请随时向我提问。