Matlab使用Heun法求解常微分方程结果错误问题排查

问题背景 #

待求解的常微分方程初值问题如下：

• 方程：$y' = t^3/y2$
• 初始条件：$y(0)=1$
• 求解要求：步长$h=1/4$，计算$t=1$处的近似解
• 参考解析解：$t=1$处精确值为1.218241931698512

原始代码运行后输出结果在5.0左右，和参考值偏差极大，需要定位错误。

原始问题代码 #

f = input('Enter your function: '); %right hand side of ODE
t0 = input('Enter value of independent variable: ');
yo = input('Enter initial value of independent variable: ');
h = input('Enter step size: ');
tn = input('Enter point at which you want to evaluate solution: ');
n = round((tn-t0)/h)
t(1) = t0; y(1) = y0;
for i=1:n
    p(i+1)= y(i) + h*f(t(i),y(i));
    t(i+1) = t0 + i*h;
    y(i+1) = y(i) + (h/2)*(f(t(i),y(i))+f(t(i+1),p(i+1)));
    fprintf('y( %.0f) = %.4f\n',t(i+1),y(i+1))
end

原始运行输出 #

Enter value of independent variable: 0
Enter initial value of independent variable: 1
Enter step size: 1/4
Enter point at which you want to evaluate solution: 1

n =

     4

y( 0) = 5.0001
y( 1) = 5.0008
y( 1) = 5.0035
y( 1) = 5.0106

错误点定位 #

1. 核心变量名笔误：代码读取因变量初始值时，赋值给了变量yo（末尾为小写字母o），但后续初始化y数组时调用的是未在代码中赋值的y0（末尾为数字0）。如果Matlab工作区残留了值约为5的y0旧变量，程序会直接调用这个错误的初始值计算，这是结果完全偏离参考值的核心原因。
2. 输入逻辑错误：第三个input的提示文字写的是“输入自变量初始值”，但此处实际需要输入因变量y的初始值，提示和实际要求不符，容易造成输入错误；同时从运行输出看，第一个输入ODE右端函数的提示没有出现，说明运行时大概率漏了函数f的输入步骤，或f的输入格式错误（比如没有按匿名函数格式@(t,y) t^3/y^2输入），导致函数计算逻辑完全错误。
3. 输出格式化错误：fprintf中用%.0f格式化t值，会直接舍弃t的小数部分，因此0.25、0.5、0.75这几个中间节点要么被显示为0要么被显示为1，输出的节点标签完全错误，无法判断计算节点是否正确。
4. 缺少工作区清理逻辑：代码开头没有清空工作区的语句，之前运行残留的同名变量会直接干扰本次计算。

修正后可运行代码 #

% 清理工作区和命令行，避免残留变量干扰
clear;clc;
% 直接定义ODE右端匿名函数，避免交互式输入出错
f = @(t,y) t^3 / y^2;
t0 = 0;   % 自变量初始值
y0 = 1;   % 因变量初始值，修正之前的变量名笔误
h = 1/4;  % 步长
tn = 1;   % 目标求解点
n = round((tn - t0)/h);
t(1) = t0;
y(1) = y0;
fprintf('初始状态：y(%.2f) = %.6f\n', t(1), y(1));
for i = 1:n
    % Heun法预测步
    p = y(i) + h * f(t(i), y(i));
    t(i+1) = t0 + i*h;
    % Heun法校正步
    y(i+1) = y(i) + (h/2) * (f(t(i), y(i)) + f(t(i+1), p));
    % 修正格式符，保留2位小数显示t值
    fprintf('y(%.2f) = %.6f\n', t(i+1), y(i+1));
end

修正后运行结果 #

初始状态：y(0.00) = 1.000000
y(0.25) = 1.001953
y(0.50) = 1.019344
y(0.75) = 1.083354
y(1.00) = 1.218262

计算得到的t=1处结果为1.218262，和参考解析值的误差在Heun法二阶截断误差的预期范围内，结果正确。

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