Delay differential equations where every variable can take on multiple delays

Hello,

I'm trying to code the Kuramoto model as a delay differential equation. I have a connectivity matrix and a delay matrix. So the equation looks like this:

The problem is the delay matrix D has 0s on the diagonal and when formulating the DDE all the time delays has to be positive or I get the error message "The lags must all be positive". One way to circumvent this is to do an if / else loop:

if i ~= j
    theta_j = Z(:, ???)
else
    theta_j = theta(j);
end

Then the problem is how to formulate lags and Z. One can do the following:

lags = D(find(triu(D, 1)));

Which would only give the non-diagonal upper triangle of D. Then one would need to find a mapping between ij and lags to put in as a Z index. I am inclined to use sub2ind but they would assume that the matrix is symmetrical and square, not only the upper triangle. Is there a way to make this work?

For context, here is the full code for the Kuramoto equation:

function dtheta = krm_dde(t, theta, Z, p)

for i = 1:length(p.omega)
    sum_coupling = 0;
    % Do the sigma term
    for j = 1:length(p.omega)
        if i ~= j
            theta_j = Z(:, ???);
        else
            theta_j = theta(j);
        end
    sum_coupling = sum_coupling + p.C(i, j) * sin(theta_j - theta(i));
    end
    % Do the RHS
    dtheta(i) = p.omega(i) + K * sum_coupling;
end
end

Tangentially related: I would also appreciate any tips to make this code vectorized. It is horribly slow in this state.

Thanks,

Yasir

Answer 1

[Torsten]

I don't understand your mathematical formulation. Should the summation be over j instead of i ? And shouldn't there be a bracket around the sin arguments ?

If this is the case, I guess the following should work if you define the lags as tau(1)=D12,tau(2)=D13, tau(3)= D14,...,tau(n-1)=D1n,tau(n) = D21,tau(n+1) = D23,... etc.

delay_number = 0;
for i = 1:length(p.omega)
    sum_coupling = 0;
    % Do the sigma term
    for j = 1:length(p.omega)
        if j ~= i
            delay_number = delay_number + 1;
            theta_j = Z(j, delay_number);
        else
            theta_j = theta(j);
        end
    sum_coupling = sum_coupling + p.C(i, j) * sin(theta_j-theta(i));
    end
    % Do the RHS
    dtheta(i) = p.omega(i) + K * sum_coupling;
end

[Yasir Çatal]

You're right about the mathematical formulation, sorry; sleep deprivation is not good for math. Editing the formula now. I will check your recommendation today.

[Yasir Çatal]

Thanks! I implemented the following and it worked:

tau = zeros(length(D) * (length(D)-1), 1);
c = 1;
for i = 1:length(D)
    for j = 1:length(D)
        if i~=j
            tau(c) = D(i, j);
            c = c + 1;
        end
    end
end