Fonction de résolution pour la partie réelle = 0 au lieu de l'imaginaire dans MATLAB

Question

J'ai une fonction qui produit un vecteur de valeurs propres complexes. Il prend un seul argument, rho. J'ai besoin de trouver un rho pour lequel les valeurs propres complexes se trouvent sur l'axe imaginaire. En d'autres termes, les parties réelles doivent être 0.

Quand je lance fzero() il jette l'erreur suivante

Opérandes au || et & & Les opérateurs doivent être convertibles en valeurs scalaires logiques.

Alors que fsolve() résout simplement pour la partie imaginaire = 0, ce qui est exactement l'opposé de ce que je veux.

Ceci est la fonction I écrit

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) 

Il délivre en sortie 3, 2 valeurs propres complexes conjugués et une valeur propre réelle (avec une partie complexe = 0). Je dois trouver rho pour lesquels les valeurs propres complexes se trouvent sur l'axe imaginaire de sorte que les parties réelles sont 0.

Answer 1

Deux problèmes:

1. fzero ne convient pour les fonctions d'une valeur scalaire (f: ℝ → ℝ)
2. Les nombres complexes sont des nombres simples, traités comme des entités signées par presque toutes les fonctions. Vous devrez forcer Matlab à diviser le nombre complexe dans ses parties imaginaires et réelles

Ainsi, une solution de contournement possible est de prendre la partie réelle du premier complexe de valeurs propres:

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 

Appel comme ceci:

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