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

1 visualización (últimos 30 días)
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?

Respuesta aceptada

John D'Errico
John D'Errico el 24 de Oct. de 2017
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.

Más respuestas (0)

Categorías

Más información sobre Ordinary Differential Equations en Help Center y File Exchange.

Etiquetas

Community Treasure Hunt

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

Start Hunting!

Translated by