Loop over initial guesses with fsolve

5 Ansichten (letzte 30 Tage)
Olimpia am 29 Jan. 2023
Kommentiert: Davide Masiello am 30 Jan. 2023
Hi everyone!
I have a system of 7 non-linear equations with both real solutions and complex solutions, depending on my initial guess. In order to find the real solutions, I want to loop over initial guesses.
I have 7 variables, I want my initial guess vector to take the values [1, 1, 1, 1, 0.1, 0.1, 0.1], my second initial guess vector to increment over the first four guesses by 0.1 and over the last 3 by 0.01 so that it takes the values [1.1, 1.1, 1.1, 1.1, 0.11, 0.11, 0.11, 0.11] until the guess [1.9, 1.9, 1.9, 1.9, 0.1, 0.1, 0.1] after ten steps.
Finally, I want to save all the guesses in a matrix 10x7 (or 7x10).
I tried the following code, and I have the error
Index exceeds the number of array elements (1).
Caused by:
Failure in initial objective function evaluation. FSOLVE cannot continue.
func=@f;
x0=[1:0.1:1.9,1:0.1:1.9,1:0.1:1.9,1:0.1:1.9,0.01:0.01:0.1,0.01:0.01:0.1,0.01:0.01:0.1]
x = zeros(10,7)
for i=1:10
x(i)=fsolve(func,x0(i))
end
function my_func=f(x)
rho = 1/2;
c1=1;
c2=11/10;
c3=12/10;
c4=7/10;
y = 20;
pall = 2;
sigma = 10;
eta = 2;
alpha=3/9;
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
2 KommentareKeine anzeigenKeine ausblenden
Alex Sha am 30 Jan. 2023
The unique real solution should be:
x1: 1.13597072587292
x2: 1.23385304815548
x3: 1.33911918335085
x4: 0.783110223199337
x5: 0.017879890424952
x6: 0.000114204714425553
x7: 1.05898294440764E-6
Olimpia am 30 Jan. 2023
Thank you @Alex Sha, but how did you obtain the real solution? When I tried the loop over the guesses I still got only complex solutions...

Melden Sie sich an, um zu kommentieren.

Akzeptierte Antwort

Davide Masiello am 29 Jan. 2023
Bearbeitet: Davide Masiello am 29 Jan. 2023
Try with
func=@f;
x0=[1:0.1:1.9;1:0.1:1.9;1:0.1:1.9;1:0.1:1.9;0.01:0.01:0.1;0.01:0.01:0.1;0.01:0.01:0.1];
x = zeros(10,7);
for i=1:10
x(i,:)=fsolve(func,x0(:,i))
end
Please note the use of ";" instead of "," in the x0 matrix and the correct recursive indexing in the loop (x0(:,i)).
4 Kommentare2 ältere Kommentare anzeigen2 ältere Kommentare ausblenden
Olimpia am 30 Jan. 2023
Thanks a lot!
Davide Masiello am 30 Jan. 2023

Melden Sie sich an, um zu kommentieren.

Kategorien

Mehr zu Loops and Conditional Statements finden Sie in Help Center und File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by