MATLAB中心差分法求解ODE：h系数为何为1而非推导的2？

用户提供的脚本代码 #

clc;
clear;
close all;
% Step size and paramteres
h = 0.0001;
tspan = 0:h:10;
initial_condition = 3;

% ode45 solver
[t,y] = ode45(@central_diff_f1,tspan,initial_condition);
plot(t,y);

% using the central difference
[t1,y1] = central_diff_f2(tspan,initial_condition,h);
hold on
plot(t1,y1)
legend('ODE 45','Central-Difference');

用户提供的函数代码 #

function [t1,y1] = central_diff_f2(tspan,initial_condition,h)

y1 = zeros(size(tspan));
y1(1) = initial_condition;
t1 = tspan;
% implementing central difference approximation
for i=2:1:length(tspan)
    y1(i) = y1(i-1) + h*(3*y1(i-1)-6);
end

用户的疑问 #

根据中心差分的数学推导，公式为$y'(t)≈\frac{y(t+h)-y(t-h)}{2h}$，整理可得$y(t+h)≈y(t-h)+2h·y'(t)$，代入待求解的ODE方程$y'=3y-6$后，得到$y(t+h)≈y(t-h)+2h·(3y-6)$。但对比代码中的实现y1(i) = y1(i-1) + h*(3*y1(i-1)-6);，为何将t-h对应为索引i-1，且h的系数取1而非推导得出的2？

解答 #

你当前的代码实现根本不是中心差分法，而是向前欧拉法（Forward Euler Method），这就是两者公式差异的核心原因：

1. 索引与步长的对应逻辑
   向前欧拉法的推导基于一阶向前差分近似导数：$y'(t)≈\frac{y(t+h)-y(t)}{h}$，整理后得到递推公式：$y(t+h) = y(t) + h·y'(t)$。
   在代码中，y1(i)对应$y(t+h)$，y1(i-1)对应$y(t)$，步长是h，所以直接用h作为系数，完全符合向前欧拉法的公式。

2. 中心差分法的真正应用场景
   你推导的中心差分公式，是用两点的函数值近似中间点的导数，这种方法通常用于：

   • 已知中间点导数，求解前后点的函数值（但需要额外的初始条件，比如已知$y(t_0)$和$y(t_1)$才能递推后续值）；
   • 数值微分（用已知函数值计算导数），而非直接用于ODE的初值问题求解。

3. 如果要实现中心差分法求解该ODE
   中心差分法属于多步法，需要至少两个初始值才能启动递推：

   • 先通过向前欧拉法计算出$y(t_1) = y(t_0) + h·(3y(t_0)-6)$；
   • 之后对$i≥3$，使用公式：$y1(i) = y1(i-2) + 2h(3*y1(i-1)-6)$（这里$y1(i-1)$对应中间点$y(t)$，$y1(i)$对应$y(t+h)$，$y1(i-2)$对应$y(t-h)$）。

   但这种方法对于初值问题来说，实际应用中并不如向前欧拉或龙格-库塔法常用，因为需要额外初始值，且稳定性也有特定限制。

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