How to Use Nested For Loops to Create a Coupled System of Equations

3 ビュー (過去 30 日間)
Cameron3332
Cameron3332 2022 年 8 月 4 日
編集済み: Torsten 2022 年 8 月 4 日
I am trying to write a system of coupled ODEs to eventually pass through ode45(). The equations are:
When I try the following code in an attempt to produce those equations and solve then for n = 3, I get the following error:
options = odeset('RelTol',3e-14,'AbsTol',1e-15,'Stats','on');
[t,y]=ode45(@hmmPD,[0 5000],[0.1 0.2 0.3],options);
figure
plot(t,y(:,1),'b')
hold on
plot(t,y(:,2),'r')
hold on
plot(t,y(:,3),'g')
title('Compare')
hold off
function dydt = hmmPD(t,y)
k=1/5;
a=1.25;
r=1/4;
n=3;
m=4;
for i = 1:n
for j = 1:n
dydt = [(k/m)*(sin(-y(i)+a)-r*sin(-2*y(i)) + sin(y(j)-y(i)+a)-r*sin(2*(y(j)-y(i))));];
end
end
end
Error using odearguments (line 95)
HMMPD returns a vector of length 1, but the length of initial conditions vector is 3. The vector
returned by HMMPD and the initial conditions vector must have the same number of elements.
Error in ode45 (line 106)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in LargeSystem (line 4)
[t,y]=ode45(@hmmPD,[0 5000],[0.1 0.2 0.3],options);
This appears, to me at least, to be a dimensionality issue. The for loops are likely not doing what I expect and either not constructing the equation I was attempting to or the size of the dydt vector is not right. Either way, I have tried quite a few things but haven't been able to fix it. Any suggestions?

採用された回答

Alan Stevens
Alan Stevens 2022 年 8 月 4 日
More like this perhaps (you've very tight tolerances!):
options = odeset('RelTol',3e-14,'AbsTol',1e-15,'Stats','on');
[t,y]=ode45(@hmmPD,[0 5000],[0.1 0.2 0.3],options);
21891 successful steps 9 failed attempts 131401 function evaluations
figure
plot(t,y(:,1),'b')
hold on
plot(t,y(:,2),'r')
hold on
plot(t,y(:,3),'g')
title('Compare')
hold off
function dydt = hmmPD(~,y)
k=1/5;
a=1.25;
r=1/4;
n=3;
m=4;
for i = 1:n
term = sin(-y(i)+a)-r*sin(-2*y(i));
summation = 0;
for j = 1:n
summation = summation + sin(y(j)-y(i)+a)-r*sin(2*(y(j)-y(i)));
end
dydt(i,1) = (k/m)*(term + summation);
end
end
  1 件のコメント
Torsten
Torsten 2022 年 8 月 4 日
編集済み: Torsten 2022 年 8 月 4 日
@Cameron3332 comment moved here:
Hey Alan, thanks so much for the answer - it did exactly what I was looking for! As far as the tolerances go, unfortunately the extended system has some pretty nasty heteroclinic orbits that have required me to knock the tolerances up a bit. That has also done the trick through.
Thanks again!

サインインしてコメントする。

その他の回答 (0 件)

製品


リリース

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by