2
$\begingroup$

I am a beginner in learning DDE and i am trying to understand how to use matrix notation to simulate basic systems of the specific form below: Let us assume a directed network with three nodes

$$\dot{y_{1}} = a_{11}y_{1} + a_{12}y_{2}(t-\tau_{12}) + a_{13}y_{3}(t-\tau_{13})$$

$$\dot{y_{2}} = a_{21}y_{1}(t-\tau_{21}) + a_{22}y_{2} + a_{23}y_{3}(t-\tau_{23})$$

$$\dot{y_{3}} = a_{31}y_{1}(t-\tau_{31}) + a_{32}y_{2}(t-\tau_{32}) + a_{33}y_{3}$$

The direct. way to code in Matlab for the above, is:

tau = [1.0 0.5 1.5];
tf = 10;
sol = dde23(@dde,tau,@history,[0 tf]);

t = linspace(0,tf,200);
y = deval(sol,t);

figure
subplot(1,2,1)
plot(t,y)
grid on;

function dydt = dde(t,y,yd)

a11 = -1;a12 = 2; a13 = 1;
a21 = 2; a22 = -1;a23 = 1;
a31 = 1; a32 = 1; a33 = -1;
 
dydt = [a11*y(1)    + a12*yd(2,1) + a13*yd(3,2)
        a21*yd(1,1) + a22*y(2)    + a23*yd(3,3)
        a31*yd(1,2) + a32*yd(2,3) + a33*y(3)];

end

function y = history(t)

y = [1;0;-1];

end

the result is this: result with the direct implementation

But for a n by n system, the process should be somehow automatized to produce the above $n(n-1)/2$ matrices plus the first one which is the diagonal:

$$\begin{multline} \begin{bmatrix} \dot{y_{1}}\\ \dot{y_{2}}\\ \dot{y_{3}}\\ \end{bmatrix} = \begin{bmatrix} a_{11} & 0 & 0\\ 0 & a_{22} & 0 \\ 0& 0 & a_{33} \end{bmatrix} \begin{bmatrix} y_{1}\\ y_{2}\\ y_{3}\\ \end{bmatrix} + \begin{bmatrix} 0 & a_{12} & 0\\ a_{21} & 0 & 0 \\ 0& 0 & 0 \end{bmatrix} \begin{bmatrix} y_{1}(t-\tau _{1})\\ y_{2}(t-\tau _{1})\\ y_{3}(t-\tau _{1})\\ \end{bmatrix} \\ + \begin{bmatrix} 0 & 0 & a_{13}\\ 0 & 0 & 0 \\ a_{31} & 0 & 0 \end{bmatrix} \begin{bmatrix} y_{1}(t-\tau _{2})\\ y_{2}(t-\tau _{2})\\ y_{3}(t-\tau _{2})\\ \end{bmatrix} + \begin{bmatrix} 0 & 0 & 0\\ 0 & 0 & a_{23} \\ 0 & a_{32} & 0 \end{bmatrix} \begin{bmatrix} y_{1}(t-\tau _{3})\\ y_{2}(t-\tau _{3})\\ y_{3}(t-\tau _{3})\\ \end{bmatrix} \end{multline}$$

The code for the general case, applied in the same 3 by 3 is the following:

tau = [1.0 0.5 1.5];
tf = 10;
a11 = -1;a12 = 2; a13 = 1;
a21 = 2; a22 = -1;a23 = 1;
a31 = 1; a32 = 1; a33 = -1;

A = [a11 a12 a13;
     a21 a22 a23;
     a31 a32 a33];

sol = dde23(@(t,y,z)(dde(t,y,z,A)),tau,@(t,y)(history(t,A)),[0 tf],A);

t = linspace(0,tf,200);
y = deval(sol,t);

figure
subplot(1,2,1)
plot(t,y)
grid on;

function dydt = dde(t,y,yd,A)

 n = length(A);
 
 A0 = diag(diag(A));
 m = n*(n-1)/2;
 P = zeros(n,n*m);
 k = 0;
 for i = 1:n
     for j = 1:n
         if j>i
             P(i,j+n*k) = A(i,j);
             P(j,i+n*k) = A(i,j);
             k = k+1;
         end
     end
 end
 
for c = 0:m-1
    
    dydt = A0*y + P(:,n*c+1:(c+1)*n)*yd(:,c+1);

end

end

function y = history(t,A)

y = [1;0;-1];

end

result with the general case

$\endgroup$
9
  • $\begingroup$ This looks good. One point that has to be handled carefully is the order of the delays between the different implementations. Can you give the code implementing the matrix system? It should be something like A0*y+A1*yd(:,1)+A2*yd(:,2)+A3*yd(:,3). $\endgroup$ – Lutz Lehmann Mar 15 at 11:15
  • $\begingroup$ Yes, I can! Please, mind the following: The problem is about the general case where i have to calculate the above for an n by n. So, the matrix system cannot be manually introduced within the function. The code I will add above, is still for a 3 by 3 but, I am trying to automate the process, and I am obviously failing- don't know what I mess up. $\endgroup$ – Alex Kps Bdc Mar 15 at 11:24
  • $\begingroup$ You can combine for j=1:n if j>i to for j=(i+1):n. Is it ensured that the reference solution gets a matching delay sequence? $\endgroup$ – Lutz Lehmann Mar 15 at 11:43
  • $\begingroup$ I will reform the post with the codes and the graphs I receive $\endgroup$ – Alex Kps Bdc Mar 15 at 11:45
  • $\begingroup$ In this case the matrix approach might also not be the most efficient. It should be simpler and slightly faster to build dydt starting from dydt=diag(A).*y in a single double loop, instead of constructing P. $\endgroup$ – Lutz Lehmann Mar 15 at 11:58
1
$\begingroup$

It is obvious once you have found it, but at a superficial glance seems to have all the correct parts:

for c = 0:m-1 
    dydt = A0*y + P(:,n*c+1:(c+1)*n)*yd(:,c+1);
end

does not what you intended it to do. It should be

dydt = A0*y;
for c = 0:m-1 
    dydt = dydt + P(:,n*c+1:(c+1)*n)*yd(:,c+1);
end
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.