Runge Kutta 4 - Oriented Object Alternative to Biological Aplication

Hello everyone!

I am trying to develop an alternative version of the Runge-Kutta 4 object-oriented for biological applications.

However, the values are not correct when comparing to the standard system.

Alternative system (TED_MICRO - attached):

%%
%==========================================================================
%
%  NECROMASS MOISTURE DYNAMIC MODEL - NMDM 
%
%                                SETUP
%
%==========================================================================
%==========================================================================

    clear all;
    close all;

    
%% PREPARATION TO INTEGRATION LOOP    

%PARAMETERS

    a = 0.5471;
    b = 0.0281;
    c = 0.8439;
    d = 0.0266;
    
% INITIAL CONDICTIONS

    PY(1) = 30.00;
    PD(1) =  4.00;
    
    
    cwd = cwd_init_state(PD,PY);  % CWD = COARSE WOODY DEBRIS
    
        
% =========================================================================       
    
    

%% START TIME INTEGRATION LOOP TO ALL SUBMODELS

    h = 0.001;
    ti = 1;       
    tf = 100;    % TOTAL DAYS
    t = ti:h:tf;
  
    
    for it=1:(length(t)-1) 
    
   
           
        cwd_rk = cwd_initial_rk(cwd,it); % CWD = COARSE WOODY DEBRIS
        
        
    for isub=1:4
%%
% =========================================================================
% =========================================================================
%                           
%                  RUNGE-KUTTA 4 ORDER INTEGRATION LOOP
%
% =========================================================================
% =========================================================================  
        
        
        % The purpose here is to update the state to the appropriate
        % partial-step position
        % k1 is calcualated at it
        % k2 is calculated at dt_sec/2, k1/2
        % k3 is calculated at dt_sec/2, k2/2
        % k4 is calculated at dt_sec, k3
        
        % This subroutine updates sub-step diagnostic variables, assuming
        % That qmoist is known
        
    
        
        cwd_rk = cwd_updates_rk_new(cwd,cwd_rk,it,isub);
        
        cwd_rk = cwd_derivatives_rk(cwd_rk,h,a,b,c,d);
       
        
    end
 
    
%%    
% =========================================================================
% =========================================================================
%                           
%                     INTEGRATION PROCESS
%
% =========================================================================
% =========================================================================    
    
    
        cwd = cwd_integrate_rk(cwd_rk,cwd,it);
    
 
  
    end
    
    ax1 = subplot(2,1,1); % top subplot
    x = linspace(0,100);
    plot(t,cwd.PY,t,cwd.PD)
    title(ax1,'Micro Dynamic:ALTERNATIVE')
    ylabel(ax1,'Population Size')
    ax2 = subplot(2,1,2); % bottom subplot
    plot(cwd.PY,cwd.PD)
    title(ax2,'Micro Dependence:ALTERNATIVE')
    ylabel(ax2,'Population Size')
    

Standard system:

    clear all;
close all;
h = 0.01;
ti = 1;
tf = 100;
t = ti:h:tf;

PY = zeros(1,length(t));
PD = zeros(1,length(t));

cwd.PY = zeros(1,length(t));
cwd.PD = zeros(1,length(t));


PY(1) = 30.00;
PD(1) =  4.00;

a = 0.5471;
b = 0.0281;
c = 0.8439;
d = 0.0266;


F_PY = @(t,PY,PD)  a*PY - b*PY*PD;

F_PD = @(t,PY,PD) -c*PD + d*PY*PD;


for i = 1:length(t)-1
    
   KPY_1 = F_PY(t(i),PY(i),PD(i)); 
   KPD_1 = F_PD(t(i),PY(i),PD(i)); 

   KPY_2 = F_PY(t(i)+h*0.5,PY(i)+h*0.5*KPY_1,PD(i)+h*0.5*KPD_1); 
   KPD_2 = F_PD(t(i)+h*0.5,PY(i)+h*0.5*KPY_1,PD(i)+h*0.5*KPD_1); 

   KPY_3 = F_PY(t(i)+h*0.5,PY(i)+h*0.5*KPY_2,PD(i)+h*0.5*KPD_2); 
   KPD_3 = F_PD(t(i)+h*0.5,PY(i)+h*0.5*KPY_2,PD(i)+h*0.5*KPD_2);
   
   KPY_4 = F_PY(t(i)+h,PY(i)+h*KPY_3,PD(i)+h*KPD_3); 
   KPD_4 = F_PD(t(i)+h,PY(i)+h*KPY_3,PD(i)+h*KPD_3);
   
   PY(i+1) = PY(i) + (h/6)*(KPY_1+2*KPY_2+2*KPY_3+KPY_4);
   PD(i+1) = PD(i) + (h/6)*(KPD_1+2*KPD_2+2*KPD_3+KPD_4);
   
   cwd.PY(i) = PY(i+1);
   cwd.PD(i) = PD(i+1);
   
   
end

ax1 = subplot(2,1,1); % top subplot
x = linspace(0,100);
plot(t,PY,t,PD)
title(ax1,'Micro Dynamic: STANDART')
ylabel(ax1,'Population Size')
ax2 = subplot(2,1,2); % bottom subplot
plot(PY,PD)
title(ax2,'Micro Dependence: STANDART')
ylabel(ax2,'Population Size')

Can you help me?

I really need to develop object-oriented Runge-Kutta 4.

[darova]

Did you try to decrease timestep?

[Paulo Araujo Filho]

Hi, yes.

The result is the same.