Really am I simulation this system with ODE45? How am I can optimizing the code?

Hi, everyone,

I new in matlab and I wrote this script that simulates a system of odes equations that interacting whit 3 parameters that varying on time:

.

The ODEs equations are:

form how will varaying the parameters is:

where , is pluv in the code and is temp in the code

this is the code that am I wrote:

%==========================================
tic;
% Random seed at the start of the simulation dependent on the computer clock 
rand('state',sum(100*clock)); %

numreps=1; % Realization numbers
%for j = 100:numreps % Repetitions numbe
    options = odeset('RelTol',1e-6,'Stats','on'); 
    
    %==================================
    %**********************************
    % ---- parameters ----
    %**********************************
    %==================================
    %
    % 
    ah = 0.36; %
    %
    % 
    bH = 0.085; % 
    bh = 0.85; %
    %
    % 
    Kbh = 9.7e-04;%
    %
    % 
    Kh = 9.7e-04;%

    %=====================================
    %*************************************
    % ---- parameters ----
    %*************************************
    %=====================================
    %
    % 
    alphaH = 10.^-5; % 
    alphah = 10.^-5; %
    %
    % 
    lambdaH = 8 * 10.^-1; %
    lambdah = 8 * 10.^-1; %
    %
    % 
    muH = 7 * 10.^-3; % 
    muh = 7 * 10.^-3; % 
    %
    muZ = 0.80; % 
    %
    %  
    betaH = 0.0125;  % 
    betah = 0.0125;  % 
    %
    phiH = 0.0017; % 
    phih = 0.0017; % 
    %
    %===============================================
    %***********************************************
    % ---- Values of temperature and pluviosity ----
    %***********************************************
    %===============================================
    % 
    tmax=27;    tmin=4;    plumax=200;     plumin=10;
    %
    S = 52; %
    %
    % Simulation time for rH, rh and ah and aH
    %tt = (0:0.1:300);  % Generate t for aH, ah, rH, rh
    %
    tspan = (0:2000);
    tt = tspan;
    %=======================================
    % Values minimum and maximum of aH
    %=======================================
    %
    aHmin = 45; % 
    %
    aHmax = 70; %  
    %
    %===================================================
    % Maximum / minimum rh and rH 
    %====================================================
    %
    rHmin = 0.01; % 
    %
    rHmax = 0.2;  % 
    %
    rhmin = 0.01; % 
    %
    rhmax = 0.2;  
    %
    %==================================
    % ---- Function of pluviosity  ----
    %==================================
    pluv = (1./2)*(2*plumax-(plumax-plumin)*(1-cos((2*pi./52)*(tt-1./2)+pi))); %
    %=================================
    % ---- Function of Temperature ----
    %==================================
    %
    temp = (1./2)*(2*tmax-(tmax-tmin)*(1-cos((2*pi./52)*(tt-1/2)+pi))); % 
    %
    te=temp; pl=pluv;
    %
    %=======================================================
    %----- Growing function due temperature ----
    %=======================================================
    rh=rhmax-((rhmax-rhmin)/(tmax-tmin))*(te-tmin); 
    %
    %=======================================================
    %----- Growing function due temperature ---- 
    %=======================================================
    rH=rHmax-((rHmax-rHmin)/(tmax-tmin))*(te-tmin);
    %
    %=================================================================
    %----- Fecundity due pluviosity----
    %=================================================================
    aH = aHmax-((aHmax-aHmin)/(S-plumin))*(pl-plumin); 
    if aH < 0
        aH = 0;
    end
    aH = round(aH);
    %================================================================================
    %-------------------- Numeric solution ode45 --------------
    %================================================================================
    x0 = [50 50 100 100 500]; % intitials conditions
    [t,xa] = ode45(@(t,x) bd_host_ode(t, x, tt, lambdaH, lambdah,muZ, aH, bh, bH, Kbh, alphaH,alphah,...
        betaH, betah, rH, rh, muH, muh,ah, Kh, phiH, phih),tspan, x0, options);
  
     %================================================================================= 
     %
     %================================================      
     %----- Solutions of model -----
     %================================================
     % 
     H = xa(:,1); % 
     h = xa(:,2); % 
     %
     % 
     ZH = xa(:,3); % 
     Zh = xa(:,4); % 
     Z = xa(:,5); % 
     %
 
    % Function that describe of ODES model
    function dxdt = bd_host_ode(t, x, tt, lambdaH, lambdah,muZ, aH, bh, bH, Kbh, alphaH,alphah,...
            betaH, betah, rH, rh, muH, muh,ah, Kh, phiH, phih)
        rh = interp1(tt,rh,t);
        rH = interp1(tt,rH,t);
        aH = interp1(tt,aH,t);
        dxdt = zeros(5,1);
        dxdt(1) = ah * exp(-Kh * x(2)) * x(2) - bH * x(1) - alphaH * x(3);
        dxdt(2) = aH * x(1) - (ah * exp(-Kh * x(2)) + bh - bh * exp(-Kbh * x(2))) * x(2) - alphah * x(4);
        dxdt(3) = x(5) * betaH * (x(1)./(x(2)+x(1))) + x(3) * (rH - bH - muH - alphaH * (1 + (x(3) * (phiH + 1)./x(1) * phiH)));
        dxdt(4) = x(5) * betah * (x(2)./(x(2)+x(1))) + x(4) * (rh - ah * exp(-Kh * x(2)) - bh + bh * exp(-Kbh * x(2)) - muh - alphah * (1 + (x(4) * (phih + 1)./x(2) * phih)));
        dxdt(5) = lambdaH * x(3) + lambdah * x(4) - muZ * x(5) - (x(5) * betaH * (x(1)./(x(2)+x(1))) + x(5) * betah * (x(2)./(x(2)+x(1))));
    end
toc;

My question is: Really am I simulation this system with ODE45? How am I can optimizing the code?

I'm not sure if the code is it fine.

I hope someone kind can check my code to see if it is correctly written.

[Torsten]

Why don't you compute rh,rH and aH from t directly in function bd_host_ode ?

This will give smooth behaviour of parameters and avoid interpolation.

[Alvaro Sepulveda]

Hi @Torsten,

Why don't you compute rh,rH and aH from t directly in function bd_host_ode ?

you would be so kind to tell me how to do this. Thank you very much for the help you can give me

Answer 1

[Torsten]

function main
options = odeset('RelTol',1e-6,'AbsTol',1.0e-6,'Stats','on'); 
tspan = (0:0.1:8.9);

%================================================================================
%-------------------- Numeric solution ode45 --------------
%================================================================================
x0 = [50 50 100 100 500]; % intitials conditions
[t,xa] = ode15s(@bd_host_ode,tspan, x0, options);

%================================================================================= 
%
%================================================      
%----- Solutions of model -----
%================================================
% 
H = xa(:,1); % 
h = xa(:,2); % 
%
% 
ZH = xa(:,3); % 
Zh = xa(:,4); % 
Z = xa(:,5); % 
plot(t,Z)%,t,h,t,ZH,t,Zh,t,Z)
%
end
% Function that describe of ODES model
function dxdt = bd_host_ode(t, x)

%**********************************
% ---- parameters ----
%**********************************
%==================================
%
% 
ah = 0.36; %
%
% 
bH = 0.085; % 
bh = 0.85; %
%
% 
Kbh = 9.7e-04;%
%
% 
Kh = 9.7e-04;%
%=====================================
%*************************************
% ---- parameters ----
%*************************************
%=====================================
%
% 
alphaH = 10.^-5; % 
alphah = 10.^-5; %
%
% 
lambdaH = 8 * 10.^-1; %
lambdah = 8 * 10.^-1; %
%
% 
muH = 7 * 10.^-3; % 
muh = 7 * 10.^-3; % 
%
muZ = 0.80; % 
%
%  
betaH = 0.0125;  % 
betah = 0.0125;  % 
%
phiH = 0.0017; % 
phih = 0.0017; % 
%
%===============================================
%***********************************************
% ---- Values of temperature and pluviosity ----
%***********************************************
%===============================================
% 
tmax=27;    tmin=4;    plumax=200;     plumin=10;
%
S = 52; %    
aHmin = 45; % 
%
aHmax = 70; %  
%
%===================================================
% Maximum / minimum rh and rH 
%====================================================
%
rHmin = 0.01; % 
%
rHmax = 0.2;  % 
%
rhmin = 0.01; % 
%
rhmax = 0.2;  
%
%==================================
% ---- Function of pluviosity  ----
%==================================
pluv = (1./2)*(2*plumax-(plumax-plumin)*(1-cos((2*pi./52)*(t-1./2)+pi))); %
%=================================
% ---- Function of Temperature ----
%==================================
%
temp = (1./2)*(2*tmax-(tmax-tmin)*(1-cos((2*pi./52)*(t-1/2)+pi))); % 
%
te=temp; pl=pluv;
%
%=======================================================
%----- Growing function due temperature ----
%=======================================================
rh=rhmax-((rhmax-rhmin)/(tmax-tmin))*(te-tmin); 
%
%=======================================================
%----- Growing function due temperature ---- 
%=======================================================
rH=rHmax-((rHmax-rHmin)/(tmax-tmin))*(te-tmin);
%
%=================================================================
%----- Fecundity due pluviosity----
%=================================================================
aH = aHmax-((aHmax-aHmin)/(S-plumin))*(pl-plumin); 
aH = max(aH,0);
dxdt = zeros(5,1);
dxdt(1) = ah * exp(-Kh * x(2)) * x(2) - bH * x(1) - alphaH * x(3);
dxdt(2) = aH * x(1) - (ah * exp(-Kh * x(2)) + bh - bh * exp(-Kbh * x(2))) * x(2) - alphah * x(4);
dxdt(3) = x(5) * betaH * x(1)./(x(2)+x(1)) + x(3) * (rH - bH - muH - alphaH * (1 + x(3)/x(1) * (phiH + 1)./ phiH));
dxdt(4) = x(5) * betah * x(2)./(x(2)+x(1)) + x(4) * rh - x(4)*(ah * exp(-Kh * x(2)) - bh + bh * exp(-Kbh * x(2)) - muh - alphah * (1 + x(4)/x(2) * (phih + 1)./ phih));
dxdt(5) = lambdaH * x(3) + lambdah * x(4) - x(5)*( muZ + betaH * x(1)./(x(2)+x(1)) + betah * x(2)./(x(2)+x(1)));
end

I corrected some discrepancies between your original code and the differential equations given in mathematical notation, too.

Further I doubt that t in the cos() expression must be given in seconds, but in weeks. Or is your end of the simulation equal to 2000 weeks ?

[Alvaro Sepulveda]

The values of the parameters are given in weeks, so the at the end of the simulation is equivalent to 2000 weeks

[Torsten]

Ok. Then I hope that units in the parameters involved are consistent with t in weeks.

Check aH, rh and rH for correctness. They also don't seem to be consistent with the formulae in mathematical notation.

[Alvaro Sepulveda]

Yes, it is correct, the solutions of the system will be given in weeks

[Alvaro Sepulveda]

Thank you very much for taking the time and reviewing my code, as well as helping me with the suggestions. And you are right about the rk and aH functions that need to be corrected.

One last question,

One last question, why do you use ode15s instead of ode45?

[Torsten]

ode15s can handle non-stiff and stiff, ode45 only non-stiff ordinary differential equations. Because I don't know about the stiffness of your system, I wanted to save an unnecessary attempt to solve your equations with ode45 if they really were stiff.

But after integration, all solvers should come up with the same result - just check which one is the most efficient for your application.

Answer 2

[Jan]

Your function to integrated contains linear intepolations of parameters. Then the outpuzt is not smooth. Matlab's ODE intergrators are deigned fpr smooth functions only. For other functions the stepsize controller can drive crazy: Either the integration stops with an error message or with less luck it reduces the stepsize dramatically, such that the accumulated rounding errors can dominate the resulting trajectory. Then the integrator is an expensive pseudo-random-number generator with a poor entropy.

[Alvaro Sepulveda]

and how can i fix this? should I write the code again? I have found it difficult trying to write this code with parameters that vary over time. I have seen very few examples that help me as a guide to write my codeI have seen very few examples that help me as a guide to write my code