Search code examples
matlabdifferential-equations

Solving DDE in Matlab

I am trying to learn how to solve DDE (delay diff. eq) on Matlab and I am using a very helpful (Youtube-tutorial) where the guy solves examples. In the case of a 3-dimensional system, the code goes as follows:

tau = [1 0.5];
tf = 10;
sol = dde23(@dde,tau,@history,[0 tf]);
t = linspace(0,tf,200);
y = deval(sol,t);

figure(1)
plot(t,y)

function y = history(t)
y = [1;0;-1];
end

The dde function according to the tutorial is:

function dydt = dde(t,y,tau)
y1tau1 = tau(:,1);
y2tau2 = tau(:,2);
dydt = [y1tau1(1)
    y(1) - y1tau1(1) + y2tau2(2)
    y(2) - y(3)];
end

, where if I get it well, we have a 3by2 matrix with the state variables and their respective delays: the first column is about the first delay tau_1 for each of the three state variables, and the second column is for the tau_2. Nevertheless, this matrix has only 2 non-zero components and they are the tau_1 delay for y_1 and the tau_2 for y_2. Given that, I thought it would be the same (at least in this case where we only have one delay for y_1 and one for y_2) as if the function were:

function dydt = dde(t,y,tau)
dydt = [tau(1)
    y(1) - tau(1) + tau(2)
    y(2) - y(3)];
end

I ran the script for both of them and the results are totally different, both qualitatively and quantitatively, and I cannot figure out why. Could someone explain the difference?

Solution

  • y1tau1(1) and y2tau2(2) are tau(1,1) and tau(2,2)

    while

    tau(1) and tau(2) are the same as tau(1,1) and tau(2,1).

    So the second one is different, there is no reason why the two different, if only slightly so, DDE systems should have the same solutions.

    As you have also a general interpretation question, the DDE system is mathematically

    dy1(t)/dt = y1(t-1)
dy2(t)/dt = y1(t) - y1(t-1) + y2(t-0.5)
dy3(t)/dt = y2(t) - y3(t)

    It is unfortunate that the delayed values are also named tau, it would be better to use something else like yd,

    function dydt = dde(t,y,yd)
  dydt = [ yd(1,1)
           y(1) - yd(1,1) + yd(2,2)
           y(2) - y(3) ];

end

    This yd has two columns, the first is the value y(t-tau(1))=y(t-1), the second column is y(t-tau(2))=y(t-0.5). The DDE function does not know the delays themselves, it only computes with the values at these delays. These are obtained from a piecewise interpolation function that continues the history function with the solution so far. The history function takes the place of the initial condition.