MATLAB: Runge-kutta 4th order Method Help

MATLABprogramming

% if use ode23 work
% but usin rk4 dont, please help!
% I use matlab to test and implement routines in fortran 90
% Tks!
h = 1; % step day
ti = 1; % first day
tf = 365; % one year
t = ti:h:tf;
y = zeros(1, length(t));
R1(1) = 1.0; % initial resource 1
R2(1) = 0.3; % initial resource 2
C1(1) = 1.0; % initial consumer 1
C2(1) = 1.0; % initial consumer 2
P(1) = 1.0; % initial predator
% prey-predator model
% P
% / \
% C1 C2
% | |
% R1 R2
% Paremeters
p = 0.50;
xc = 0.15;
xp = 0.08;
yc = 2.30;
yp = 1.70;
r0 = 0.25;
c0 = 0.50;
%% Equations
% R1
F_R1 = (R1.*(1-R1))-(xc*yc*((C1.*R1)./(R1+r0)));
% R2
F_R2 = (R2.*(1-R2))-(xc*yc*((C2.*R2)./(R2+r0)));
% C1
F_C1 = (xc*yc*((C1.*R1)./(R1+r0)))-((p.*C1./(p.*C1+(1-p).*C2)).*xp.*yp.*((P.*C1)/(c0+C1)))-xc.*C1;
% C2
F_C2 = (xc*yc*((C2.*R2)./(R2+r0)))-((((1-p).*C2)./(p.*C1+(1-p).*C2)).*xp.*yp.*((P.*C2)./(c0+C2)))-xc.*C2;
% P
F_P = ((p*C1./(p*C1+(1-p)*C2)).*xp.*yp.*((P.*C1)/(c0+C1)))+((((1-p)*C2)./(p*C1+(1-p)*C2)).*xp.*yp.*((P.*C2)./(c0+C2)))-xp.*P;
for i=1:(length(t)-1)
t(i+1)=t(i)+h;
kR1_1 = F_R1(t(i),R1(i));
kR1_2 = F_R1(t(i)+0.5*h,R1(i)+0.5*h*kR1_1);
kR1_3 = F_R1((t(i)+0.5*h),(R1(i)+0.5*h*kR1_2));
kR1_4 = F_R1((t(i)+h),(R1(i)+kR1_3*h));
R1(i+1) = R1(i) + (1/6)*(kR1_1+2*kR1_2+2*kR1_3+kR1_4)*h; % main equation
kR2_1 = F_R2(t(i),R2(i));
kR2_2 = F_R2(t(i)+0.5*h,R2(i)+0.5*h*kR2_1);
kR2_3 = F_R2((t(i)+0.5*h),(R2(i)+0.5*h*kR2_2));
kR2_4 = F_R2((t(i)+h),(R2(i)+kR2_3*h));
R2(i+1) = R2(i) + (1/6)*(kR2_1+2*kR2_2+2*kR1_3+kR2_4)*h;
kC1_1 = F_C1(t(i),C1(i));
kC1_2 = F_C1(t(i)+0.5*h,C1(i)+0.5*h*kC1_1);
kC1_3 = F_C1((t(i)+0.5*h),(C1(i)+0.5*h*kC1_2));
kC1_4 = F_C1((t(i)+h),(C1(i)+kC1_3*h));
C1(i+1) = C1(i) + (1/6)*(kC1_1+2*kC1_2+2*kC1_3+kC1_4)*h;
kC2_1 = F_C2(t(i),C2(i));
kC2_2 = F_C2(t(i)+0.5*h,C2(i)+0.5*h*kC2_1);
kC2_3 = F_C2((t(i)+0.5*h),(C2(i)+0.5*h*kC2_2));
kC2_4 = F_C2((t(i)+h),(C2(i)+kC2_3*h));
C2(i+1) = C2(i) + (1/6)*(kC2_1+2*kC2_2+2*kC2_3+kC2_4)*h;
kP_1 = F_P(t(i),P(i));
kP_2 = F_P(t(i)+0.5*h,P(i)+0.5*h*kP_1);
kP_3 = F_P((t(i)+0.5*h),(P(i)+0.5*h*kP_2));
kP_4 = F_P((t(i)+h),(P(i)+kP_3*h));
P(i+1) = P(i) + (1/6)*(kP_1+2*kP_2+2*kP_3+kP_4)*h;
end

Best Answer

All of your F's need to be functions that operate on the current full states. E.g., take this one
F_C1 = (xc*yc*((C1.*R1)./(R1+r0)))-((p.*C1./(p.*C1+(1-p).*C2)).*xp.*yp.*((P.*C1)/(c0+C1)))-xc.*C1;
It is obvious that the C1 derivative depends on the current states of C1, C2, R1, and P. So you need a function that takes these as inputs to compute the derivative. E.g., something like this
F_C1 = @(t,C1,C2,R1,P)(xc*yc*((C1.*R1)./(R1+r0)))-((p.*C1./(p.*C1+(1-p).*C2)).*xp.*yp.*((P.*C1)/(c0+C1)))-xc.*C1;
So modify your code to turn all of those F's into function handles as shown above.
Then you need to realize that all of your variables are interdependent, so using predictor methods will need to predict all of the states simultanously, not one at a time as you are coding. E.g., take these lines in your current code for C1:
    kC1_1 = F_C1(t(i),C1(i));
    kC1_2 = F_C1(t(i)+0.5*h,C1(i)+0.5*h*kC1_1);
    kC1_3 = F_C1((t(i)+0.5*h),(C1(i)+0.5*h*kC1_2));
    kC1_4 = F_C1((t(i)+h),(C1(i)+kC1_3*h));
The first line would change to this using our function handle from above:
    kC1_1 = F_C1( t(i), C1(i), C2(i), R1(i), P(i) );
But what about the next line? It would have to change to something like this:
    kC1_2 = F_C1( t(i)+0.5*h, C1(i)+0.5*h*kC1_1, C2(i)+0.5*h*kC2_1, R1(i)+0.5*h*kR1_1, P(i)+0.5*h*kP_1 );
where the C1, C2, R1, and P predictions are using their respective kC1_1 and kC2_1 and kR1_1 and kP_1 values. But some of these values haven't been computed yet for the current states, such as kP_1. Similar problems will exist in your other derivative calculations. Soooo ... you need to rearrange your code to compute all of the k_1 values first, then compute all of the k_2 values next, etc. E.g.,
kR1_1 = etc.
kR2_1 = etc.
kC1_1 = etc.
kC2_1 = etc.
kP_1 = etc.
kR1_2 = etc.
kR2_2 = etc.
kC1_2 = etc.
kC2_2 = etc.
kP_2 = etc.
   :
 etc.
Give a shot at making these changes and then come back to us with any further problems you might have.
Finally, at some point you will probably get tired of making lots of copy & paste of the RK4 code for each variable. One thing to look into is making your state a 5-element vector instead of using five different variables. That way you only need to write one set of RK4 code that will apply to this vector. This also makes it so that you can use the built-in functions such as ode45( ) to compare with.
Related Question