MATLAB Answers

ODE solving using RK4 method (Predator prey)

조회 수: 9(최근 30일)
Diya Saha
Diya Saha 3 Apr 2020
댓글: Rena Berman 14 May 2020 16:10
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

표시 이전 댓글 수: 2
Stephen Cobeldick
Stephen Cobeldick 3 Apr 2020
Original question from Google Cache:
"ODE solving using RK4 method (Predator prey)"
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
Diya Saha
Diya Saha 3 Apr 2020
@Stephen thankyou so much for doing that.
Rena Berman
Rena Berman 14 May 2020 16:10
(Answers Dev) Restored edit

로그인 to comment.

채택된 답변

James Tursa
James Tursa 3 Apr 2020
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( ).

  댓글 수: 1

Diya Saha
Diya Saha 3 Apr 2020
Thankyou so much! The of k1 and all the other k's being a 2-element vector really helped me edit my code and I got the answer. Thakyou again!

로그인 to comment.

추가 답변(0개)

이 질문에 답변하려면 로그인을(를) 수행하십시오.

태그

아직 태그를 입력하지 않았습니다.


Translated by