How to solve a system of distributed delay equations?

I have a code, which gives a solution of a system of discrete delay equations.

This is how I run it

lags=1;
tspan=[0 600];
sol=ddesd(@ddefunc,lags,[0.2; 0.08],tspan);
p=plot(sol.x,sol.y);
set(p,{'LineWidth'},{2;2})
title('y(t)')
xlabel('Time(days)'), ylabel('populations')
legend('x','y')

and this is the function

function yp = ddefunc(~,y,Z)

a=0.1;
b=0.05;
c=0.08;
d=0.02;

yl1=Z(:,1);

yp = [a*y(1)-b*y(1)*yl1(2);
    c*y(1)*y(2)-d*y(2)];
end

Now, instead of one discrete delay value, I would like to consider a continuous delay values. That is, instead of , . Would it be possible to do this? Thanks!

Answer 1

[Torsten]

As a start, you could define three delays, namely delay(1) = tau-gamma, delay(2) = tau and delay(3) = tau+gamma, and approximate the integral as "gamma * ( Z(:,1)/2 + Z(:,2) + Z(:,3)/2 )" (trapezoidal rule with three points).

In principle, you can approximate the integral arbitrarily close by choosing a sufficient number of delays:

tau = 1;
gamma = 0.5;
number_of_delays = 11;   % should be odd
lags = linspace(tau-gamma,tau+gamma,number_of_delays);
tspan=[0 600];
sol=ddesd(@(t,y,Z)ddefunc(t,y,Z,lags),lags,[0.2; 0.08],tspan);
p=plot(sol.x,sol.y);
set(p,{'LineWidth'},{2;2})
title('y(t)')
xlabel('Time(days)'), ylabel('populations')
legend('x','y')

function yp = ddefunc(~,y,Z,lags)
  a=0.1;
  b=0.05;
  c=0.08;
  d=0.02;

  yl1 = trapz(lags,Z(2,:));

  yp = [a*y(1)-b*y(1)*yl1;
        c*y(1)*y(2)-d*y(2)];
end

[Fares]

Thank you Torsten!

[José]

Dear Torsten, thanks for the solution. I have two questions:

1) What if we integrate with negative values, e.g., if the distributed delay is: where . How do we adapt the code?

2) I'm trying to distribute the delay with a normalized truncated Gaussian-type kernel, specifically the distributed delay: where and . How do we adapt the code?

Thank you very much. Best regards.

[Torsten]

The interval you integrate over is [t-2*r,t]. Thus the lags should be defined as

lags = linspace(0,2*r,number_of_delays)

And for the second question, "trapz" should be applied as

yl1 = trapz(lags,Z(2,:).*exp(-((-lags+r).^2)/2));

to approximate

integral_{s=-2*r}^{s=0} exp(-((s+r)^2)/2)*u(t+s) ds = integral_{s'=0}^{s'=2*r}exp(-((-s'+r)^2)/2)*u(t-s') ds'

[Walter Roberson]

You cannot solve with negative delays. Negative delays are predictions about the future.

[Torsten]

@Walter Roberson

The delays @José uses are non-negative. Note that his notation is u(t+s) instead of u(t-s) under the integral for s being non-positive.

[José]

@Walter Roberson I'm using the standard notation for dde. Thank you very much Torsten, I just have a last question:

1) How do we put several different distributed delays into the equation? (or even a list of others delays?), i.e., by using:

lags1 = linspace(0,2r_1,number_of_delays_1);

lags2 = linspace(0,2r_2,number_of_delays_2);

lags3= [t-cos(t)

t-pi/2

t-1];

It is possible to do this with ddesd? Thanks in advance.

[Torsten]

The delay t-cos(t) is dubious because cos(t) easily becomes negative.

If you define

lags1 = linspace(0,2r_1,number_of_delays_1);
lags2 = linspace(0,2r_2,number_of_delays_2);
lags3= @(t,y)[t-cos(t),t-pi/2,t-1];
lags = @(t,y)[lags1,lags2,lags3(t,y)]

it should be no problem to use "ddesd" for it (assuming cos(t) remains non-negative).

But you will have to think about the "trapz" integration for time-dependent delays.

[José]

Thanks a lot Torsten. The list lags3 is generated to include both constant or variable delays (as you say must be non-negative deviations) into the model. I guess that to do this I should consider the order of the new list as:

lags = @(t,y)[lags3(t,y),lags1,lags2]

If we call Z(:,1) this will choose the first delay of the list ,i.e., y(t-cos(t)) (assuming cos(t) remains non-negative). So there is no need of think about "trapz" integration for time-dependent delays (we aboid this by separating "normal" delays from distributed delays)

Thank you Torsten. Best regards.

[Torsten]

I made a mistake on how to pass the delay lists to "ddesd". It should read

lags1 = linspace(0,2r_1,number_of_delays_1);
lags2 = linspace(0,2r_2,number_of_delays_2);
lags3 = @(t)[t-cos(t),t-pi/2,t-1];  % For your example, there is no dependency on y
lags = @(t,y)[t-lags1,t-lags2,lags3(t)] % t must be included for lags1 and lags2

So there is no need of think about "trapz" integration for time-dependent delays (we aboid this by separating "normal" delays from distributed delays)

So distributed delays in your notation only refer to a list of constant delays like lags1 and lags2 ?

[José]

lags1 and lags2 refers to the linspaces to define the distributed delays into the equation (I think it is ok as you wrote before). On the other hand, lags3 refer to a list of constant delays (not necessarily generated by a linspace) or well defined time-dependent delays to include into the model. For example (if the model is scalar):

where stands for some kernels. The delays and are selected from lags3, while lags1 and lags2 define the distributed delays into the model.

[Torsten]

lags1 and lags2 refers to the linspaces to define the distributed delays into the equation (I think it is ok as you wrote before).

Since you have to use ddesd, the lags are not passed as the (constant) values, but as the complete argument to y, thus as t - linspace(...).

[José]

Got it. Nevertheless, the lists lags1 and lags2 will be used just to define the distributed delays with the "trapz" rule. Not to be called as delays into the equation (to do this I will use lags3). I understand that your last correction is to define properly: yp = ddefunc(~,y,Z,lags)? It is necessary that @Fares consider this in his code?

[Torsten]

I understand that your last correction is to define properly: yp = ddefunc(~,y,Z,lags)? It is necessary that @Fares consider this in his code?

I don't understand what you mean here.

For "ddesd", you need to define a function handle in which you define the delays. For the delays given,

lags = @(t,y)[t-lags1,t-lags2,lags3(t)]

can serve as this function handle.

For use in the delay function, it might be better to pass each array of delays separately because you have to access them separately to compute the integrals over the distributed delays:

yp = ddefunc(t,y,Z,lags1,lags2,lags3)

[José]

I mean that there is no need in create: lags = @(t,y)[t-lags1,t-lags2,lags3(t)], since we don't use lags1 and lags2 as a list of delays in the model (as you say is better to access them separately). For example (if the model is scalar):

sol = ddesd(@(t,y,Z)ddefunc(t,y,Z,lags1,lags2,@lags3),@lags3,initialcondition,tspan);

dydt = ddefunc(t,y,Z,lags1,lags2,lags3)

yl1= trapz(lags1,Z(:)); % First distributed delay

yl2 = trapz(lags2,Z(:)); % Second distributed delay

ylag1 = Z(1); % Evaluate the solution in the first delay of the list lags3

ylag2 = Z(2); % Evaluate the solution in the second delay of the list lags3

dydt =-ylag1+yl1-ylag2+yl2;

My question is that if we put lags3 at the end (to define sol and dydt): Z(1), Z(2), etc. Will evaluate the solution in the list of delays lags3?

[Torsten]

It's necessary to know the values of y(t-lags1) and y(t-lags2) to evaluate yl1 and yl2. So it does not suffice to only pass lags3 to "ddesd".

Why else should we need

number_of_delays = 11;   % should be odd
lags = linspace(tau-gamma,tau+gamma,number_of_delays);
sol=ddesd(@(t,y,Z)ddefunc(t,y,Z,lags),lags,[0.2; 0.08],tspan);

in the code above ?

If it were as you said, we could just use

number_of_delays = 11;   % should be odd
lags = linspace(tau-gamma,tau+gamma,number_of_delays);
sol=ddesd(@(t,y,Z)ddefunc(t,y,Z,lags),[],[0.2; 0.08],tspan);

[José]

Dear Torsten, to clarifiy all, the code should be corrected as (omitting details for a scalar equation):

sol = ddesd(@(t,y,Z)ddefunc(t,y,Z,@lags),@lags,initialcondition,tspan);
lags1 = linspace(0,2r_1,number_of_delays_1);
lags2 = linspace(0,2r_2,number_of_delays_2); 
lags3= @(t,y)[t-(cos(t)+2),t-pi/2,t-1];
lags = @(t,y)[t-lags1,t-lags2,lags3(t,y)];
¿dydt = ddefunc(t,y,Z,@lags); or dydt = ddefunc(t,y,Z,lags1,lags2,@lags3)? 
yl1= trapz(lags1,Z(:)); % First distributed delay
yl2 = trapz(lags2,Z(:)); % Second distributed delay
ylag1 = Z(1); % ¿Evaluate the solution in the first delay of the list lags3?
ylag2 = Z(2); % ¿Evaluate the solution in the second delay of the list lags3?
dydt =-ylag1+yl1-ylag2+yl2;

How do we call the delays from lags3 into the model (to define ylag1 and ylag2)? Thanks in advance. Best regards.

[Torsten]

lags1 = linspace(0,2r_1,number_of_delays_1);
lags2 = linspace(0,2r_2,number_of_delays_2); 
lags3= @(t)[t-(cos(t)+2),t-pi/2,t-1];
lags = @(t,y)[t-lags1,t-lags2,lags3(t)];
sol = ddesd(@(t,y,Z)ddefunc(t,y,Z,lags1,number_of_delays_1,lags2,number_of_delays_2,lags3),lags,initialcondition,tspan);

function dydt = ddefunc(t,y,Z,lags1,number_of_delays_1,lags2,number_of_delays_2,lags3)
  yl1 = trapz(lags1,Z(:,1:number_of_delays_1)); % First distributed delay
  yl2 = trapz(lags2,Z(:,number_of_delays_1+1:number_of_delays_1+number_of_delays_2)); % Second distributed delay
  ylag1 = Z(:,number_of_delays_1+number_of_delays_2+1); % ¿Evaluate the solution in the first delay of the list lags3?
  ylag2 = Z(:,number_of_delays_1+number_of_delays_2+2); % ¿Evaluate the solution in the second delay of the list lags3?
  dydt =-ylag1+yl1-ylag2+yl2;
end

[José]

Thanks a lot Torsten.