Norm of step in fsolve

Hi,

I have yet another fsolve question. Norm of step in fsolve seems to be able to show me false roots (which i have far too many). Is it possible to adjust the code so that it will continue to keep searching for roots until norm of step is<eps. I have tried to explain the problem more in detail below;

myfun = @det_sine_strip_Epol; % function
format long
r=3.5 

options=optimset('Display','iter','MaxIter',3000,'MaxFunEvals',3000, 'TolFun', 1.0e-16, 'TolX',1.0e-16); 
fsolve(myfun,r,options)

returns this;

which is correct as

myfun(3.819451575438322 - 0.002164123209083i)= 1.784467423482892e-15 + 2.718518249093796e-16i

now if i change initial guess to 3.3, i have this;

which is not correct as

myfun(3.782640978692568 + 7.326316917663921i)=-6.194760250156455e-09 - 3.142629031397849e-09i

Thanks a lot.

[Walter Roberson]

In your previous question, the behavior of fsolve() with respect to complex-valued functions was discussed, but that behavior is relatively new. You did not indicate which MATLAB version you are using; if you are using an older version then it is going to give inaccurate results.

[ttopal]

hi Walter, sorry I must have missed the question about the version. I am using R2016a. This question is sort of follow up. Since I couldn't get to find the roots in a range, i decided to check them manually and realized that I have been seeing false roots because of tolerances. I increased the accuracy and some disappeared. Yet, it is still far from perfect. The issue is, i have a good estimate on the real part of the roots and not so much about the imaginary part. That is why, i give a real valued initial guess and look for complex root. If I can omit the false roots with the help of norm of step, that should solve my case. ( at the moment i check them manually)

Answer 1

[ttopal]

So I have found that my own root finder with improvement based on Newton-Raphson method with numerical approximations to the derivative works way better. Here is the code;

%
% Newton-Raphson method with numerical approximations to the derivative.
%
function [x0,flag] = newton_root(myfun,x0)
maxit = 15;
tol = 1.0e-15;
err = 100;
icount = 0;
xold =x0;
%--- search rooth from x0-range to x0+range---
check=x0;
range=0.2;
while (err > tol && icount <= maxit)
 icount = icount + 1;
 flag=0;
 f = myfun(xold);
 h = min(0.0001*xold,0.0001);
 fp = myfun(xold+h);
 fn = myfun(xold-h);
 df = (fp - fn)/(2*h);
 xnew = xold - f/df;
 if ((check-range-xnew)*(xnew-check-range) < 0.0)
%      fprintf(1,'Sorry. You jumped out of interval');
%      fprintf(1,'The final value of x was %e \n', x0);
     flag=1;
     break
 end
 err = abs((xnew-xold)/xnew);

   %fprintf(1,'icount = %i xold = %e f = %e df = %e xnew = %e err = %e \n',icount, xold, f, df, xnew, err);
   %fprintf(1,'%i %e %e %e %e %e \n',icount, xold, f, df, xnew, err);
   xold = xnew;
  end
  %
  x0 = xnew;
  if (icount >= maxit)
   % you ran out of iterations

   % fprintf(1,'Sorry. You did not converge in %i iterations.\n',maxit);
   % fprintf(1,'The final value of x was %e \n', x0);
   flag=2;
  end

Answer 2

[David Ding]

Hello Turker,

In your second case, the result is not so much of a "false root" as instead a less accurate answer than the first case (1e-9 vs 1e-15). What happened was the solver determined that the "first-order optimality" was small enough and thus terminated the iterations. Since first-order optimality is a necessary but not sufficient condition for solving the equation, you might not obtain a result that you expected.

A possible workaround to avoid seeing this discrepancy is to raise both the "TolFun" and "TolX" parameters to above 1e-14. Setting small tolerances does not guarantee accurate results and can affect solver convergence. Note that the default "TolFun" and "TolX" values are 1e-4.

More information about optimization settings can be found here:

https://www.mathworks.com/help/releases/R2016b/matlab/math/setting-options.html

Thanks,

David

[ttopal]

Hi David,

Thank you for your answer, but I had already set both TolFun and TolX values to 1e-16.

[Alan Weiss]

As documented, when you set a tolerance to below eps ~ 2e-16, then you disable the tolerance; it is as if you set it to zero. This can cause convergence issues, as the solver can fail to recognize when it should stop.

Basically, you seem to want an answer that is exactly zero, but floating point math may prevent you from ever achieving this result. I would say that fsolve found a correct root at 3.782640978692568 + 7.326316917663921i.

Alan Weiss

MATLAB mathematical toolbox documentation

[ttopal]

Thank you Alan. I see what you mean with the tolerance set. And, I agree that 3.782640978692568 + 7.326316917663921i seems valid root. Nevertheless, I am still seeking for an answer to my initial question. Do i have an option to limit results based on norm of step?