Matlab中显式欧拉法与二阶Runge-Kutta求解ODE系统问题咨询

让我们一步步拆解你的问题，先从方程组本身说起，再分析数值方法的特性，最后给出可对比的实现代码帮你排查问题：

当a=0、b=1时，你的方程组简化为：

x' = -y
y' = x

结合初值x(0)=0、y(0)=1，精确解是x(t) = sin(t)，y(t) = cos(t)——这是一个标准的单位圆周运动，周期为2π，非常适合用来验证数值方法的稳定性和精度。

显式欧拉法的局限性 #

显式欧拉是一阶数值方法，全局截断误差为O(h)，而且它不满足哈密顿系统的辛条件。对于这种周期性的圆周运动，欧拉法的误差会随时间累积，导致轨迹逐渐偏离圆周（甚至向外发散）：

• h越大，每一步的误差越大，累积起来偏离就越明显；
• h过小虽然能降低误差，但会大幅增加计算量。

步长h与总步数N的关联 #

你提到的N应该是总模拟步数，它和步长h的关系是 N = 总模拟时间T / h——两者是强关联的，不能独立调整。如果你的代码里把h和N当成独立参数设置，很可能会导致模拟时间不一致，或者步长计算错误，这也会让结果看起来“过度依赖”这两个参数。

二阶RK法（比如Heun法或中点法）的全局误差为O(h²)，精度比欧拉法高一阶，对步长的敏感性会显著降低。即使步长稍大，也能保持较好的轨迹精度。

下面是两种方法的完整可运行代码，你可以和自己的代码对比，排查问题：

显式欧拉法实现 #

% 参数配置
a = 0;
b = 1;
x0 = 0;
y0 = 1;
T = 2*pi; % 模拟一个完整周期
h_options = [0.1, 0.01, 0.001]; % 不同步长测试

figure;
hold on;
title('显式欧拉法：不同步长结果对比');
xlabel('x');
ylabel('y');
grid on;

for h = h_options
    N = round(T / h);
    t = linspace(0, T, N+1);
    x = zeros(1, N+1);
    y = zeros(1, N+1);
    x(1) = x0;
    y(1) = y0;
    
    for n = 1:N
        % 计算当前点的导数
        dx_dt = a*x(n) - b*y(n);
        dy_dt = b*x(n) + a*y(n);
        % 欧拉法更新
        x(n+1) = x(n) + h*dx_dt;
        y(n+1) = y(n) + h*dy_dt;
    end
    
    plot(x, y, 'DisplayName', ['步长h=', num2str(h)]);
end
% 绘制精确解作为参考
t_exact = linspace(0, T, 1000);
x_exact = sin(t_exact);
y_exact = cos(t_exact);
plot(x_exact, y_exact, 'k--', 'DisplayName', '精确解');
legend;
hold off;

二阶Runge-Kutta（Heun法）实现 #

% 参数配置
a = 0;
b = 1;
x0 = 0;
y0 = 1;
T = 2*pi;
h_options = [0.1, 0.01, 0.001];

figure;
hold on;
title('二阶RK(Heun法)：不同步长结果对比');
xlabel('x');
ylabel('y');
grid on;

for h = h_options
    N = round(T / h);
    t = linspace(0, T, N+1);
    x = zeros(1, N+1);
    y = zeros(1, N+1);
    x(1) = x0;
    y(1) = y0;
    
    for n = 1:N
        % 第一步：欧拉预测步
        dx1 = a*x(n) - b*y(n);
        dy1 = b*x(n) + a*y(n);
        x_pred = x(n) + h*dx1;
        y_pred = y(n) + h*dy1;
        
        % 第二步：校正步
        dx2 = a*x_pred - b*y_pred;
        dy2 = b*x_pred + a*y_pred;
        x(n+1) = x(n) + h*(dx1 + dx2)/2;
        y(n+1) = y(n) + h*(dy1 + dy2)/2;
    end
    
    plot(x, y, 'DisplayName', ['步长h=', num2str(h)]);
end
% 绘制精确解作为参考
t_exact = linspace(0, T, 1000);
x_exact = sin(t_exact);
y_exact = cos(t_exact);
plot(x_exact, y_exact, 'k--', 'DisplayName', '精确解');
legend;
hold off;

对照上面的代码，你可以检查这几个常见错误点：

• 步长与N的计算错误：确保N是总模拟时间T除以h，而不是随便设置的数值；
• 导数表达式错误：比如把x'=-y写成x'=y，这会直接导致轨迹完全错误；
• 数组索引错误：Matlab数组从1开始，确保更新x(n+1)时没有越界，初值赋值正确；
• 未对比精确解：没有精确解作为参考的话，很难判断误差是否在合理范围内。

运行上面的代码你会看到：

• 显式欧拉法在h=0.1时，轨迹已经明显偏离圆周；h=0.01时接近精确解；h=0.001时几乎重合；
• 二阶RK法在h=0.1时就已经和精确解非常接近，步长敏感性远低于欧拉法。

内容的提问来源于stack exchange，提问作者macck