2017-01-06 2 views
1

У меня есть функция, которая выводит вектор сложных собственных значений. Требуется один аргумент: rho. Мне нужно найти rho, для которого комплексные собственные значения лежат на мнимой оси. Другими словами, вещественные части должны быть 0.Решите функцию для действительной части = 0 вместо мнимой в MATLAB

Когда я бегу fzero() он бросает следующую ошибку

операнды к || и & Операторы & должны быть преобразованы в логические скалярные значения.

Принимая во внимание, что fsolve() просто решает для мнимой части = 0, что является абсолютно опосредованным из того, что я хочу.

Это функция я написал

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) 

Он выводит 3 собственных значений, 2 сложные конъюгаты и 1 вещественное собственное значение (со сложной части = 0). мне нужно найти rho, для которых комплексных собственных лежат на мнимой оси, так что действительные части равны 0.

+0

Каков допустимый диапазон ρ? Поскольку просто заполнение некоторых случайных значений не дает 1-real-and-2-conjugates, которые вы описываете ... это действительно требование? –

+0

@RodyOldenhuis Допустимый диапазон равен rho> 0. Решение rho = 24.737, но теперь мне нужен способ получить это, чтобы работать в MATLAB – Ortix92

+0

2 конъюгата + реальное значение начинается только при ρ> 1.34561 ... –

ответ

1

две проблемы:

  1. fzero является только подходит для скалярных функций (F: ℝ → ℝ)
  2. комплексные числа - это одиночные числа, которые рассматриваются как объекты привязки почти всеми функциями. Вы должны заставить MATLAB разделить комплексное число в его мнимые и действительные части

Таким образом, один можно обойти это взять действительную часть первого комплексного собственного значения:

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 

вызовов например:

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

I как раз собирался прокомментировать, что вам нужно получить 'lambda (1)', но вы избили меня;) – Ortix92

+0

... да, но это не очень здорово; это работает только до тех пор, пока ρ> 1.34561 ... Я все еще не счастлив! Будет еще несколько изменений :) –

+0

@ Ortix92 там вы идете; это работает для -∞ <ρ <+ ∞. –

Смежные вопросы