ODE solving using RK4 method (Predator prey)
1 visualización (últimos 30 días)
Mostrar comentarios más antiguos
Diya Saha
el 3 de Abr. de 2020
Comentada: Rena Berman
el 14 de Mayo de 2020
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 comentarios
Respuesta aceptada
James Tursa
el 3 de Abr. de 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 comentario
Más respuestas (0)
Ver también
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!