Solving N non-linear equations using fsolve. How do I pass these equations into my function without typing them out individually?
3 次查看(过去 30 天)
显示 更早的评论
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?
15 个评论
Ethan Goode
2020-6-30
编辑:Ethan Goode
2020-6-30
@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
2020-6-30
Build a cell array of function handles.
fsolve(@(x) arrayfun(@(F)F(x), CellOfHandles), x0)
Ethan Goode
2020-6-30
@ 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
2020-6-30
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
2020-6-30
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
2020-7-1
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
2020-7-1
编辑:Ethan Goode
2020-7-1
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
2020-7-1
编辑:Walter Roberson
2020-7-1
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
2020-7-1
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
2020-7-1
编辑:Walter Roberson
2020-7-1
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
2020-7-1
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
2021-8-21
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
2021-8-21
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
2021-8-21
Thanks Walter. Let me check with him then. Hello Ethan Goode, did the code worked for you.
回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Surrogate Optimization 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!发生错误
由于页面发生更改,无法完成操作。请重新加载页面以查看其更新后的状态。
您也可以从以下列表中选择网站:
如何获得最佳网站性能
选择中国网站(中文或英文)以获得最佳网站性能。其他 MathWorks 国家/地区网站并未针对您所在位置的访问进行优化。
美洲
- América Latina (Español)
- Canada (English)
- United States (English)
欧洲
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom(English)
亚太
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)