# How to store a variable in a vector?

21 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!

Pedro Sandoval on 17 Jul 2019
Thank you very much Jan for your help.