四阶龙格-库塔法(RK4)在几次迭代后失效

Question

我想解决:

x' = 60*x - 0.2*x*y;
y' = 0.01*x*y - 100* y;

使用四阶龙格-库塔算法。

起点:x(0) = 8000, y(0) = 300范围:[0,15]

完整的函数如下:

function [xx yy time r] = rk4_m(x,y,step)
A = 0;
B = 15;
h = step;
iteration=0;
t = tic;
xh2 = x;
yh2 = y;

rr = zeros(floor(15/step)-1,1);
xx = zeros(floor(15/step)-1,1);
yy = zeros(floor(15/step)-1,1);
AA = zeros(1, floor(15/step)-1);
while( A < B)

    A = A+h;
    iteration = iteration + 1;
    xx(iteration) = x;
    yy(iteration) = y;
    AA(iteration) = A;
    [x y] = rkstep(x,y,h);

    for h2=0:1
        [xh2 yh2] = rkstep(xh2,yh2,h/2);
    end
    r(iteration)=abs(y-yh2);

end
time = toc(t);
xlabel('Range');
ylabel('Value');     
hold on
plot(AA,xx,'b');
plot(AA,yy,'g');
plot(AA,r,'r');
fprintf('Solution:n');
fprintf('x: %fn', x);
fprintf('y: %fn', y);
fprintf('A: %fn', A);
fprintf('Time: %fn', time);
end
function [xnext, ynext] = rkstep(xcur, ycur, h)
    kx1 = f_prim_x(xcur,ycur);
    ky1 = f_prim_y(xcur,ycur);
    kx2 = f_prim_x(xcur+0.5*h,ycur+0.5*h*kx1);
    kx3 = f_prim_x(xcur+0.5*h,ycur+0.5*h*kx2);
    kx4 = f_prim_x(xcur+h,ycur+h*kx3);
    ky2 = f_prim_y(xcur+0.5*h*ky1,ycur+0.5*h);
    ky3 = f_prim_y(xcur+0.5*h*ky2,ycur+0.5*h);
    ky4 = f_prim_y(xcur+h*ky2,ycur+h);
    xnext = xcur + (1/6)*h*(kx1 + 2*kx2 + 2*kx3 + kx4);
    ynext = ycur + (1/6)*h*(ky1 + 2*ky2 + 2*ky3 + ky4);
end
function [fx] = f_prim_x(x,y)
fx = 60*x - 0.2*x*y;
end
function [fy] = f_prim_y(x,y)
fy = 0.01*x*y - 100*y;
end

我通过执行[xx yy time] = rk4_m(8000,300,10)

来运行它

问题是，在2-3次迭代返回无用的结果后，所有内容都崩溃了。我做错了什么?还是这种方法不适用于这类方程?

有意省略分号。

---

看起来我没有注意实际的h大小。它现在工作了!谢谢!

Answer 1

看起来像是某种形式的Lotka-Volterra方程?

我不确定如果你的初始条件是[300;8000]还是[8000;300](你在上面两种方式都指定了)，但不管怎样，你有一个振荡系统，你试图用一个大的固定时间步长来积分，这个时间步长(远远)大于振荡周期。这就是为什么你的错误会爆发。如果你尝试增加n(比如1e6)，你会发现最终你会得到一个稳定的解决方案(假设你的龙格-库塔实现是正确的)。

是否有理由不使用Matlab的内置ODE求解器，例如ode45或ode15s ?

function ode45demo
[t,y]=odeode45(@f,[0 15],[300;8000]);
figure;
plot(t,y);
function ydot=f(t,y)
ydot(1,1) = 60*y(1) - 0.2*y(1)*y(2);
ydot(2,1) = 0.01*y(1)*y(2) - 100*y(2);

你会发现自适应步长求解器对于这些类型的振荡问题更有效。因为你的系统有如此高的频率，看起来相当僵硬，我建议你也看看ode15s给出的和/或用odeset调整'AbsTol'和'RelTol'选项。

Answer 2

直接的问题是RK4代码没有完全从标量情况进化到两个耦合方程的情况。注意，导数函数中没有时间参数。x和y都是因变量，因此得到每一步导数函数定义的斜率更新。然后xcur得到kx的更新，ycur得到ky的更新。

function [xnext, ynext] = rkstep(xcur, ycur, h)
    kx1 = f_prim_x(xcur,ycur);
    ky1 = f_prim_y(xcur,ycur);
    kx2 = f_prim_x(xcur+0.5*h*kx1,ycur+0.5*h*ky1);
    ky2 = f_prim_y(xcur+0.5*h*kx1,ycur+0.5*h*ky1);
    kx3 = f_prim_x(xcur+0.5*h*kx2,ycur+0.5*h*ky2);
    ky3 = f_prim_y(xcur+0.5*h*kx2,ycur+0.5*h*ky2);
    kx4 = f_prim_x(xcur+h*kx3,ycur+h*ky3);
    ky4 = f_prim_y(xcur+h*kx3,ycur+h*ky3);
    xnext = xcur + (1/6)*h*(kx1 + 2*kx2 + 2*kx3 + kx4);
    ynext = ycur + (1/6)*h*(ky1 + 2*ky2 + 2*ky3 + ky4);
end