Effacer les filtres
Effacer les filtres

Loop over initial guesses with fsolve

3 vues (au cours des 30 derniers jours)
Olimpia
Olimpia le 29 Jan 2023
Commenté : Davide Masiello le 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.
Any idea? Thanks in advance!
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 commentaires
Alex Sha
Alex Sha le 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
Olimpia le 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...

Connectez-vous pour commenter.

Réponse acceptée

Davide Masiello
Davide Masiello le 29 Jan 2023
Modifié(e) : Davide Masiello le 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 commentaires
Olimpia
Olimpia le 30 Jan 2023
Thanks a lot!
Davide Masiello
Davide Masiello le 30 Jan 2023
No problem. If it helped, please consider to Accept the answer!

Connectez-vous pour commenter.

Plus de réponses (0)

Catégories

En savoir plus sur Loops and Conditional Statements dans Help Center et File Exchange

Tags

Produits

Community Treasure Hunt

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

Start Hunting!

Translated by