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

2021 12 12
2 个回答
回答已采纳
更新时间：2021 12 12
13 次查看（30 天）

0 个投票

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 个评论
显示 无 隐藏 无

Torsten
Torsten 2021-12-12
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
Alvaro Sepulveda 2021-12-12
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

请先登录，再进行评论。

请先登录，再回答此问题。

 采纳的回答

Torsten
Torsten 2021-12-12
编辑：Torsten 2021-12-12

0 个投票

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 个评论
显示 3更早的评论 隐藏 3更早的评论

Alvaro Sepulveda
Alvaro Sepulveda 2021-12-12
The values of the parameters are given in weeks, so the at the end of the simulation is equivalent to 2000 weeks
Torsten
Torsten 2021-12-12
编辑：Torsten 2021-12-12
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
Alvaro Sepulveda 2021-12-12
Yes, it is correct, the solutions of the system will be given in weeks
Alvaro Sepulveda
Alvaro Sepulveda 2021-12-12
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
Torsten 2021-12-12
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.

请先登录，再进行评论。

更多回答（1 个）

Jan
Jan 2021-12-12

0 个投票

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 个评论
显示 -1更早的评论 隐藏 -1更早的评论

Alvaro Sepulveda
Alvaro Sepulveda 2021-12-12
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

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Linear Algebra 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by