четвертого порядка методом Рунге-Кутта (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: %f\n', x); 
fprintf('y: %f\n', y); 
fprintf('A: %f\n', A); 
fprintf('Time: %f\n', 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

Похож на форму волны уравнения Лотки-Вольтерра?

Я не уверен, что если ваше начальное условие [300;8000] или [8000;300] (вы указываете его в обоих направлениях выше), но независимо от того, у вас есть колебательная система, которую вы пытаетесь интегрировать с большим фиксированным шагом времени, который (намного) больше периода колебаний. Вот почему ваша ошибка взрывается. Если вы попытаетесь увеличить n (скажем, 1e6), вы обнаружите, что в итоге вы получите стабильное решение (при условии, что ваша реализация Runge-Kutta в противном случае правильна).

Есть ли причина, по которой вы не используете встроенные решатели ODE от Matlab, например. 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 и/или регулировать параметры 'AbsTol' и 'RelTol' с помощью odeset.