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 CenterFile Exchange에서 Loops and Conditional Statements에 대해 자세히 알아보기

태그

Community Treasure Hunt

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

Start Hunting!

Translated by