change the value of a parameter Lambda at each time step

1 vue (au cours des 30 derniers jours)
ahmadou
ahmadou le 3 Avr 2024
Modifié(e) : ahmadou le 3 Avr 2024
Hello everyone,
I have a code that works well when the parameters are constant. Now I'm trying to make the Lambda parameter time-dependent, but I can't do it. Can someone please help me?
function Z = fun_model(p,t,mu,Lambda)
% Unknown parameters
beta1 = p(1);
beta2 = p(2);
sigma = p(3);
eps = p(4);
gamma = p(5);
eta = p(6);
alpha = p(7);
Delta = p(8);
theta = p(9);
y0 = p(10:17); % initial conditions
a = p(18);
b= p(19);
% systems
f = @(t,y) [ mu*(y(1)+y(2)+y(3)+y(4)+y(5)) - mu*y(1) - beta1*y(1)*y(7) - beta2*y(1)*y(8); ...
beta1*y(1)*y(7) + beta2*y(1)*y(8) - sigma*y(2) - mu*y(2); ...
sigma*y(2) - eps*y(3) - mu*y(3); ...
eps*y(3) - gamma*y(4) - mu*y(4); ...
gamma*y(4) - mu*y(5); ...
eta*(y(6)+y(7)) - alpha*y(6)*y(8) - eta*y(6); ...
alpha*y(6)*y(8) - eta*y(7); ...
(a+b*Lambda)*y(8) - Delta*y(8) - theta*eta*y(7) ];
% Resolution with ode45
options = odeset('RelTol',1e-6,'AbsTol',1e-6); % Adjust tolerances for performance
[~,ypred] = ode45(f,t,y0,options);
% Return the desired output
Z = [ypred(:,3), ypred(:,7), ypred(:,8)];
end
% Experiemental data on 12 months
time = 1:12 ; % experimental time 12 months
Yc = [20 18 19 24 28 26 25 24 20 19 21 19]' ; % Bu cases
W = [0.14 0.08 0.18 0.08 0 0.37 0.24 0.18 0 0.1 0.08 0.05]'; % MU prevalence in vectors
U = [0.05 0.02 0.05 0.07 0 0.09 0.1 0.16 0.05 0.13 0.07 0.06]'; %MU positivity
beta1 = 0.1; beta2 = 0.1; sigma = 0.2; eps = 0.16; gamma = 0.17;
eta = 0.7; alpha = 0.8; Delta = 0.7; theta = 0.9;
y0 = [6500; 0; 20; 0; 0; 0.63; 0.37; 1] ; a=1; b=1;
mu = 0.00216; % Mortality/fertility rate 2,6% by year
Lambda = 0.29;
p0 = [beta1; beta2; sigma; eps; gamma; eta; alpha; Delta; theta; y0; a; b] ; % Initial guess for parameters
%Parameter upper, lower, and Using lsqcurvefit to adjust model parameters
lb = [0, 0, 1/5, 1/6, 1/6, 0, 0, 0, 0, 0,0,0,0,0,0,0,0, -inf,-inf]; % lower bound
ub = [1, 1, 1/3, 1.0, 1/2, 1, 1, 1, 1, 10^5,50,50,50,50,2,2,3, inf,inf]; % upper bound
opts = optimoptions('lsqcurvefit');
opts.MaxFunEvals = 10000;
[p_fit,resnorm,residual,exitflag,output,jacobian] = lsqcurvefit(@(p,t) fun_model(p,t,mu,Lambda),p0,time,[Yc,W,U],lb,ub,opts);

Réponse acceptée

Torsten
Torsten le 3 Avr 2024
Modifié(e) : Torsten le 3 Avr 2024
% Experiemental data on 12 months
time = 1:12 ; % experimental time 12 months
Yc = [20 18 19 24 28 26 25 24 20 19 21 19]' ; % Bu cases
W = [0.14 0.08 0.18 0.08 0 0.37 0.24 0.18 0 0.1 0.08 0.05]'; % MU prevalence in vectors
U = [0.05 0.02 0.05 0.07 0 0.09 0.1 0.16 0.05 0.13 0.07 0.06]'; %MU positivity
beta1 = 0.1; beta2 = 0.1; sigma = 0.2; eps = 0.16; gamma = 0.17;
eta = 0.7; alpha = 0.8; Delta = 0.7; theta = 0.9;
y0 = [6500; 0; 20; 0; 0; 0.63; 0.37; 1] ; a=1; b=1;
mu = 0.00216; % Mortality/fertility rate 2,6% by year
Lambda = @(t) ??? % Your time-dependent function for Lambda
p0 = [beta1; beta2; sigma; eps; gamma; eta; alpha; Delta; theta; y0; a; b] ; % Initial guess for parameters
%Parameter upper, lower, and Using lsqcurvefit to adjust model parameters
lb = [0, 0, 1/5, 1/6, 1/6, 0, 0, 0, 0, 0,0,0,0,0,0,0,0, -inf,-inf]; % lower bound
ub = [1, 1, 1/3, 1.0, 1/2, 1, 1, 1, 1, 10^5,50,50,50,50,2,2,3, inf,inf]; % upper bound
opts = optimoptions('lsqcurvefit');
opts.MaxFunEvals = 10000;
[p_fit,resnorm,residual,exitflag,output,jacobian] = lsqcurvefit(@(p,t) fun_model(p,t,mu,Lambda),p0,time,[Yc,W,U],lb,ub,opts);
function Z = fun_model(p,t,mu,Lambda)
% Unknown parameters
beta1 = p(1);
beta2 = p(2);
sigma = p(3);
eps = p(4);
gamma = p(5);
eta = p(6);
alpha = p(7);
Delta = p(8);
theta = p(9);
y0 = p(10:17); % initial conditions
a = p(18);
b= p(19);
% systems
f = @(t,y) [ mu*(y(1)+y(2)+y(3)+y(4)+y(5)) - mu*y(1) - beta1*y(1)*y(7) - beta2*y(1)*y(8); ...
beta1*y(1)*y(7) + beta2*y(1)*y(8) - sigma*y(2) - mu*y(2); ...
sigma*y(2) - eps*y(3) - mu*y(3); ...
eps*y(3) - gamma*y(4) - mu*y(4); ...
gamma*y(4) - mu*y(5); ...
eta*(y(6)+y(7)) - alpha*y(6)*y(8) - eta*y(6); ...
alpha*y(6)*y(8) - eta*y(7); ...
(a+b*Lambda(t))*y(8) - Delta*y(8) - theta*eta*y(7) ];
% Resolution with ode45
options = odeset('RelTol',1e-6,'AbsTol',1e-6); % Adjust tolerances for performance
[~,ypred] = ode45(f,t,y0,options);
% Return the desired output
Z = [ypred(:,3), ypred(:,7), ypred(:,8)];
end
  1 commentaire
ahmadou
ahmadou le 3 Avr 2024
Modifié(e) : ahmadou le 3 Avr 2024
Thank you so much @Torsten. Youve been a great help.

Connectez-vous pour commenter.

Plus de réponses (0)

Catégories

En savoir plus sur Audio I/O and Waveform Generation dans Help Center et File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!

Translated by