Runge Kutta 4 - Biological Application

Hi,

My model of decomposition of debris (D) and loss of carbon (C1 and C2) and nitrogen (N) for the new microorganism biomass (B) works.

However, I would like some help to understand if I did it right.

Thank`s.

%%======================================================
clear all;
close all;
h = 1;
ti = 1;
tf = 1000;
t = ti:h:tf;

%% PARAMETERS

%% PARAMETERS HANNA PARNAS (1975)
% Gmax = maximal growth rate ofa mixture of decomposers 0.01 day-1
    Gmax = 0.01;

% kc = constant equivalent in its meaning to the Michaclis constant. but
% has area concentration units instead of volumetric concentration units. 
% Means the concentration of carbon which gives half the maximal 
% growth rate 20 g * m-2 or 0.020 kg * m-2
    kc = 0.020;

% kn = means the concentration of nitrogen which gives the half maximal 
% growth rate 0.4 g * m-2 or 0.0004 kg * m-2
    kn = 0.0004;

% F = ratio of the carbon assimilated into decomposers to the total amount 
% of carbon released from detritus by decomposers
% F = efficiency of carbon assimilation = 0.4 
% FUTURE CUE
    F = 0.40;

% fc = average fraction of carbon in the decomposer's cells
% fc = fraction of carbon in the microbial biomass = 0.5
    fc = 0.50;

% fn = average fraction of nitrogen in the decomposer's cells.
% fn = fraction of nitrogen In the microbial biomass = 0.05
    fn = 0.05;

% alfa = is the total carbon used by the decomposer organisms per unit 
% decomposer biomass increment.
% alfa = fc./F;
    alfa = 1.25;





%% Initial Values


    N1D(1) = 6.25;

    NH4(1) = 0.01;

    ND(1) = 1.00;
 
    C1D(1) = 25.00;

    C2D(1) = 0.00;

    CD(1) = C1D(1) + C2D(1);

    B(1) = 0.10;

    CTB(1) = fc.*B(1);

    NTB(1) = fn.*B(1);
           
    G(1) = (Gmax.*CD(1).*ND(1))./((kc+CD(1)).*(kn+ND(1)));      
 
    in(1) = 0.00;
 
    beta(1) = 4.00;



    %F_N1D = @(G,B) - (fn.*G.*B - in);

    %F_C1D = @(G,B) - beta.*(fn.*G.*B - in);

    %F_C2D = @(G,B) - (alfa.*G.*B - beta.*(fn.*G.*B - in));

    %F_CTB = @(G,B) F.*alfa.*G.*B;

    %F_NTB = @(G,B) fn.*G.*B;



    F_N1D = @(t,N1D) - (fn.*G.*B - in);

    F_C1D = @(t,C1D) - beta.*(fn.*G.*B - in);

    F_C2D = @(t,C2D) - (alfa.*G.*B - beta.*(fn.*G.*B - in));

    F_CTB = @(t,CTB) F.*alfa.*G.*B;

    F_NTB = @(t,NTB) fn.*G.*B;



    N1D     = zeros(1,length(t));
    NH4     = zeros(1,length(t));
    ND      = zeros(1,length(t));
    C1D     = zeros(1,length(t));
    C2D     = zeros(1,length(t));
    CD      = zeros(1,length(t));
    B       = zeros(1,length(t));
    CTB     = zeros(1,length(t));
    NTB     = zeros(1,length(t));
    G       = zeros(1,length(t));
    in      = zeros(1,length(t));
    beta    = zeros(1,length(t));




  for i=1:(length(t)-1) 
       t(i+1)=t(i)+h;  
  
    kN1D_1 = F_N1D( t(i), N1D(i) );
    kC1D_1 = F_C1D( t(i), C1D(i) );
    kC2D_1 = F_C2D( t(i), C2D(i) );
    kCTB_1 = F_CTB( t(i), CTB(i) );
    kNTB_1 = F_NTB( t(i), NTB(i) );   
    
    
    kN1D_2 = F_N1D( t(i)+0.5*h, N1D(i)+0.5*h*kN1D_1 );
    kC1D_2 = F_C1D( t(i)+0.5*h, C1D(i)+0.5*h*kC1D_1 );
    kC2D_2 = F_C2D( t(i)+0.5*h, C2D(i)+0.5*h*kC2D_1 );
    kCTB_2 = F_CTB( t(i)+0.5*h, CTB(i)+0.5*h*kCTB_1 );
    kNTB_2 = F_NTB( t(i)+0.5*h, NTB(i)+0.5*h*kNTB_1 );
    
    kN1D_3 = F_N1D( t(i)+0.5*h, N1D(i)+0.5*h*kN1D_2 );
    kC1D_3 = F_C1D( t(i)+0.5*h, C1D(i)+0.5*h*kC1D_2 );
    kC2D_3 = F_C2D( t(i)+0.5*h, C2D(i)+0.5*h*kC2D_2 );
    kCTB_3 = F_CTB( t(i)+0.5*h, CTB(i)+0.5*h*kCTB_2 );
    kNTB_3 = F_NTB( t(i)+0.5*h, NTB(i)+0.5*h*kNTB_2 );
    
    kN1D_4 = F_N1D( t(i)+h, N1D(i)+kN1D_3*h );
    kC1D_4 = F_C1D( t(i)+h, C1D(i)+kC1D_3*h );
    kC2D_4 = F_C2D( t(i)+h, C2D(i)+kC2D_3*h );
    kCTB_4 = F_CTB( t(i)+h, CTB(i)+kCTB_3*h );
    kNTB_4 = F_NTB( t(i)+h, NTB(i)+kNTB_3*h );
    
    
    N1D(i+1) = N1D(i) + (1/6)*(kN1D_1+2*kN1D_2+2*kN1D_3+kN1D_4)*h;
    C1D(i+1) = C1D(i) + (1/6)*(kC1D_1+2*kC1D_2+2*kC1D_3+kC1D_4)*h;
    C2D(i+1) = C2D(i) + (1/6)*(kC2D_1+2*kC2D_2+2*kC2D_3+kC2D_4)*h;
    CTB(i+1) = CTB(i) + (1/6)*(kCTB_1+2*kCTB_2+2*kCTB_3+kCTB_4)*h;   
    NTB(i+1) = NTB(i) + (1/6)*(kNTB_1+2*kNTB_2+2*kNTB_3+kNTB_4)*h;
     
 
    
    NH4(i) = 0.01;
     
    ND(i) = N1D(i) + NH4(i);
 
    CD(i) = C1D(i) + C2D(i);
  
    B(i) = (CTB(i)./fc);
           
    G(i) = (Gmax.*CD(i).*ND(i))./((kc+CD(i)).*(kn+ND(i)));      
 
    in(i) = (fn.*G(i).*B(i).*NH4(i))./(N1D(i) + NH4(i));     
    
    beta(i) = C1D(i)./N1D(i);
 
 
 
  end
       

Answer 1

[Jan]

I would like some help to understand if I did it right.

I cannot estimate what you consider as right or not right.

For a real scientific appraoch you would never wirte the integrator by your own and mix it with the actaul application, because you do not have a chance to test the integrator with some test functions. Using a well tested function as ODE45 would be much smarter, especially due to the step size control. If your ODE is stiff, the fixed stepwidth does not give you any warning about this.

clear all; deletes all loaded functions from the memory, such that running any code needs to relaod the files from the slow disk. There is not advantage in doing this, but this is a waste of time only. Prefer to store the code in function, then there is no reason for clearing.

This line

t = ti:h:tf;

defines t already. Then you can omit:

t(i+1)=t(i)+h;

It is not wrong, but only a little bit confusing to use .* and ./ for scalar values. The standard operators * and / are sufficient.