# Obtaining a matrix of solutions with fsolve

6 views (last 30 days)
Olimpia on 31 Jan 2023
Edited: Torsten on 31 Jan 2023
Hello everyone !
I am trying to solve a system of 7 non-linear equations and 7 variables.
I have 4 parameters (c1, c2, c3, c4) and I want to solve my system for different (specific) values of these four parameters. Concretely, I am trying to solve the seven equations over my 210 possible vectors (c1i, c2i, c3i, c4i), that I upload from the file « grid_wide_4firms ». Instead of obtaining a vector of solution (1x7), I then wish to obtain a matrix (210x7).
I have therefore created a loop outside that replaces the vector c at each step, and store the results each step. But have the error:
Function definition are not supported in this context. Functions can only be created as local or nested functions in code files
Any idea ?
PS : I attach the file grid_wide_4firms with the different values of c1, c2, c3, c4 and the code that works over one specific vector of ci (so that hopefully it is clearer)
final_matrix=zeros(210,7);
for i=1:210
clearvars c1 c2 c3 c4 x my_func
func=@f;
x0=[1.99;1.99;1.99;1.99;0.0109;0.0109;0.0109]
x = zeros(1,7);
x=fsolve(func,x0),
display(x)
final_matrix((i),1)=x
function my_func=f(x)
% import the data
opts = delimitedTextImportOptions("NumVariables", 4);
opts.DataLines = [2, Inf];
opts.Delimiter = ",";
opts.VariableNames = ["cost1", "cost2", "cost3", "cost4"];
opts.VariableTypes = ["double", "double", "double", "double"];
opts.ExtraColumnsRule = "ignore";
costs= table2array(gridwide4firms2);
% gen c1, c2, c3, c4 so they are vectors
c1=costs((i),1)
c2=costs((i),2)
c3=costs((i),3)
c4=costs((i),4)
% parameters
rho = 1/2;
y = 20;
pall = 2;
sigma = 10;
eta = 2;
alpha=3/4;
my_func(1)= (sigma*(1 - (x(1)/(x(1)^(1 - sigma) + x(2)^(1 - sigma) + x(3)^(1 - sigma) + (x(4)*(1 + (x(5)^rho + x(6)^rho + x(7)^rho)^(alpha/rho)))^(1 - sigma))^(1/(1 - sigma)))^(1 - sigma)) + eta*(x(1)/(x(1)^(1 - sigma) + x(2)^(1 - sigma) + x(3)^(1 - sigma) + (x(4)*(1 + (x(5)^rho + x(6)^rho + x(7)^rho)^(alpha/rho)))^(1 - sigma))^(1/(1 - sigma)))^(1 - sigma))/(sigma*(1 - (x(1)/(x(1)^(1 - sigma) + x(2)^(1 - sigma) + x(3)^(1 - sigma) + (x(4)*(1 + (x(5)^rho + x(6)^rho + x(7)^rho)^(alpha/rho)))^(1 - sigma))^(1/(1 - sigma)))^(1 - sigma)) + eta*(x(1)/(x(1)^(1 - sigma) + x(2)^(1 - sigma) + x(3)^(1 - sigma) + (x(4)*(1 + (x(5)^rho + x(6)^rho + x(7)^rho)^(alpha/rho)))^(1 - sigma))^(1/(1 - sigma)))^(1 - sigma) - 1)*c1 -x(1);
my_func(2)= (sigma*(1 - (x(2)/(x(1)^(1 - sigma) + x(2)^(1 - sigma) + x(3)^(1 - sigma) + (x(4)*(1 + (x(5)^rho + x(6)^rho + x(7)^rho)^(alpha/rho)))^(1 - sigma))^(1/(1 - sigma)))^(1 - sigma)) + eta*(x(2)/(x(1)^(1 - sigma) + x(2)^(1 - sigma) + x(3)^(1 - sigma) + (x(4)*(1 + (x(5)^rho + x(6)^rho + x(7)^rho)^(alpha/rho)))^(1 - sigma))^(1/(1 - sigma)))^(1 - sigma))/(sigma*(1 - (x(2)/(x(1)^(1 - sigma) + x(2)^(1 - sigma) + x(3)^(1 - sigma) + (x(4)*(1 + (x(5)^rho + x(6)^rho + x(7)^rho)^(alpha/rho)))^(1 - sigma))^(1/(1 - sigma)))^(1 - sigma)) + eta*(x(2)/(x(1)^(1 - sigma) + x(2)^(1 - sigma) + x(3)^(1 - sigma) + (x(4)*(1 + (x(5)^rho + x(6)^rho + x(7)^rho)^(alpha/rho)))^(1 - sigma))^(1/(1 - sigma)))^(1 - sigma) - 1)*c2 -x(2);
my_func(3)= (sigma*(1 - (x(3)/(x(1)^(1 - sigma) + x(2)^(1 - sigma) + x(3)^(1 - sigma) + (x(4)*(1 + (x(5)^rho + x(6)^rho + x(7)^rho)^(alpha/rho)))^(1 - sigma))^(1/(1 - sigma)))^(1 - sigma)) + eta*(x(3)/(x(1)^(1 - sigma) + x(2)^(1 - sigma) + x(3)^(1 - sigma) + (x(4)*(1 + (x(5)^rho + x(6)^rho + x(7)^rho)^(alpha/rho)))^(1 - sigma))^(1/(1 - sigma)))^(1 - sigma))/(sigma*(1 - (x(3)/(x(1)^(1 - sigma) + x(2)^(1 - sigma) + x(3)^(1 - sigma) + (x(4)*(1 + (x(5)^rho + x(6)^rho + x(7)^rho)^(alpha/rho)))^(1 - sigma))^(1/(1 - sigma)))^(1 - sigma)) + eta*(x(3)/(x(1)^(1 - sigma) + x(2)^(1 - sigma) + x(3)^(1 - sigma) + (x(4)*(1 + (x(5)^rho + x(6)^rho + x(7)^rho)^(alpha/rho)))^(1 - sigma))^(1/(1 - sigma)))^(1 - sigma) - 1)*c3 -x(3);
my_func(4)= (sigma*(1 - (x(4)*(1+(x(5)^rho + x(6)^rho + x(7)^rho)^(alpha/rho))/(x(1)^(1 - sigma) + x(2)^(1 - sigma) + x(3)^(1 - sigma) + (x(4)*(1 + (x(5)^rho + x(6)^rho + x(7)^rho)^(alpha/rho)))^(1 - sigma))^(1/(1 - sigma)))^(1-sigma)*(1 + (x(5)^rho + x(6)^rho + x(7)^rho)^(alpha/rho))^(1 - sigma)) + eta*(x(4)*(1+(x(5)^rho + x(6)^rho + x(7)^rho)^(alpha/rho))/(x(1)^(1 - sigma) + x(2)^(1 - sigma) + x(3)^(1 - sigma) + (x(4)*(1 + (x(5)^rho + x(6)^rho + x(7)^rho)^(alpha/rho)))^(1 - sigma))^(1/(1 - sigma)))^(1-sigma)*(1 + (x(5)^rho + x(6)^rho + x(7)^rho)^(alpha/rho))^(1 - sigma))/(sigma*(1 - (x(4)*(1+(x(5)^rho + x(6)^rho + x(7)^rho)^(alpha/rho))/(x(1)^(1 - sigma) + x(2)^(1 - sigma) + x(3)^(1 - sigma) + (x(4)*(1 + (x(5)^rho + x(6)^rho + x(7)^rho)^(alpha/rho)))^(1 - sigma))^(1/(1 - sigma)))^(1-sigma)*(1 + (x(5)^rho + x(6)^rho + x(7)^rho)^(alpha/rho))^(1 - sigma)) + eta*(x(4)*(1+(x(5)^rho + x(6)^rho + x(7)^rho)^(alpha/rho))/(x(1)^(1 - sigma) + x(2)^(1 - sigma) + x(3)^(1 - sigma) + (x(4)*(1 + (x(5)^rho + x(6)^rho + x(7)^rho)^(alpha/rho)))^(1 - sigma))^(1/(1 - sigma)))^(1-sigma)*(1 + (x(5)^rho + x(6)^rho + x(7)^rho)^(alpha/rho))^(1 - sigma) - 1)*c4 -x(4);
my_func(5)= pall^(eta - 1)*(x(1)^(1 - sigma) + x(2)^(1 - sigma) + x(3)^(1 - sigma) + (x(4)*(1 + (x(5)^rho + x(6)^rho + x(7)^rho)^(alpha/rho)))^(1 - sigma))^(1/(1 - sigma))^(sigma - eta)*y*x(1)^(-sigma)*(x(1) - c1)*(sigma - eta)*(x(4)*(1+(x(5)^rho + x(6)^rho + x(7)^rho)^(alpha/rho))/(x(1)^(1 - sigma) + x(2)^(1 - sigma) + x(3)^(1 - sigma) + (x(4)*(1 + (x(5)^rho + x(6)^rho + x(7)^rho)^(alpha/rho)))^(1 - sigma))^(1/(1 - sigma)))^(1-sigma)*(x(5)^rho + x(6)^rho + x(7)^rho)^((1 - rho)/rho)*x(5)^(rho - 1) - (1 + (x(5)^rho + x(6)^rho + x(7)^rho)^(alpha/rho))^sigma/(x(5)^(rho - 1)*(x(5)^rho + x(6)^rho + x(7)^rho)^((rho - 1)/rho)*alpha*(x(5)^rho + x(6)^rho + x(7)^rho)^(1/rho))^(alpha-1);
my_func(6)= pall^(eta - 1)*(x(1)^(1 - sigma) + x(2)^(1 - sigma) + x(3)^(1 - sigma) + (x(4)*(1 + (x(5)^rho + x(6)^rho + x(7)^rho)^(alpha/rho)))^(1 - sigma))^(1/(1 - sigma))^(sigma - eta)*y*x(2)^(-sigma)*(x(2) - c2)*(sigma - eta)*(x(4)*(1+(x(5)^rho + x(6)^rho + x(7)^rho)^(alpha/rho))/(x(1)^(1 - sigma) + x(2)^(1 - sigma) + x(3)^(1 - sigma) + (x(4)*(1 + (x(5)^rho + x(6)^rho + x(7)^rho)^(alpha/rho)))^(1 - sigma))^(1/(1 - sigma)))^(1-sigma)*(x(5)^rho + x(6)^rho + x(7)^rho)^((1 - rho)/rho)*x(6)^(rho - 1) - (1 + (x(5)^rho + x(6)^rho + x(7)^rho)^(alpha/rho))^sigma/(x(6)^(rho - 1)*(x(5)^rho + x(6)^rho + x(7)^rho)^((rho - 1)/rho)*alpha*(x(5)^rho + x(6)^rho + x(7)^rho)^(1/rho))^(alpha-1);
my_func(7)= pall^(eta - 1)*(x(1)^(1 - sigma) + x(2)^(1 - sigma) + x(3)^(1 - sigma) + (x(4)*(1 + (x(5)^rho + x(6)^rho + x(7)^rho)^(alpha/rho)))^(1 - sigma))^(1/(1 - sigma))^(sigma - eta)*y*x(3)^(-sigma)*(x(3) - c3)*(sigma - eta)*(x(4)*(1+(x(5)^rho + x(6)^rho + x(7)^rho)^(alpha/rho))/(x(1)^(1 - sigma) + x(2)^(1 - sigma) + x(3)^(1 - sigma) + (x(4)*(1 + (x(5)^rho + x(6)^rho + x(7)^rho)^(alpha/rho)))^(1 - sigma))^(1/(1 - sigma)))^(1-sigma)*(x(5)^rho + x(6)^rho + x(7)^rho)^((1 - rho)/rho)*x(7)^(rho - 1) - (1 + (x(5)^rho + x(6)^rho + x(7)^rho)^(alpha/rho))^sigma/(x(7)^(rho - 1)*(x(5)^rho + x(6)^rho + x(7)^rho)^((rho - 1)/rho)*alpha*(x(5)^rho + x(6)^rho + x(7)^rho)^(1/rho))^(alpha-1);
end
end

Torsten on 31 Jan 2023
Edited: Torsten on 31 Jan 2023
Don't solve for all x values at once.
Call "fsolve" in a loop and solve for each pair of combinations for (c1,c2,c3,c4) separately.
And something like x(1;i) does not exist. Indexing a 2d array works as x(1,i).
Torsten on 31 Jan 2023
Edited: Torsten on 31 Jan 2023
There is no general rule that generates good initial guesses.
If your inputs for c don't change much from call to call, you can use the result of a call as initial guess for the next call. If this doesn't work: tough luck.

### Categories

Find more on Time Series Objects in Help Center and File Exchange

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!