How do I use a loop variable in the initial conditions in ODE45?
1 回表示 (過去 30 日間)
古いコメントを表示
I am trying to use loops variable in the initial conditions in ODE45.
function Poincaresection
clear
Ke=0.80;
omega=2*pi/6;
options=odeset('RelTol',1e-15,'AbsTol',1e-15);
for x0=0.1:0.1:0.9;
for y0=0.1:0.1:0.9;
[t,X]=ode45(@Preypred,0:(2/omega)*pi:(4000/omega)*pi,[1.5;x0;y0]);
%only plot last half of solutions to avoid messy transient dynamics
sh=round(length(t)*.05); %sh=secondhalf of solutions
figure(2)
subplot(2,1,1)
plot(X(sh:end,1),X(sh:end,3),'k.','MarkerSize',1)
%plot(X(:,1),X(:,3),'k.','MarkerSize',1)
fsize=15;
axis([0 2 0 2])
xlabel('K','FontSize',fsize)
ylabel('y','FontSize',fsize)
title('Poincare Section of the LVE System')
hold on
subplot(2,1,2)
plot(X(sh:end,2),X(sh:end,3),'k.','MarkerSize',1)
%plot(X(:,2),X(:,3),'k.','MarkerSize',1)
fsize=15;
axis([0 2 0 2])
xlabel('x','FontSize',fsize)
ylabel('y','FontSize',fsize)
title('Poincare Section of the LVE System')
hold on
end
end
function dX=Preypred(t,X)
dX(1)=0.5*(Ke-X(1))+(Ke*omega*cos(omega*t)/2);
dX(2)=1.2*X(2)*(1-X(2)/(min(X(1),(0.025-0.03*X(3))/0.0038)))-0.81*X(2)*X(3)/(0.25+X(2));
dX(3)=0.8*min(1,((0.025-0.03*X(3))/X(2))/0.03)*0.81*X(2)*X(3)/(0.25+X(2))-0.25*X(3);
dX=[dX(1);dX(2);dX(3)];
end
end
The code gives me the following error:
Warning: Failure at t=2.503199e-01. Unable to meet integration tolerances without reducing the step size below the smallest value
allowed (8.881784e-16) at time t.
> In ode45 (line 308)
In Poincaresection4 (line 16)
Subscript indices must either be real positive integers or logicals.
Error in Poincaresection4 (line 26)
plot(X(sh:end,1),X(sh:end,3),'k.','MarkerSize',1)
Could you please help me to fix it?
0 件のコメント
採用された回答
John D'Errico
2017 年 10 月 24 日
No. As you did it, it is not possible. You are trying to solve multiple problems in one call to ODE45. Instead, just learn to use a for loop. Solve each problem one at a time.
0 件のコメント
その他の回答 (0 件)
参考
カテゴリ
Help Center および File Exchange で Ordinary Differential Equations についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!