How to store a variable in a vector?

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.

Answer 1

[Jan]

iVec = 0.001:0.001:0.999999;
result = zeros(1, numel(iVec));
for k = 1:numel(iVec)
    i = iVec(k);
    ... your code
    result(k) = t(end);
end

[Pedro Sandoval]

Hello Jan, thanks for your response. When you put "my code" is all my code, or is it a part of my code?

Thanks.

[Jan]

%Función para el cálculo del periodo de un sistema L-V
function calculo_periodo(~)
     
%************************************************************%
%Ciclo anidado para el cálculo del periodo, donde j=C e i=eta1
iVec = 0.001:0.001:0.999999;
result = zeros(1, numel(iVec));
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))/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

    result(k) = t(end);
end
A
end

%***********************************************************%
%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

[Pedro Sandoval]

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]

@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]

Hello Jan, again, thank you for yout help!

Answer 2

[Pedro Sandoval]

Thank you very much Jan for your help.