필터 지우기
필터 지우기

Delay differential equations where every variable can take on multiple delays

조회 수: 3 (최근 30일)
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

채택된 답변

Torsten
Torsten 2023년 10월 3일
편집: Torsten 2023년 10월 3일
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
  댓글 수: 2
Yasir Çatal
Yasir Çatal 2023년 10월 3일
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
Yasir Çatal 2023년 10월 3일
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

댓글을 달려면 로그인하십시오.

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Programming에 대해 자세히 알아보기

제품


릴리스

R2023a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by