# How to store a variable in a vector?

15 views (last 30 days)
Pedro Sandoval on 15 Jul 2019
Commented: Pedro Sandoval on 19 Jul 2019
Hola, he tenido una serie de problemas con le siguiente código:
%Función para el cálculo del periodo de un sistema L-V
function calculo_periodo(~)
%***********************************************************%
%Termina cuando y retorna al punto incial (donde su ángulo
%con mu es el mismo que el ángulo entre eta y mu)
function [g,isterm,dir] = pits(t,y)
sig = sign(eta(1)-mu(1)+realmin);
theta1 = atan2(sig*(y(2)-mu(2)),sig*(y(1)-mu(1)));
theta0 = atan2(sig*(eta(2)-mu(2)),sig*(eta(1)-mu(1)));
g = theta1 - theta0;
isterm = t > 1;
dir = 1;
end
%************************************************************%
%Ciclo anidado para el cálculo del periodo, donde j=C e i=eta1
%for j=205:220
for i=0.001:0.001:0.999999
mu = [1 1]; % Punto de Equilibrio.
eta = [i 1]; % Condiciones Iniciales.
F = @(t,y) [(1-(y(2))/mu(2))*(y(1));
-220*(1-(y(1))/mu(1))*(y(2))];
opts = odeset('reltol',1.e-8,'events',@pits);
[t,~,~] = ode45(F,[0 Inf],eta,opts);
A = [i,t(end)]
%************************************************************%
%Generación de gráficos
hold on
plot(i,t(end),'.b')
hold off;
%end
end
A
end
the idea of the code is to find the period function of a lotka volterra system. For this, it is necessary to store the value of t (end) inside a vecto. How can I save t (end) in a vector?
If someone can help me, I appreciate it.
Greetings.
Pedro.

Jan on 15 Jul 2019
iVec = 0.001:0.001:0.999999;
result = zeros(1, numel(iVec));
for k = 1:numel(iVec)
i = iVec(k);
result(k) = t(end);
end

Pedro Sandoval on 17 Jul 2019
Excuse me Jan, I have tried to parameterize (with exponentials) the differential equation with m and n=0.5, I show you:
function periodo
function [g,isterm,dir] = pits(t,y)
sig = sign(eta(1)-mu(1)+realmin);
theta1 = atan2(sig*(y(2)-mu(2)),sig*(y(1)-mu(1)));
theta0 = atan2(sig*(eta(2)-mu(2)),sig*(eta(1)-mu(1)));
g = theta1 - theta0;
isterm = t > 1;
dir = 1;
end
iVec = 0.001:0.001:0.999999;
result = zeros(1, numel(iVec));
for m =0.1:0.1:0.9
for k = 1:numel(iVec)
i = iVec(k);
mu = [1 1]; % Punto de Equilibrio.
eta = [i 1]; % Condiciones Iniciales.
F = @(t,y) [(1-(y(2))^(m)/mu(2))*(y(1))^(0.5);
-10.3*(1-(y(1))^(0.5)/mu(1))*(y(2))^(m)];
opts = odeset('reltol',1.e-8,'events',@pits);
[t,~,~] = ode45(F,[0 Inf],eta,opts);
result(k) = t(end);
%************************************************************%
end
hold on
plot(0.001:0.001:0.999999,result)
hold off
end
%A(:,1) = 0.001:0.001:0.999999;
%A(:,2) = result;
%dlmwrite('data_2.txt',A,',');
end
but, I receive the following error message:
Error using atan2
Inputs must be real.
Error in periodo/pits (line 7)
theta1 = atan2(sig*(y(2)-mu(2)),sig*(y(1)-mu(1)));
Error in odezero (line 28)
[vnew,isterminal,direction] = feval(eventfun,tnew,ynew,eventargs{:});
Error in ode45 (line 406)
odezero(@ntrp45,eventFcn,eventArgs,valt,t,y,tnew,ynew,t0,h,f,idxNonNegative);
Error in periodo (line 26)
[t,~,~] = ode45(F,[0 Inf],eta,opts);
the function spits is for determinate the time of the revolution (the period). How can I improve my code to do not have de error mensaje?
Again, thank you for yout help!
Jan on 18 Jul 2019
@Pedro: Use the debugger to find out, why the values are not real. Insert some code to check, if the trajectories are real in the function to be integrated. Then set a breakpoint to stop Matlab.
The debugging is impeded massively, if you use anonymous functions. Prefer functions instead:
function dy = F(t,y)
dy = [(1-(y(2))^(m)/mu(2))*(y(1))^(0.5);
-10.3*(1-(y(1))^(0.5)/mu(1))*(y(2))^(m)];
if ~isreal(dy)
disp('not real') % Set a break point here
end
end
Pedro Sandoval on 19 Jul 2019
Hello Jan, again, thank you for yout help!