Using identical delay and length of integration domain for delay model (dde23)
조회 수: 1 (최근 30일)
이전 댓글 표시
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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/928269/image.png)
댓글 수: 0
답변 (1개)
Karan Singh
2024년 1월 30일
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!
댓글 수: 0
참고 항목
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!