MATLAB Как мне зациклить якобиана и построить собственные значения?

Question

У меня есть система трех уравнений: инфляция, выпуск и процент. Все переменные сильно взаимосвязаны в этих уравнениях, однако уравнение интереса является вспомогательным.

Первым шагом было определение устойчивых комбинаций инфляции и потребления (см. Первое изображение) и построение этой зависимости. Следующий код должен был это сделать, но он выдает ошибку о 'preallocation' c_sol.

Выше код должен дать мне 50 комбинаций инфляции и выхода; проценты затем рассчитываются как инфляция, деленная на бета.

Следующим шагом является оценка якобиана по каждому из этих значений, определение собственных значений и построение собственных значений с инфляцией на оси x. Для устойчивости требуется, чтобы собственные значения находились внутри единичной окружности, а в случае комплексных чисел - модуля. Вот код, который у меня есть, но он требует гораздо большего редактирования:

Помощь очень ценится!

https://imgur.com/a/Gx3hU

Обновление: это код, с которым я работаю, но мне не хватает части после 'for'. Мне нужно рассчитать числовую матрицу для каждого из этих 50 значений. 'c2' и 'p2' являются отстающими переменными и равны 'c' и 'p'.

syms c p R c2 p2

outp = (c2*p2)/(((99*R)/100 - 99/100)*(p2/p)^((2*R)/(R - 1))*(c2/c)^(R/(2*(R - 1))) + 99/100);

infl = (((2*c2*p2)/(35*(((99*R)/100 - 99/100)*(p2/p)^((2*R)/(R - 1))*(c2/c)^(R/(2*(R - 1))) + 99/100)) + 2/175)/(((99*R)/100 - 99/100)*(p2/p)^((2*R)/(R - 1))*(c2/c)^(R/(2*(R - 1))) + 99/100) - (99*p2)/100 + (99*p2^2)/100 + (3*((c2*p2)/(((99*R)/100 - 99/100)*(p2/p)^((2*R)/(R - 1))*(c2/c)^(R/(2*(R - 1))) + 99/100) + 1/5)^(20/7))/35 + 1/4)^(1/2) + 1/2;

J = [jacobian(outp, c2), jacobian(outp, p2) ; jacobian(infl, c2), jacobian(infl, p2)];

%%%%%%%%

beta = 0.99 ; v = 21 ; gamma = 350 ; eps = 1 ; g = 0.2 ; sigma = 1; alpha = 0.7;

c_sol = zeros(1,50);
p = linspace(0.8,1.2,50);
for     i = 1 : numel(p)
pi_actual = p(i);
fun = @(c)  (1 - beta)*pi_actual*(pi_actual-1) - (v/(alpha*gamma))*(c+g).^((1+eps)/alpha) - ((1-v)/gamma)*(c+g)*c.^(-sigma);
c_sol(i) = fzero(fun,0.5);
end

c = c_sol; c2 = c_sol; p2 = p;
R = p/beta; % interest 

%%%%%%%


eigVals = zeros(2,numel(p));

for k = 1:numel(p)


    eigVals(:,k) = eig(B);
end

Answer 1

Предполагая, что ваш вектор x имеет ту же длину, что и ширина матрицы eigVals, вы можете использовать следующий код для построения собственных значений для значения x:

figure;
hold on;
for i = 1:size(eigVals, 1)
    plot(x, eigVals(i, :), '*');
end
title('Eigenvalues against x');
xlabel('x');
ylabel('Eigenvalue');