ODE solving using RK4 method (Predator prey)
조회 수: 1 (최근 30일)
이전 댓글 표시
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
채택된 답변
James Tursa
2020년 4월 3일
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( ).
추가 답변 (0개)
참고 항목
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!