RK4 代码未产生预期输出

Question

我正在尝试编写自己的 ode4 方法来求解微分方程组，即：

x′ = σy x(t0 = 0) = 0
y′ = ρy + βx y(t0 = 0) = 0.01

代码没有产生预期的输出，我看不出错误来自哪里。 在这种情况下，tout = tspan。 我转置了 k 值，因为尺寸不匹配，并且如果没有它们，代码将无法工作。 我已经绘制了使用另外 2 种方法 ode45 和 ode23 应该获得的结果。

这是代码（2 个文件，一个用于参数和系统，一个用于 ode4 函数）。

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
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('Solution of a system of ODEs by Runge-Kutta methods of different order','fontsize',24);
    xlabel('x',"fontsize", 22);
    ylabel('y',"fontsize", 22);
    legend("4th order RK (ode4, my own)", "3rd order (ode23)","5th order (ode45)","fontsize",18)

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

这是输出：

Answer 1

Matlab/Octave 知道一维数组，并将行向量和列向量视为具有 2D 形状的对象。然而，这些 2D 向量也可以作为数组来寻址，对于矩阵，您必须注意使用 Fortran 列优先数据排列。

Matlab/Octave 和 numpy 的传统被设计为“有帮助”。在某些情况下，尺寸不匹配可以通过称为广播的技术来纠正。这允许添加行和列向量，输入向量只是被复制，直到结果矩阵具有兼容的维度。

这就是这里发生的事情。当你纠正了不兼容被暴露为错误的情况时，你留下了类似的不匹配，它们被默默地“纠正”。

回顾一下：

y(i,:)

是行向量，而

odefun

生成列向量，使用

dtdy=zeros(2,1)

显式构造。因此

v(i,:)+(h/2)*k1

通过作为 2x2 矩阵进行广播而得到“纠正”，这至少给出了一个错误的值。

您可以在构造或操作级别更正此问题，但最好保留其他求解器也使用的函数不变，更改为

  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');

或者也许更系统，也就是说，换位更少

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

这给出了视觉上无法区分的数值解的预期图。