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

1 view (last 30 days)
Alvaro Sepulveda
Alvaro Sepulveda on 12 Dec 2021
Commented: Torsten on 12 Dec 2021
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.
  2 Comments
Alvaro Sepulveda
Alvaro Sepulveda on 12 Dec 2021
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

Sign in to comment.

Accepted Answer

Torsten
Torsten on 12 Dec 2021
Edited: Torsten on 12 Dec 2021
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 ?
  5 Comments
Torsten
Torsten on 12 Dec 2021
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.

Sign in to comment.

More Answers (1)

Jan
Jan on 12 Dec 2021
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.
  1 Comment
Alvaro Sepulveda
Alvaro Sepulveda on 12 Dec 2021
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

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!