ODE solving using RK4 method (Predator prey)

Question

Diya Saha 2020-4-3

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/515005-ode-solving-using-rk4-method-predator-prey

评论： Rena Berman 2020-5-14

function Q3
clc
clear all
close all 
 
y = zeros(1,1); 
x = zeros(1,1);
%Initial conditions 
t = 0.0; 
x(1) = 2; 
y(1) = 1;
%Parameters
h = 0.001; 
tfinal = 30.0; 
nsteps=(tfinal-t)/h; 
 
%store intial values 
tout(1) = t; 
xout(1) = x(1); 
yout(1) = y(1); 
 
for i=2:nsteps
    k1 = getf(t,y,x); 
    k2 = getf(t+h/2,y+h/2*k1,x+h/2*k1);
    k3 = getf(t+h/2,y+h/2*k2, x+h/2*k2);
    k4 = getf(t+h,y+h*k3, x+h/2*k3);
    y = y+h/6*(k1+2*k2 + 2*k3 + k4);
    x = x+h/6*(k1+2*k2 + 2*k3 + k4);
 
    t = t+h; 
    xout(i) = x(1); 
    yout(i) = y(1); 
    tout(i) = t; 
end 
 
plot(tout,yout,"x",tout,xout,"-r")
 
function f = getf(t,y,x)
alpha = 1.2;
beta = 0.6;
gamma = 0.8;
delta = 0.3;
f(1) = alpha*x(1)-beta*x(1)*y(1); 
f(2) = -gamma*x(1)+delta*y(1)*x(1);

There is something wrong with this because the solution isnt right

5 个评论
显示 3更早的评论隐藏 3更早的评论

Diya Saha 2020-4-3

@Stephen thankyou so much for doing that.

Rena Berman 2020-5-14

(Answers Dev) Restored edit

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

James Tursa 2020-4-3

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/515005-ode-solving-using-rk4-method-predator-prey#answer_423696

在 MATLAB Online 中打开

Your basic problem is that you have two states, x and y, but your function arguments are inconsistent with this. Take this code:

k1 = getf(t,y,x);

getf( ) returns a 2-element vector. So k1 is a 2-element vector. Yet you turn around and treat it like it is a scalar in the very next line:

k2 = getf(t+h/2,y+h/2*k1,x+h/2*k1);

I.e., you end up passing vectors into getf in those 2nd and 3rd input aguments when they should be scalars.

It is unclear to me which element is supposed to be x and which one is y. Since y is first in your argument list, let's assume f(1) is y and f(2) is x. Then the calculation for k2 should be something like this instead:

k2 = getf(t+h/2,y+h/2*k1(1),x+h/2*k1(2));

So all of your function calls will have to be fixed for this.

Having said all of that, I would advise carrying one state vector variable that contains both y and x, instead of carrying two separate y and x variables. That would simplify your code and also put your functions in a form that could be used by calls to the MATLAB supplied integrators such as ode45( ).