Why does fsolve seem not iterate towards the solution?
Afficher commentaires plus anciens
Hello all,
I am trying to sovle a two non-linear equation system using fsolve and dogleg method. My objective function along with its jacobian is like this
function [F jacF]= objective(x)
F(:,1) = ((((x(:,2)./10).*k).*(x(:,1)./100)).^2).*(rZ - Rs) +(( Cmax .* ( x(:,1)./100 ) ).^2).*( w.^2.*(rZ - Rs) ) - (((x(:,2)./10).*k).*(x(:,1)./100));
F(:,2) = (x(:,2).*k).^2.*(iZ - w.*Ls) + (x(:,2).*k).^2.*x(:,1).*((w.*Ls)./200) + x(:,1).*((w.*Ls)/200).*(w.*Cmax).^2 + (w.*Cmax).^2 .*(iZ -(w.*Ls));
if nargout > 1 % need Jacobian
jacF = [- k - (k.^2.*x(:,2).*x(:,1).*(Rs - rZ))./50, - (k.^2.*x(:,2).^2.*(Rs - rZ))./100 - (Cmax.^2.*w.^2.*(Rs - rZ))./100;
2.*k.^2.*x(:,2).*(iZ - Ls.*w) + (k.^2.*Ls.*x(:,2).*w.*x(:,1))./100,(Ls.*Cmax.^2.*w.^3)./200 + (Ls.*k.^2.*x(:,2).^2.*w)./200];
end
end
Then my configuration for fsolve looks like this
options = optimoptions('fsolve','Display','iter-detailed','PlotFcn',@optimplotfirstorderopt);
% options.StepTolerance = 1e-13;
options.OptimalityTolerance = 1e-12;
options.FunctionTolerance = 6e-11;
options.MaxIterations = 100000;
options.MaxFunctionEvaluations = 400;%*400;
options.Algorithm = 'trust-region-dogleg';%'trust-region'%'levenberg-marquardt';%
% options.FiniteDifferenceType= 'central';
options.SpecifyObjectiveGradient = true;
fun= @objective;
x0 = [x1',x2'];
% Solve the function fun
[gwc,fval,exitflag,output,jacobianEval] =fsolve(fun,x0,options);
Being the values of the equations
Rs =
0.1640
Ls =
1.1000e-07
Cmax =
7.0000e-11
w =
1.7040e+08
rZ =
12.6518
iZ =
14.5273
K =
0.1007
x0 =
70.56 0.0759
My problem comes because I don't understand why fsolve seems not to iteratate over x(:,1) as i was expecting. I do know that the solution for the above system and parameters should be x1=58.8 and x2=0.0775.
In order to test the convergence of the method I am setting as initial guess x0 = [x1*(1+20/100) 0.0759] = [70.56 0.0759] ( 20 percent error in x1 and a higer value on x2), but the solution given by fsolve is the initial point, why is this? Am I doing something incorrect in my settings?
Thanks in advance
7 commentaires
Torsten
le 12 Sep 2019
What do you get for "res" if you insert the line
res = objective([58.8,0.0775])
before calling fsolve ?
Torsten
le 12 Sep 2019
If you divide the first equation by x1, then solve for x2 in terms of x1, insert the expression in equation 2 and multiply by the denominator, you arrive at a cubic polynomial in x2 that can be solved using "roots".
Hi Torsten,
I don't think I am coming to the same expression that you for x2 in terms of x1, so far I have simplified the first equation as you mentioned, calling the constants as A,B,C 
and I reach to the right equation by dividing first equation by x1
Do you want me to solve the cuadratic above like this and insert it in F2?
Torsten
le 12 Sep 2019
Ok, with some constants I get
a1*x1*x2^2 + a2*x1 + a3*x2 = 0 (1st equation after division by x1)
b1*x2^2 + b2*x1*x2^2 + b3*x1 + b4 = 0 (2nd equation)
Now you can either solve the first or the second equation for x1 and insert the expression in the other one.
In the end, you'll get a polynomial of degree 4 in x2 that can be solved using "roots".
Javier
le 12 Sep 2019
Torsten
le 12 Sep 2019
If you have a good starting guess ...
Réponse acceptée
Torsten
le 12 Sep 2019
0 votes
k = 0.1007;
rZ = 12.6518;
Rs = 0.164;
Cmax = 7.0e-11;
W = 1.704e8;
iZ = 14.5273;
Ls = 1.1e-7;
a1 = (k/1000)^2*(rZ-Rs);
a2 = (Cmax/100*W)^2*(rZ-Rs);
a3 = -k/1000;
b1 = k^2*(iZ-W*Ls);
b2 = k^2*W*Ls/200;
b3 = W*Ls/200*(Cmax*W)^2;
b4 = (Cmax*W)^2*(iZ-W*Ls);
A = (a2*b1-b2*a3+a1*b4)/(a1*b1);
B = (-a3*b3+a2*b4)/(a1*b1);
disk = A^2/4-B;
if disk >=0
x21squared = -A/2+sqrt(disk);
x22squared = -A/2-sqrt(disk);
end
solx1 = zeros(4,1);
solx2 = zeros(4,1);
iflag1 = 0;
if x21squared >= 0
iflag1 = 1;
solx2(1) = sqrt(x21squared);
solx2(2) = -sqrt(x21squared);
end
iflag2 = 0;
if x22squared >= 0
iflag2 = 1;
solx2(3) = sqrt(x22squared);
solx2(4) = -sqrt(x22squared);
end
solx1 = zeros(4,1);
if iflag1 == 1
solx1(1) = -a3/(a1*solx2(1)^2+a2);
solx1(2) = -a3/(a1*solx2(2)^2+a2);
end
if iflag2 == 1
solx1(3) = -a3/(a1*solx2(3)^2+a2);
solx1(4) = -a3/(a1*solx2(4)^2+a2);
end
if iflag1 == 1
solx1(1)
solx2(1)
solx1(2)
solx2(2)
a1*solx1(1)*solx2(1)^2+a2*solx1(1)+a3
b1*solx2(1)^2+b2*solx1(1)*solx2(1)^2+b3*solx1(1)+b4
a1*solx1(2)*solx2(2)^2+a2*solx1(2)+a3
b1*solx2(2)^2+b2*solx1(2)*solx2(2)^2+b3*solx1(2)+b4
end
if iflag2 == 1
solx1(3)
solx2(3)
solx1(4)
solx2(4)
a1*solx1(3)*solx2(3)^2+a2*solx1(3)+a3
b1*solx2(3)^2+b2*solx1(3)*solx2(3)^2+b3*solx1(3)+b4
a1*solx1(4)*solx2(4)^2+a2*solx1(4)+a3
b1*solx2(4)^2+b2*solx1(4)*solx2(4)^2+b3*solx1(4)+b4
end
5 commentaires
Javier
le 12 Sep 2019
Torsten
le 12 Sep 2019
I solved the two equations analytically. The results are stored in solx1(1:4), solx2(1:4).
If you follow my advice to build the polynomial of degree 4 in x2, you'll notice that it only has even powers of x2:
x2^4 + A*x2^2 + B = 0
where A and B are given as in my code.
This equation can be solved for x2^2 by the quadratic formula.
Torsten
le 13 Sep 2019
This is likely to break the bi-quadratic equation that you made, then should I go again for an interative process right?
No. Insert the expression from F1 for x1 in F2, multiply by the denominator, order according to powers of x2 and use MATLAB's "roots" to solve for x2. This will come out much more stable than using "fsolve".
Javier
le 13 Sep 2019
Hi Torsten, thanks again.
I have followed your advice, since F1 hasn't change, the relationship between x1 and x2 based on first equation is the same as you mentioned.
If I now replace the value of x1 into F2 as you suggests, I will get an polynom in order fith instead of fourth, which I will resolve with roots
Plus de réponses (0)
Catégories
En savoir plus sur Systems of Nonlinear Equations dans Centre d'aide et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
