Solve function for real part = 0 instead of imaginary in MATLAB

Question

I have a function which outputs a vector of complex eigenvalues. It takes a single argument, rho. I need to find a rho for which the complex eigenvalues lie on the imaginary axis. In other words, the real parts have to be 0.

When I run fzero() it throws the following error

Operands to the || and && operators must be convertible to logical scalar values.

Whereas fsolve() simply solves for the imaginary part = 0, which is exactly the oposite of what I want.

This is the function I wrote

function lambda = eigLorenz(rho)
beta = 8/3;
sigma = 10;
eta = sqrt(beta*(rho-1));
A = [ -beta 0 eta;0 -sigma sigma;-eta rho -1];
y = [rho-1; eta; eta];

% Calculate eigenvalues of jacobian
J = A + [0 y(3) 0; 0 0 0; 0 -y(1) 0]

lambda = eig(J)

It outputs 3 eigenvalues, 2 complex conjugates and 1 real eigenvalue (with complex part = 0). I need to find rho for which the complex eigenvalues lie on the imaginary axis so that the real parts are 0.

Answer 1

Two problems:

1. fzero is only suited for scalar-valued functions (f: ℝ → ℝ)
2. complex numbers are single numbers, treated as signle entities by almost all functions. You'll have to force MATLAB to split up the complex number into its imaginary and real parts

So, one possible workaround is to take the real part of the first complex eigenvalue:

function [output, lambda] = eigLorenz(rho)

    % Constants
    beta  = 8/3;
    sigma = 10;

    eta = sqrt(beta*(rho-1));
    A = [-beta        0     eta
             0   -sigma   sigma
          -eta      rho      -1];

    y = [rho-1
         eta
         eta];

    % Calculate eigenvalues of jacobian
    J = A + [0 y(3)  0
             0    0  0
             0 -y(1) 0];

    lambda = eig(J);

    % Make it all work for all rho with FZERO(). Check whether:
    % - the complex eigenvalues are indeed each other's conjugates   
    % - there are exactly 2 eigenvalues with nonzero imaginary part
    complex = lambda(imag(lambda) ~= 0);
    if numel(complex) == 2 && ...
            ( abs(complex(1) - conj(complex(2))) < sqrt(eps) )

        output = real(complex(1));   

    else
        % Help FZERO() get out of this hopeless valley
        output = -norm(lambda);
    end

end

Call like this:

rho = fzero(@eigLorenz, 0);
[~, lambda] = eigLorenz(rho);