CODEWITHSUNDEEP
matlabnumerical-methodsnumerical-analysisnumerical-computing
Raul Guarini Riva
Raul Guarini Riva

Reputation: 795

Function fzero in Matlab is not converging

I am solving a problem in my Macroeconomics class. Consider the following equation:

enter image description here

Here, k is fixed and c(k) was defined through the ```interp1''' function in Matlab. Here is my code:

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

When I run this code, for some indexes $i$, it runs fine and fzero can find the solution. But for indexes it just returns NaN and exits without finding the root. This is a somewhat well-behaved problem in Economics and the solution we are looking indeed exists and the algorithm I tried to implement is guaranteed to work. But I don't have much experience with solving this in MATLAB and I guess I have a silly mistake somewhere. Any ideas on how to procede?

This is the typical error message:

Exiting fzero: aborting search for an interval containing a sign change because complex function value encountered during search. (Function value at -2.61092 is 0.74278-0.30449i.) Check function or try again with a different starting value.

Thanks a lot in advance!

Upvotes: 1

Views: 311

Answers (1)

Lutz Lehmann
Lutz Lehmann

Reputation: 25992

The only term that can produce complex numbers is 

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

You probably want the result according to the real variant of the cube root, which you could get as 

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

or more generally and avoiding divisions by zero

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

Upvotes: 1

Related Questions