自制RK4(ode4)求解ODE系统输出异常，寻求错误排查

尝试编写自定义四阶龙格-库塔（RK4）方法求解如下微分方程组：

x′ = σy，x(t₀=0) = 0
y′ = ρy + βx，y(t₀=0) = 0.01

但自制RK4方法输出结果异常：初始值正确，但后续轨迹与Matlab内置的ode45、ode23方法偏差明显，无法定位错误。

主脚本代码 #

clear all
clc

function dtdy=kernel(t,y,sigma,beta,rho)
    dtdy=zeros(2,1);
    dtdy(1)=sigma*y(2);
    dtdy(2)=rho*y(2)+beta*y(1);
end

sigma=1;
rho=-0.9;
beta=-25;
h=0.08;
t=0:h:10;
y0=[0 0.01];

[tRk4,vRk4]=ode4(@(t,y)kernel(t,y,sigma,beta,rho),t,y0);
[t23,v23]=ode23(@(t,y)kernel(t,y,sigma,beta,rho),t,y0);
[t45,v45]=ode45(@(t,y)kernel(t,y,sigma,beta,rho),t,y0);

% 绘图
plot(vRk4(:,1),vRk4(:,2),"b-o",v23(:,1),v23(:,2),'g-^',v45(:,1),v45(:,2),'r-')
set(gca, "linewidth", 0.8, "fontsize", 20)
title('不同阶数龙格-库塔方法求解微分方程组','fontsize',24);
xlabel('x',"fontsize", 22);
ylabel('y',"fontsize", 22);
legend("自制四阶RK(ode4)", "ode23(三阶RK)","ode45(五阶RK)","fontsize",18)

自定义ode4函数代码 #

function[tout,v]=ode4(odefun,tspan,v0)
    v=zeros(length(tspan),length(v0));
    v(1,:)=v0; % 初始条件
    h=tspan(2)-tspan(1); % 步长
    tout=tspan;

    for i=1:length(tspan)-1
        k1=odefun(tspan(i),v(i,:));
        k2=odefun(tspan(i)+(h/2),v(i,:)+(h/2)*k1);
        k3=odefun(tspan(i)+(h/2),v(i,:)+(h/2)*k2);
        k4=odefun(tspan(i)+h,v(i,:)+h*k3);
        v(i+1,:)=v(i,:)+(h/6)*(k1'+(2*k2')+(2*k3')+k4');
    end
end

错误原因分析 #

问题核心是维度不匹配：

• v(i,:)是1×2的行向量，而kernel函数返回的k1/k2/k3/k4是2×1的列向量
• 执行v(i,:)+(h/2)*k1时，Matlab自动广播维度导致计算逻辑错误，误差不断积累最终使轨迹偏离

修正后的ode4函数 #

将状态向量的维度统一为列向量计算，避免广播错误：

function[tout,v]=ode4(odefun,tspan,v0)
    v=zeros(length(tspan),length(v0));
    v(1,:)=v0; % 初始条件
    h=tspan(2)-tspan(1); % 步长
    tout=tspan;

    for i=1:length(tspan)-1
        % 将当前状态转为列向量，与k的维度统一
        y_current = v(i,:)';
        k1=odefun(tspan(i),y_current);
        k2=odefun(tspan(i)+h/2, y_current + h/2*k1);
        k3=odefun(tspan(i)+h/2, y_current + h/2*k2);
        k4=odefun(tspan(i)+h, y_current + h*k3);
        % 更新后的状态转回行向量存入v
        v(i+1,:) = y_current' + (h/6)*(k1 + 2*k2 + 2*k3 + k4)';
    end
end

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