MATLAB: Complex numbers and Nan’s as output using ‘fsolve’

complex numbersfsolve

 Hi,
I hope someone will spare a few minutes, I am getting complex numbers (NaN's sometimes) as results and wondering if I have incorrect code for the function and/or my use of fsolve is not apporriate here.
clc
clear all
format compact
global a b tau A
I=200;
tau=0.07;
a=0.2;
b=1.5;
Ac=zeros(I+1,1);
Aw=zeros(I+1,1);
Nc=zeros(I,1);
Nw=zeros(I,1);
Theta=zeros(I,1);
Ac0=100;
Aw0=500;
Ac(1,1)=Ac0;
Aw(1,1)=Aw0;
 x0=[0.7,0.3,0.1];
   for i=1:I 
    A=[Ac(i,1),Aw(i,1)];  
    fun=@(x)productivity3(x);
    result=fsolve(fun,x0);
    Nc(i,1)=result(1);
    Nw(i,1)=result(2);
    Theta(i,1)=result(3);
    Ac(i+1,1)=(1+tau*Theta(i,1))*Ac(i,1);
    Aw(i+1,1)=(1+tau*Theta(i,1))*Aw(i,1);
     end
     out=[Nc(I-5:I,1), Nw(I-5:I,1), Ac(I-5:I,1), Aw(I-5:I,1)]
Function
   function F=productivity3(x)      
      global a b tau A
   F=[x(1)+x(2)-1;
       x(1)-(x(3)/a)^(1/b)*(1+tau)*A(1);
       x(2)-(x(3)/a)^(1/b)*(1+tau)*A(2)];
Maybe it is due to choice of parameter values? Thank you. R

Best Answer

Replace your function by
      function F=productivity3(x) 
      global a b tau A
      F=[x(1)+x(2)-1;
       x(1)-(x(3)^2/a)^(1/b)*(1+tau)*A(1);
       x(2)-(x(3)^2/a)^(1/b)*(1+tau)*A(2)];
This way, it is not possible for fsolve to take the power (^1/b) of a negative number.
Best wishes
Torsten.
Related Question