How do I use a loop variable in the initial conditions in ODE45?

1 回表示 (過去 30 日間)
LALE ASIK
LALE ASIK 2017 年 10 月 23 日
回答済み: John D'Errico 2017 年 10 月 24 日
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?

採用された回答

John D'Errico
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 件)

カテゴリ

Help Center および File ExchangeOrdinary Differential Equations についてさらに検索

タグ

Community Treasure Hunt

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

Start Hunting!

Translated by