Solving N non-linear equations using fsolve. How do I pass these equations into my function without typing them out individually?

Hello all,

I am currently working with the Eaton-Kortum Trade Model in MATLAB. In this model we have N countries, and wish to solve 2N + N^2 non-linear equations for equilibrium outcomes. I am working currently on an example with four countries, which means I will need to use fsolve to solve 24 equations. I understand that I could type all 24 equations individually, but what happens when we allow N to grow in the model (to better reflect what the world looks like)? If I wanted to consider trade between 10 countries I would have to type 120 equations seperatley! Luckily these equations take one of three forms.

N of the equations take the form: (gam.*((sum(Ti.*(dni.*(w(i).^(beta)).*(p(i))).^(1-beta)).^(-theta)).^(-1./theta))) - p(n);

N^2 of the equations take the form: Ti.*(gam.*dni.*(w(i).^(beta))*(p(i).^(1-beta))*(1./p(i))).^(theta) - (x(i));

N of the equations take the form: ((beta).*(sum(Ti.*(gam.*dni.*(w(i).^(beta))*(p(i).^(1-beta))*(1./p(i))).^(theta)*x(i)))) - (w(i).*Li)

Where our unknowns are w's, p's, and x's and everything else is given.

Is there a way for me to iteratively feed these equations into fsolve?

For example:

F(1) = (gam.*((sum(Ti.*(dni.*(w(i).^(beta)).*(p(i))).^(1-beta)).^(-theta)).^(-1./theta))) - p(1);
F(2) = (gam.*((sum(Ti.*(dni.*(w(i).^(beta)).*(p(i))).^(1-beta)).^(-theta)).^(-1./theta))) - p(2);
F(3) = (gam.*((sum(Ti.*(dni.*(w(i).^(beta)).*(p(i))).^(1-beta)).^(-theta)).^(-1./theta))) - p(3);
F(4) = (gam.*((sum(Ti.*(dni.*(w(i).^(beta)).*(p(i))).^(1-beta)).^(-theta)).^(-1./theta))) - p(4);

If there is not a way to do what I suggest how should I attempt to implement this?

[darova]

What about for loop?

[Ethan Goode]

@darova

That was my initial thought, but I am not sure how to implement it. Would it be along the lines of:

function F=ekmodel(p, w, x)
     n = 4
     for i = 1:n
        F(i) = (gam.*((sum(Ti.*(dni.*(w(n).^(beta)).*(p(n))).^(1-beta)).^(-theta)).^(-1./theta))) - p(i)
     end
end

x0 =

And so on?

I'm having trouble making sure the function above is then summing across all of the n observations, because in reality the function is indexed by both i and n. Where i is a certain country, and n is all other countries. The above function should represent 4 equations.

How can I better implement the for loop? Can I call a for loop like this in my function, and if so what is the appropriate way to do it?

[Walter Roberson]

Build a cell array of function handles.

fsolve(@(x) arrayfun(@(F)F(x), CellOfHandles), x0)

[Ethan Goode]

@ Walter Roberson

Could you provide a little more commentary?

Here is a valid way to solve my problem:

%I write a function s.t.

function F= ekmodelv2(x)

Ti= 1;
Li= 1;
dni = 2;
theta = 4;
gam = 1;
beta = 0.5;

%p(1:4) = x(1:4) These are the four price equations.
%x(1:16) = x(5:20) These are the sixteen trade share equations.
%w(1:4) = x(21:24) These are the four wage equations.

F(1) = (gam.*(((Ti.*(dni.*(x(21).^(beta))*(x(1).^(beta))).^(-theta)))...
    + ((Ti.*(dni.*(x(22).^(beta))*(x(2).^(beta))).^(-theta)))...
    + ((Ti.*(dni.*(x(23).^(beta))*(x(3).^(beta))).^(-theta)))...
    + ((Ti.*(dni.*(x(24).^(beta))*(x(4).^(beta))).^(-theta)))))-x(1);
F(2) = (gam.*(((Ti.*(dni.*(x(21).^(beta))*(x(1).^(beta))).^(-theta)))...
    + ((Ti.*(dni.*(x(22).^(beta))*(x(2).^(beta))).^(-theta)))...
    + ((Ti.*(dni.*(x(23).^(beta))*(x(3).^(beta))).^(-theta)))...
    + ((Ti.*(dni.*(x(24).^(beta))*(x(4).^(beta))).^(-theta)))))-x(2);
F(3) = (gam.*(((Ti.*(dni.*(x(21).^(beta))*(x(1).^(beta))).^(-theta)))...
    + ((Ti.*(dni.*(x(22).^(beta))*(x(2).^(beta))).^(-theta)))...
    + ((Ti.*(dni.*(x(23).^(beta))*(x(3).^(beta))).^(-theta)))...
    + ((Ti.*(dni.*(x(24).^(beta))*(x(4).^(beta))).^(-theta)))))-x(3);
F(4) = (gam.*(((Ti.*(dni.*(x(21).^(beta))*(x(1).^(beta))).^(-theta)))...
    + ((Ti.*(dni.*(x(22).^(beta))*(x(2).^(beta))).^(-theta)))...
    + ((Ti.*(dni.*(x(23).^(beta))*(x(3).^(beta))).^(-theta)))...
    + ((Ti.*(dni.*(x(24).^(beta))*(x(4).^(beta))).^(-theta)))))-x(4);

F(5) = Ti.*(gam.*dni.*(x(21).^(beta))*(x(1).^(1-beta))*(1./x(1))).^(theta) - (x(5));
F(6) = Ti.*(gam.*dni.*(x(21).^(beta))*(x(1).^(1-beta))*(1./x(2))).^(theta) - (x(6));
F(7) = Ti.*(gam.*dni.*(x(21).^(beta))*(x(1).^(1-beta))*(1./x(3))).^(theta) - (x(7));
F(8) = Ti.*(gam.*dni.*(x(21).^(beta))*(x(1).^(1-beta))*(1./x(4))).^(theta) - (x(8));
F(9) = Ti.*(gam.*dni.*(x(22).^(beta))*(x(2).^(1-beta))*(1./x(1))).^(theta) - (x(9));
F(10) = Ti.*(gam.*dni.*(x(22).^(beta))*(x(2).^(1-beta))*(1./x(2))).^(theta) - (x(10));
F(11) = Ti.*(gam.*dni.*(x(22).^(beta))*(x(2).^(1-beta))*(1./x(3))).^(theta) - (x(11));
F(12) = Ti.*(gam.*dni.*(x(22).^(beta))*(x(2).^(1-beta))*(1./x(4))).^(theta) - (x(12));
F(13) = Ti.*(gam.*dni.*(x(23).^(beta))*(x(3).^(1-beta))*(1./x(1))).^(theta) - (x(13));
F(14) = Ti.*(gam.*dni.*(x(23).^(beta))*(x(3).^(1-beta))*(1./x(2))).^(theta) - (x(14));
F(15) = Ti.*(gam.*dni.*(x(23).^(beta))*(x(3).^(1-beta))*(1./x(3))).^(theta) - (x(15));
F(16) = Ti.*(gam.*dni.*(x(23).^(beta))*(x(3).^(1-beta))*(1./x(4))).^(theta) - (x(16));
F(17) = Ti.*(gam.*dni.*(x(24).^(beta))*(x(4).^(1-beta))*(1./x(1))).^(theta) - (x(17));
F(18) = Ti.*(gam.*dni.*(x(24).^(beta))*(x(4).^(1-beta))*(1./x(2))).^(theta) - (x(18));
F(19) = Ti.*(gam.*dni.*(x(24).^(beta))*(x(4).^(1-beta))*(1./x(3))).^(theta) - (x(19));
F(20) = Ti.*(gam.*dni.*(x(24).^(beta))*(x(4).^(1-beta))*(1./x(4)).^(theta) - (x(20)));

F(21) = ((x(5).*x(21).*Li)+(x(6).*x(22).*Li)+(x(7).*x(23).*Li)+(x(8).*x(24).*Li))...
    -(x(21).*Li);
F(22) = ((x(9).*x(21).*Li)+(x(10).*x(22).*Li)+(x(11).*x(23).*Li)+(x(12).*x(24).*Li))...
    -(x(22).*Li);
F(23) = ((x(13).*x(21).*Li)+(x(14).*x(22).*Li)+(x(15).*x(23).*Li)+(x(16).*x(24).*Li))...
    -(x(23).*Li);
F(24) = ((x(17).*x(21).*Li)+(x(18).*x(22).*Li)+(x(19).*x(23).*Li)+(x(20).*x(24).*Li))...
    -(x(24).*Li);
end
%And then,
x0 = ones(1, 24);
[x, fval] = fsolve(@ekmodelv2, x0);

This solves 24 equations for 24 unknowns. This code, in my opinion, is like opening a paint can with a battle axe. I typed out every single equation. I can build a cell array of function handles, but it is unclear to me how to then use fsolve the system of non linear equations. I do not have a lot of experience using fsolve in MATLAB. Previously, I have only ever needed to use it to solve three equations at most. How should I go about generalizing this to solve n equations in your opinion?

[Walter Roberson]

F(1) = (gam.*(((Ti.*(dni.*(x(21).^(beta))*(x(1).^(beta))).^(-theta)))...
    + ((Ti.*(dni.*(x(22).^(beta))*(x(2).^(beta))).^(-theta)))...
    + ((Ti.*(dni.*(x(23).^(beta))*(x(3).^(beta))).^(-theta)))...
    + ((Ti.*(dni.*(x(24).^(beta))*(x(4).^(beta))).^(-theta)))))-x(1);
F(2) = (gam.*(((Ti.*(dni.*(x(21).^(beta))*(x(1).^(beta))).^(-theta)))...
    + ((Ti.*(dni.*(x(22).^(beta))*(x(2).^(beta))).^(-theta)))...
    + ((Ti.*(dni.*(x(23).^(beta))*(x(3).^(beta))).^(-theta)))...
    + ((Ti.*(dni.*(x(24).^(beta))*(x(4).^(beta))).^(-theta)))))-x(2);
F(3) = (gam.*(((Ti.*(dni.*(x(21).^(beta))*(x(1).^(beta))).^(-theta)))...
    + ((Ti.*(dni.*(x(22).^(beta))*(x(2).^(beta))).^(-theta)))...
    + ((Ti.*(dni.*(x(23).^(beta))*(x(3).^(beta))).^(-theta)))...
    + ((Ti.*(dni.*(x(24).^(beta))*(x(4).^(beta))).^(-theta)))))-x(3);
F(4) = (gam.*(((Ti.*(dni.*(x(21).^(beta))*(x(1).^(beta))).^(-theta)))...
    + ((Ti.*(dni.*(x(22).^(beta))*(x(2).^(beta))).^(-theta)))...
    + ((Ti.*(dni.*(x(23).^(beta))*(x(3).^(beta))).^(-theta)))...
    + ((Ti.*(dni.*(x(24).^(beta))*(x(4).^(beta))).^(-theta)))))-x(4);

The only difference I can see between those four is what is being subtracted at the very end. I re-checked several times, but that is the only difference I can find. If I have not missed anything, then the only way those can be true is if x(1) and x(2) and x(3) and x(4) are all equal to each other. If that is what is desired, then just force them to be equal by substitution.

[Walter Roberson]

Example:

K3 = 5;
for K1 = 21:24
    for K2 = 1:4
        F{K3} = @(x) Ti.*(gam.*dni.*(x(K1).^(beta))*(x(K2).^(1-beta))*(1./x(K2))).^(theta) - (x(K3));
        K3 = K3 + 1;
    end
end

for K = 1 : 4
  F{20+K} = @(x) ((x(4*K+1).*x(21).*Li)+(x(4*K+2).*x(22).*Li)+(x(4*K+3).*x(23).*Li)+(x(4*K+4).*x(24).*Li))...
    -(x(20+K).*Li);
end

solution = fsolve(@(x) arrayfun(@(f)f(x), F), x0);

[Ethan Goode]

This is a first attempt at solving the model with all of the interesting variables constant across the n countries. This is why the only difference between many equations is what is being subtracted at the very end. In the end, I will want to make dni, Ti, and even Li vary. So, I don't want to let all the x's equal each other (even though letting matlab solve the system gives you exactly that) in the code.

Here is what I have tried

function F = price_equations(x)

%%Put parameters here
Ti= 1;
Li= 1;
dni = 2; 
theta = 4;
gam = 1;
beta = 0.5;
n = 4;
for i = 1:n
    for j = 1:n
        for k = n+1:n+n^2
        F(i,1) = {@(x) (gam.*(sum(Ti.*(dni.*((x(k).^(beta).*(x(j).^(1-beta)))).^(-theta)))))-x(i)};  
        end
    end
end

My understanding is that the above command will give me a nx1 cell with each cell containing a function as defined above.

Similarly I repreat this in seperate files for the other types of functions. I'll put them below just to keep everything well documented.

function F = wage_equations(x)
%%Put parameters here
Ti= 1;
Li= 1;
dni = 2; 
theta = 4;
gam = 1;
beta = 0.5;
n = 4;
for i = (n + n^2+1) : (2*n + n^2)
    for k =(n+n^2+1) : (2*n + n^2)
        for j = n+1:n + n^2
     F(i, 1) = {@(x)(sum(x(j).*x(k).*Li))-(x(i).*Li)};
        end
    end
end

And,

function F = tradeshare_equations(x)
%%Put parameters here
Ti= 1;
Li= 1;
dni = 2; 
theta = 4;
gam = 1;
beta = 0.5;
n = 4;
for i = n+1:n + n^2
    for j = (n+n^2+1) : (2*n + n^2)
        for k = 1:n
    F(i, 1)= {@(x) Ti.*(gam.*dni.*(x(j).^(beta))*(x(k).^(1-beta))*(1./x(l))).^(theta) - (x(i))};
        end
    end
end

[Ethan Goode]

After defining the three types of functions I should be able to create a cell array of function handles:

handles = {@price_equations, @wage_equations, @tradeshare_equations};

Which gives me a 1x3 cell array.

I should then be able to define a guess x0 and then,

sol = fsolve(@(x) arrayfun(@(F)F(x), handles), x0);

Is my understanding of your initial comment incorrect? I have never worked with a cell array of function handles before. Ideally I want fsolve to output a nx(n+2) matrix of solutions. In the case where n = 4, I want to get a 4x6 matrix, where column one is the price equation solutions, column two is the wage equation solutions, and the rest of the columns are the trade share equations.

How do I define my guess?

In my previous code I can just have a row vector with 2n + n^2 elements. For four countries this is a 1x24. When fsolve solves an array function how is it expecting the intial guess?

Is there anything "off" with the code I provide? I am very new to using fsolve and the cell data type.

Thanks,

Ethan

[Walter Roberson]

F(i,1) = {@(x) (gam.*(sum(Ti.*(dni.*((x(k).^(beta).*(x(j).^(1-beta)))).^(-theta)))))-x(i)};  

That code is overwriting F(i,1) for each for j for k iteration, so the result will be n cell array in which only the last j and last k value "count". The wage equations have the same problem.

You could consider having the code generate a 3D cell array of only that kind of equation, and then reshape it into a vector at the end, and then

handles = [price_equations; wage_equations; tradeshare_equations];
sol = fsolve(@() arrayfun(@(F)F(x), handles), x0);

[Ethan Goode]

Thank you Walter,

Could you provide an example as to how I might generate a 3D cell array of only that kind of equation, and then reshape it? Or maybe link to an example that could help me make some progress?

Is there no way to use the for loop to create the 4x1 cell array of functions that I want in my earlier comment for wage_equations (as an example)?

Also, what does an initial guess look like for using fsolve on an to solve an array function?

[Walter Roberson]

function F = price_equations()

%%Put parameters here
Ti= 1;
Li= 1;
dni = 2; 
theta = 4;
gam = 1;
beta = 0.5;
n = 4;
F = cell(n, n, n^2);
for i = 1:n
    for j = 1:n
        for k = 1 : n^2
          F{i,j,k} = @(x) (gam.*(sum(Ti.*(dni.*((x(k+n).^(beta).*(x(j).^(1-beta)))).^(-theta)))))-x(i);
        end
    end
end

F = F(:);

Notice by the way that you do not pass x to this, as you are dynamically generating functions that will have some x passed to them.

[Walter Roberson]

function F = wage_equations()
%%Put parameters here
Ti= 1;
Li= 1;
dni = 2; 
theta = 4;
gam = 1;
beta = 0.5;
n = 4;
F = cell(n, n, n);
ibase = n + n^2;
kbase = n + n^2;
jbase = n;
for i = 1 : n
    for k = 1 : n
        for j = 1 : n
          F{i, j, k} = @(x)(sum(x(j+jbase).*x(k+kbase).*Li))-(x(i+ibase).*Li);
        end
    end
end
F = F(:);

This is not 4 (n) equations, this is n^3 equations.. provided that you coded correctly earlier.

[Celestine Siameh]

Please was this question answered, as I am solving a similar code but in my case is multi country and multi sector, armington model. The above as guided me but I want to know if the final code is correct. Thanks.

[Walter Roberson]

It looks to me as if perhaps the poster deleted at least one of their comments, so I am not sure whether the code ended up working for them.

[Celestine Siameh]

Thanks Walter. Let me check with him then. Hello Ethan Goode, did the code worked for you.