Using identical delay and length of integration domain for delay model (dde23)

1 Ansicht (letzte 30 Tage)
Hi,
I've been trying to contsruct a biological process that is influenced by a drug every week, and causes few components of the immune system to react only 7-14 days later. Thus, the model consists a for loop with index of the number of weeks, and lags array for the delay. The results do differ from the model without delay, but if I change delay of 7 days to 14 days I get the same solution.
Is there any problem with my code? I just can't think of another option to model treatment without using for loop, but to me it seems that in this code there is something wrong with "retrieving" the history 7 or 14 days back. Thank you.
clc
a=6;
b=10;
c=9;
lags = [7 ];
y0 = [ 200 9.5*5*10^3 0.5*5*10^3 10^3 100 100];
tspan =[0,6/24];
n_weeks =6;
opt = odeset('AbsTol',1e-9,'RelTol',1e-6);
for i=1:3
sol =dde23(@(t,y,Z)ddefun(t,y,Z,a,b,c),lags,y0,tspan,opt);
semilogy(sol.x,sol.y(1,:),'r',sol.x,sol.y(2,:),'k',sol.x,sol.y(3,:),'b',sol.x,sol.y(4,:),'m',sol.x,sol.y(5,:),'g',sol.x,sol.y(6,:),'y')
legend('1','2','3','4','5','6')
xlabel('Time t');
ylabel('Solution y');
hold on
y0=sol.y(:,end);
tspan = [7*(i-1)+6/24,7*(i-1)+22/24];
m1=0;
sol =dde23(@(t,y,Z)ddefun(t,y,Z,a,b,c),lags,y0,tspan,opt);
semilogy(sol.x,sol.y(1,:),'r',sol.x,sol.y(2,:),'k',sol.x,sol.y(3,:),'b',sol.x,sol.y(4,:),'m',sol.x,sol.y(5,:),'g',sol.x,sol.y(6,:),'y')
legend('1','2','3','4','5','6')
xlabel('Time t');
ylabel('Solution y');
hold on
y0=sol.y(:,end);
mu1=0;
m1=0;
tspan = [7*(i-1)+22/24,7*i];
sol =dde23(@(t,y,Z)ddefun(t,y,Z,a,b,c),lags,y0,tspan,opt);
semilogy(sol.x,sol.y(1,:),'r',sol.x,sol.y(2,:),'k',sol.x,sol.y(3,:),'b',sol.x,sol.y(4,:),'m',sol.x,sol.y(5,:),'g',sol.x,sol.y(6,:),'y')
legend('1','2','3','4','5','6')
xlabel('Time t');
ylabel('Solution y');
hold on
y0=sol.y(:,end);
tspan = [7*i,7*i+6/24];
end
hold off
%-------------------------------------------
function dydt = ddefun(t,y,Z,a,b,c) %equation being solved
ylag1 = Z(:,1);
dydt = [ a*y(1);
b*ylag1(2)*y(1);
c*y(1)+a*ylag1(3);%
y(4);
ylag1(5)*ylag1(6);
y(6)];%
end

Antworten (1)

Karan Singh
Karan Singh am 30 Jan. 2024
Hi Marom,
The main issue with your code is that it does not properly account for the delay when updating "y0" and "tspan" within the loop. The delay in the biological response is not being modeled as intended, and this could be why changing the delay from 7 to 14 days does not affect the results.
When using "dde23" for delay differential equations, it's important to ensure that the history function (or the initial state y0) accounts for the delay period. However, in your loop, you're updating y0 with the end state of the previous iteration without considering the delay.
Here's a revised version of your loop that attempts to correct this by using a history function that can handle the delay:
clc
a = 6;
b = 10;
c = 9;
lags = [7]; % This should be in units consistent with tspan (e.g., days)
y0 = [200, 9.5*5*10^3, 0.5*5*10^3, 10^3, 100, 100];
tspan = [0, 6/24];
n_weeks = 6;
opt = odeset('AbsTol', 1e-9, 'RelTol', 1e-6);
% History function for the initial condition
history = @(t) y0;
for i = 1:n_weeks
% Define the time span for this iteration
tstart = 7 * (i - 1);
tend = 7 * i;
tspan = [tstart, tend];
% Solve the DDE
sol = dde23(@(t,y,Z) ddefun(t,y,Z,a,b,c), lags, history, tspan, opt);
% Plot the results
semilogy(sol.x, sol.y(1,:), 'r', sol.x, sol.y(2,:), 'k', sol.x, sol.y(3,:), 'b', sol.x, sol.y(4,:), 'm', sol.x, sol.y(5,:), 'g', sol.x, sol.y(6,:), 'y')
legend('1', '2', '3', '4', '5', '6')
xlabel('Time t');
ylabel('Solution y');
hold on
% Update the history function to include the latest solution
% This creates an anonymous function that interpolates the solution
% for times within the last week (or whatever the delay period is)
history = @(t) deval(sol, t);
end
hold off
%-------------------------------------------
function dydt = ddefun(t, y, Z, a, b, c) % equation being solved
ylag1 = Z(:,1);
dydt = [a*y(1);
b*ylag1(2)*y(1);
c*y(1)+a*ylag1(3);
y(4);
ylag1(5)*ylag1(6);
y(6)];
end
Hope this resolves your doubt!

Kategorien

Mehr zu Programming finden Sie in Help Center und File Exchange

Produkte


Version

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by