Fsolve initial condition accuracy

how to solve this ?

i am trying to solve the following four equations using Fsolve and the problem is that when i use zeros for initial guess the result is wrong

X = fsolve(myfun, [0;0;0;0])

X =

-0.0000

-0.0000

0.0137

0.0103

but when i use an initial near too the solution the result is corrct

X = fsolve(myfun, [0.009095;0.0005304;0;0])

X =

0.0090

0.0005

0.0229

0.0107

[Walter Roberson]

You do not appear to have posted the four equations?

[Ahmed hanafy]

myfun = @(x)[
    x(1)*x(3)*real(1-m)-x(2)*x(4)*real(1-m)-imag(1-m)*x(2)*x(3)-imag(1-m)*x(1)*x(4)+x(1)*real(zab(2,2)*(1-m)-m*zab(1,1)-h)-x(2)*imag(zab(2,2)*(1-m)-m*zab(1,1)-h)-x(3)*real(h+zab(1,1))+x(4)*imag(h+zab(1,1))-real(m)*x(1)*x(1)+real(m)*x(2)*x(2)+imag(m)*2*x(1)*x(2)-real(zab(1,1)*zab(2,2)+h*zab(1,1)+h*zab(2,2));
    x(1)*x(3)*imag(1-m)-x(2)*x(4)*imag(1-m)+real(1-m)*x(2)*x(3)+real(1-m)*x(1)*x(4)+x(1)*imag(zab(2,2)*(1-m)-m*zab(1,1)-h)+x(2)*real(zab(2,2)*(1-m)-m*zab(1,1)-h)-x(3)*imag(h+zab(1,1))-x(4)*real(h+zab(1,1))-imag(m)*x(1)*x(1)+imag(m)*x(2)*x(2)-real(m)*2*x(1)*x(2)-imag(zab(1,1)*zab(2,2)+h*zab(1,1)+h*zab(2,2));
    x(1)*x(3)*real(1-m)-x(2)*x(4)*real(1-m)+imag(m)*x(2)*x(3)+x(1)*x(4)*imag(m)+x(1)*real(zab(2,2)-c-m*(zab(1,1)+zab(2,2)))+x(3)*(real(zab(1,1)-c))-x(2)*imag(zab(2,2)-c-m*(zab(1,1)+zab(2,2)))-x(4)*imag(zab(1,1)-c)-real(m)*x(1)*x(1)+real(m)*x(2)*x(2)+2*x(1)*x(2)*imag(m)+real(zab(1,1)*zab(2,2)-c*zab(1,1)-c*zab(2,2));
    -1*x(1)*x(3)*imag(m)+x(2)*x(4)*imag(m)+real(1-m)*x(2)*x(3)+x(1)*x(4)*real(1-m)+x(1)*imag(zab(2,2)-c-m*(zab(1,1)+zab(2,2)))+x(3)*(imag(zab(1,1)-c))+x(2)*real(zab(2,2)-c-m*(zab(1,1)+zab(2,2)))+x(4)*real(zab(1,1)-c)-imag(m)*x(1)*x(1)+imag(m)*x(2)*x(2)-real(m)*x(1)*x(2)*2+imag(zab(1,1)*zab(2,2)-c*zab(1,1)-c*zab(2,2));
    ]
m = 0.5049 - 0.0057i
c = 0.0183 + 0.2868i
h =-0.0183 - 0.2868i

[Walter Roberson]

The first set of results, that you said are wrong, are in fact correct results. The second set that you said are correct, are not true roots and not especially close to true roots, and would not appear if you gave fsolve a stricter tolerance.

There are four roots to the equations, two of which involve imaginary components. x(4) for the second of the two is about -0.568

Answer 1

[Star Strider]

Nonlinear solvers of all sorts (fsolve, lsqcurvefit, etc.) are very sensitive to the initial estimates, since they will search for the solution nearest the initial parameter estimates. That is what you are seeing.

As a general rule, it is best not to use any zero values in your initial parameter estimates, largely because in some instances those can lead to NaN or Inf initial results, and the solver will stop.

[Ahmed hanafy]

thank you for replying

i understand your anser but i need the initial solution to be general .

is there any thing to do ?

[Star Strider]

My pleasure.

I am not certain what you intend by ‘general’. If you want all the solutions of your system, (since having four parameters you cannot easily plot it to see them), you simply have to search for them. This likely requires either a loop with random starting values, then using the uniquetol function with the 'rows' option to eliminate multiple similar results, or one of the Global Optimization Toolbox functions that will search the entire parameter space for solutions.

[Ahmed hanafy]

myfun = @(x)[
    x(1)*x(3)*real(1-m)-x(2)*x(4)*real(1-m)-imag(1-m)*x(2)*x(3)-imag(1-m)*x(1)*x(4)+x(1)*real(zab(2,2)*(1-m)-m*zab(1,1)-h)-x(2)*imag(zab(2,2)*(1-m)-m*zab(1,1)-h)-x(3)*real(h+zab(1,1))+x(4)*imag(h+zab(1,1))-real(m)*x(1)*x(1)+real(m)*x(2)*x(2)+imag(m)*2*x(1)*x(2)-real(zab(1,1)*zab(2,2)+h*zab(1,1)+h*zab(2,2));
    x(1)*x(3)*imag(1-m)-x(2)*x(4)*imag(1-m)+real(1-m)*x(2)*x(3)+real(1-m)*x(1)*x(4)+x(1)*imag(zab(2,2)*(1-m)-m*zab(1,1)-h)+x(2)*real(zab(2,2)*(1-m)-m*zab(1,1)-h)-x(3)*imag(h+zab(1,1))-x(4)*real(h+zab(1,1))-imag(m)*x(1)*x(1)+imag(m)*x(2)*x(2)-real(m)*2*x(1)*x(2)-imag(zab(1,1)*zab(2,2)+h*zab(1,1)+h*zab(2,2));
    x(1)*x(3)*real(1-m)-x(2)*x(4)*real(1-m)+imag(m)*x(2)*x(3)+x(1)*x(4)*imag(m)+x(1)*real(zab(2,2)-c-m*(zab(1,1)+zab(2,2)))+x(3)*(real(zab(1,1)-c))-x(2)*imag(zab(2,2)-c-m*(zab(1,1)+zab(2,2)))-x(4)*imag(zab(1,1)-c)-real(m)*x(1)*x(1)+real(m)*x(2)*x(2)+2*x(1)*x(2)*imag(m)+real(zab(1,1)*zab(2,2)-c*zab(1,1)-c*zab(2,2));
    -1*x(1)*x(3)*imag(m)+x(2)*x(4)*imag(m)+real(1-m)*x(2)*x(3)+x(1)*x(4)*real(1-m)+x(1)*imag(zab(2,2)-c-m*(zab(1,1)+zab(2,2)))+x(3)*(imag(zab(1,1)-c))+x(2)*real(zab(2,2)-c-m*(zab(1,1)+zab(2,2)))+x(4)*real(zab(1,1)-c)-imag(m)*x(1)*x(1)+imag(m)*x(2)*x(2)-real(m)*x(1)*x(2)*2+imag(zab(1,1)*zab(2,2)-c*zab(1,1)-c*zab(2,2));
    ]
m =

   0.5049 - 0.0057i
   c =

   0.0183 + 0.2868i
   h =

  -0.0183 - 0.2868i

[Star Strider]

What is ‘zab’?

[Ahmed hanafy]

it a constant

zab(1,1)=zab(2,2)=0.0299 + 0.5684i

[Star Strider]

You appear to have a very large number of acceptable solutions.

If you have the Global Optimization Toolbox, run this (change the ‘k2’ limit to get more or fewer initial estimates):

for k2 = 1:25
m = 0.5049 - 0.0057i;
c = 0.0183 + 0.2868i;
h = -0.0183 - 0.2868i;
zab = eye(2)*(0.0299 + 0.5684i);
myfun = @(x)[x(1)*x(3)*real(1-m)-x(2)*x(4)*real(1-m)-imag(1-m)*x(2)*x(3)-imag(1-m)*x(1)*x(4)+x(1)*real(zab(2,2)*(1-m)-m*zab(1,1)-h)-x(2)*imag(zab(2,2)*(1-m)-m*zab(1,1)-h)-x(3)*real(h+zab(1,1))+x(4)*imag(h+zab(1,1))-real(m)*x(1)*x(1)+real(m)*x(2)*x(2)+imag(m)*2*x(1)*x(2)-real(zab(1,1)*zab(2,2)+h*zab(1,1)+h*zab(2,2));
    x(1)*x(3)*imag(1-m)-x(2)*x(4)*imag(1-m)+real(1-m)*x(2)*x(3)+real(1-m)*x(1)*x(4)+x(1)*imag(zab(2,2)*(1-m)-m*zab(1,1)-h)+x(2)*real(zab(2,2)*(1-m)-m*zab(1,1)-h)-x(3)*imag(h+zab(1,1))-x(4)*real(h+zab(1,1))-imag(m)*x(1)*x(1)+imag(m)*x(2)*x(2)-real(m)*2*x(1)*x(2)-imag(zab(1,1)*zab(2,2)+h*zab(1,1)+h*zab(2,2));
    x(1)*x(3)*real(1-m)-x(2)*x(4)*real(1-m)+imag(m)*x(2)*x(3)+x(1)*x(4)*imag(m)+x(1)*real(zab(2,2)-c-m*(zab(1,1)+zab(2,2)))+x(3)*(real(zab(1,1)-c))-x(2)*imag(zab(2,2)-c-m*(zab(1,1)+zab(2,2)))-x(4)*imag(zab(1,1)-c)-real(m)*x(1)*x(1)+real(m)*x(2)*x(2)+2*x(1)*x(2)*imag(m)+real(zab(1,1)*zab(2,2)-c*zab(1,1)-c*zab(2,2));
    -1*x(1)*x(3)*imag(m)+x(2)*x(4)*imag(m)+real(1-m)*x(2)*x(3)+x(1)*x(4)*real(1-m)+x(1)*imag(zab(2,2)-c-m*(zab(1,1)+zab(2,2)))+x(3)*(imag(zab(1,1)-c))+x(2)*real(zab(2,2)-c-m*(zab(1,1)+zab(2,2)))+x(4)*real(zab(1,1)-c)-imag(m)*x(1)*x(1)+imag(m)*x(2)*x(2)-real(m)*x(1)*x(2)*2+imag(zab(1,1)*zab(2,2)-c*zab(1,1)-c*zab(2,2));];


ftns = @(theta) norm(myfun(theta));
PopSz = 500;
Parms = 4;
opts = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-3, 'MaxGenerations',2E3, 'PlotFcn','gaplotbestf');
t0 = clock;
fprintf('\nStart Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t0)
% [B,fval,exitflag,output,population,scores] = ga(ftns, 53, [],[],[],[],zeros(53,1),Inf(53,1),[],[],opts)
% [theta,fval,exitflag,output] = ga(ftns, Parms, [],[],[],[],zeros(Parms,1),Inf(Parms,1),[],[],opts)
[theta,fval,exitflag,output] = ga(ftns, Parms, [],[],[],[],[],[],[],[],opts)
t1 = clock;
fprintf('\nStop Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t1)
GA_Time = etime(t1,t0)
fprintf('\nElapsed Time: %23.15E\t\t%02d:%02d:%02d.%03d\n', GA_Time, hmsdv(GA_Time))

fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(theta)
    fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1))
end

[X,fsfval] = fsolve(myfun, theta)
OutMtx(k2,:) = X;
end
Xu = uniquetol(OutMtx, 0.5, 'ByRows',1)

and using uniquetol (to eliminate similar results) with the last loop I ran:

Xu =
   -0.0869   -0.0805   -0.0768   -0.0698
   -0.0296   -0.0970   -0.0187   -0.0879
   -0.0267   -0.5729   -0.0267   -0.5730
    0.0196   -0.0219    0.0331   -0.0124
    0.0680   -0.0497    0.0818   -0.0419
     

and fsolve converged successfully on all of these. You will probably find even more with a more loop iterations.

[Walter Roberson]

After converting the floating point coefficients to rational, using a symbolic solver you can find the four exact solutions:

x1 = 0, x2 = 0, x3 = 30518438974403069375/2235843560205474332672, x4 = 185072313798727454375/17886748481643794661376

x1 = -2990894712572022840244004951736921950/100029923498226993240745359851966741761, x2 = -56857008516399963015166917259264175100/100029923498226993240745359851966741761, x3 = -2033176198195073422954255136858330487554884173791875/67999203953587065427434178402871152642257991411695616, x4 = -38650747527224644283355248355941704068570010837901875/67999203953587065427434178402871152642257991411695616

x1 = -1495447356286011420122002475868460975/100029923498226993240745359851966741761 + 28428504258199981507583458629632087550*1i/100029923498226993240745359851966741761, x2 = -28428504258199981507583458629632087550/100029923498226993240745359851966741761 - 1495447356286011420122002475868460975*1i/100029923498226993240745359851966741761, x3 = -548591028899807326465951446601601742778468623538171351875/67517361594371947481815379814708407654634951284552340865024 + 9768865850400612492690517613022042225167452498265998546875*1i/33758680797185973740907689907354203827317475642276170432512, x4 = -588723019670172045893198682501481771975163085767890788125/2109917549824123358806730619209637739207342227642260652032 - 1470178082754855806213249838356950614941601640998214421875*1i/67517361594371947481815379814708407654634951284552340865024

The second of those has a duplicate copy for a total of 4.

[Star Strider]

@Walter — Thank you.