Функция fzero в Matlab не сходится - PullRequest
Новая
       76

Функция fzero в Matlab не сходится

1 голос
/ 17 апреля 2020

Я решаю проблему в моем классе макроэкономики. Рассмотрим следующее уравнение:

enter image description here

Здесь k фиксировано и c (k) определено с помощью функции `` `interp1 '' ' в Matlab. Вот мой код: 

beta = 0.98;
delta = 0.13;
A = 2;
alpha = 1/3;
n_grid = 1000;   % Number of points for capital
k_grid = linspace(5, 15, n_grid)';
tol = 1e-5;
max_it = 1000;

c0 = ones(n_grid, 1);
new_k = zeros(n_grid, 1);
dist_c = tol + 1;
it_c = 0;

while dist_c > tol && it_c < max_it
   c_handle = @(k_tomorrow) interp1(k_grid, c0, k_tomorrow, 'linear', 'extrap');

   for i=1:n_grid
       % Solve for k'
       euler = @(k_tomorrow) (1/((1-delta)* k_grid(i) + A * k_grid(i)^alpha - k_tomorrow)) - beta*(1-delta + alpha*A*k_tomorrow^(alpha - 1))/c_handle(k_prime);
       new_k(i) = fzero(euler, k_grid(i));    % What's a good guess for fzero?    
   end

   % Compute new values for consumption
   new_c = A*k_grid.^alpha + (1-delta)*k_grid - new_k;

   % Check convergence
   dist_c = norm(new_c - c0);
   c0 = new_c;
   it_c = it_c + 1;

end

Когда я запускаю этот код, для некоторых индексов $ i $ он работает нормально, и fzero может найти решение. Но для индексов он просто возвращает NaN и выходит, не найдя root. В экономике это проблема с хорошим поведением, и решение, которое мы ищем, действительно существует, и алгоритм, который я пытался реализовать, гарантированно сработает. Но у меня нет большого опыта в решении этой проблемы в MATLAB, и я думаю, что где-то глупая ошибка. Любые идеи о том, как действовать?

Это типичное сообщение об ошибке:

Выход из fzero: прерывание поиска интервала, содержащего изменение знака, поскольку во время поиска возникла сложная функция. (Значение функции в -2.61092 равно 0.74278-0.30449i.) Проверьте функцию или попробуйте снова с другим начальным значением.

Заранее большое спасибо!

1 Ответ

0 голосов
/ 17 апреля 2020

Единственный термин, который может образовывать комплексные числа: 

k'^(alpha - 1) = k'^(-2/3)

Вы, вероятно, хотите получить результат в соответствии с реальным вариантом куба root, который можно получить как 

sign(k') * abs(k')^(-2/3)

или более в целом и избегая деления на ноль 

k' * (1e-16+abs(k'))^(alpha - 2)
...