change the value of a parameter Lambda at each time step

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);

Answer 1

[Torsten]

% 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 

[ahmadou]

Thank you so much @Torsten. Youve been a great help.