Not getting the right output using 'fsolve'

Hi all,

I am trying to practice with 'fsolve' and have not figured out fully what's been going on with the following code. Can anyone please shed some light on it?

function N=productivity1(N,Ac,Aw)
global Thetac Thetaw tau a b

N=[N(1),N(2)];
N=[N(1)-(Thetac/a)^(1/b)*(1+tau)*Ac; 
  N(2)-(Thetaw/a)^(1/b)*(1+tau)*Aw;
  N(1)+N(2)-1]; %this meant to be a constraint...
end

N0=[0.7,0.3]; %initial guess for x
option=optimset('Display','iter');
result=fsolve(@(N)productivity1(N,Ac0,Aw0),N0,option);*

I must be coding it incorrectly. Much appreciated! Thank you

[Matt J]

We cannot know what you don't like about the results unless you tell us. We also can't test the code ourselves without knowing

Thetac Thetaw tau a b Ac Aw

[mrrox]

Thanks Matt,

Here are the details:

• Ac=2
• Aw=20

Ac and Aw should take values such that N(1)+N(2)=1, the starting values of Ac and Aw are not fixed

Parameter values

• Thetac=0.6
• Thetaw=0.6
• tau=0.07
• a=0.01
• b=0.8

The output is the same as the vector of initial values, regardless of value of the above parameters.

[mrrox]

With the last sentence, I meant to say that the output should not be the same as the inputs. This is the issue I am facing.

Answer 1

[Matt J]

R,

Below is what I think you're trying to implement. Note that, according to your own description, Ac and Aw are unknowns so you shouldn't be treating them as constants, as you are in your original post. Further, it makes little sense to "solve" for N1 and N2 when they depend so trivially on Ac and Aw. I've removed them accordingly and simply solve for Ac and Aw.

Your practice example is a doubtful one, however, because there are more unknowns than equations, leading to infinite solutions. That would be true whether you include N1 and N2 as unknowns or not. However, fsolve does succeed in finding one of the infinite solutions.

    q1=(Thetac/a)^(1/b)*(1+tau);
    q2=(Thetaw/a)^(1/b)*(1+tau);

    option=optimset('Display','iter');

    fun=@(x)productivity1(x,q1,q2);
    x0=[Ac,Aw]; %initial guess for x

    result=fsolve(fun,x0,option)

    function F=productivity1(x,q1,q2)

      N1=q1*x(1);
      N2=q2*x(2);

      F=N1+N2-1;

[mrrox]

Many thanks Matt! That's great, this simplified the thing that i was trying to understand. I will rewrite the code and see if it works this time round. My knowledge of matlab is limited to econometrics applications.

[mrrox]

Matt,

it works flawless, I try to use the code on a similar function now. Just a quick question, how will I retrieve the values for N?

[Matt J]

just post-calculate them using the equations defining them,

        N1=q1*result(1);
        N2=q2*result(2);

or in vectorized form

      N=result.*[q1,q2]

[mrrox]

Thanks a ton Matt, it is going well now!

[mrrox]

Matt,

Is there a way to insert the following constraint in the above function?

0<N1<1; 
0<N2<1;
N1+N2<=1;

Thanks!

[Matt J]

You could make the change of variables

x(1)=sin(t)^2/q1;
x(2)=cos(t)^2/q2;

However, linear constraints like these are more generally handled using fmincon.

[mrrox]

Matt, you answers made me realise where I made the error in coding. you have simplied the problem greately this helped to learn how to code neatly. Will you please have a look at the new version of the code below? Thanks.

Answer 2

[mrrox]

Hi Guys,

I've managed to get it right, I've been practicing ignoring a number of important issues you raised in your posts. I just did not realise that I did not have equal no of unknowns and equations. I made one parameter an unkonwn and it worked as it should. Here is the code. The additional variable is x(3) and x(1)+x(2)=1 condition is still there.

clc
clear all

format longE
global a b tau A
I=100;
tau=0.07;
a=0.02;
b=1.2;

Ac=zeros(I+1,1); 
Aw=zeros(I+1,1);
Nc=zeros(I,1);
Nw=zeros(I,1);
Theta=zeros(I,1);

Ac0=3;
Aw0=2;
Ac(1,1)=Ac0;
Aw(1,1)=Aw0;
x0=[0,0,0.5];

for i=1:I 
    A=[Ac(i,1),Aw(i,1)];  
    option=optimset('Display','iter');    
    fun=@(x)productivity3(x);
    result=fsolve(fun,x0,option);

      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 

 *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)];

I still did not get the desired result, I thought as A(1) (or Ac) increases, x(1) which is Nc would increase. The x(1) and x(2) remain constant for some reason with little or no variation. I am just wondering if I have a wrong code or incorrectly coded it again.