# Solving DDE in Matlab using matrix form

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:

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


• 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). Mar 15, 2021 at 11:15
• 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. Mar 15, 2021 at 11:24
• 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? Mar 15, 2021 at 11:43
• I will reform the post with the codes and the graphs I receive Mar 15, 2021 at 11:45
• 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. Mar 15, 2021 at 11:58

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