四阶龙格-库塔（Runge-Kutta 4）Octave代码报错排查请求

问题背景 #

使用四阶龙格-库塔方法求解一阶微分方程组，Octave环境下代码无法运行，命令窗口无错误提示。调试确认k1、k2、k3、k4计算正常，但执行v(i+1,:)赋值时程序崩溃。

原始代码 #

主脚本 #

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];

Rk4=ode4(@(t,y)kernel(t,y,sigma,beta,rho),t,y0);

ode4函数实现 #

function[tout,v]=ode4(odefun,tspan,v0)
v=zeros(length(tspan),length(v0));
v(1,:)=v0 %initial condition for y(0)
h=tspan(2)-tspan(1); % extract h value from tspan
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

错误原因 #

核心问题是向量维度不匹配：

• kernel函数返回的dtdy是2x1列向量，但主脚本中初始值y0是1x2行向量，ode4函数中v(i,:)也是行向量
• 行向量与列向量进行逐元素运算时，Octave会触发维度不兼容错误，导致v(i+1,:)赋值失败

修复方案 #

统一所有状态向量的维度，推荐使用行向量匹配v的矩阵结构（每行对应一个时间点的状态），同时补充分号避免冗余输出：

修复后主脚本 #

clear all
clc

function dtdy=kernel(t,y,sigma,beta,rho)
  dtdy=zeros(1,2); % 修改为1x2行向量，匹配v的行维度
  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]; % 保持行向量，与kernel返回值维度一致

Rk4=ode4(@(t,y)kernel(t,y,sigma,beta,rho),t,y0);

修复后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

验证说明 #

修复后所有向量维度统一为行向量，迭代过程中运算逻辑兼容，代码可正常运行并输出RK4求解结果。若偏好使用列向量，可将初始值改为y0=[0; 0.01]，同时将kernel返回的dtdy保持列向量，在ode4中修改v(1,:)=v0';即可。

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