fixed point Iterative method for finding root of an equation

The code below gives the root and the iteration at which it occur. The code goes into an infinite loop when the function contains any logarithmic or exponential function. But works fine with algebraic and trignometric function.

Can anyone point out the error?

Thank you in advance.

% let the equation whose root we want to find be x^3-5*x-7 = 0
% simplified eqation example:-  f = @(x) (5x+7)^(1/3)
function [root,iteration] = fixedpoint(a,f) %input intial approiximation and simplified form of function  
if nargin<1 % check no of input arguments and if input arguments is less than one then puts an error message
    fprintf('Error! Atleast one input argument is required.');
    return;
end
if nargin<2 % check no of input arguments and if input arguments is less than two then puts a=0
    a=0;
end
x(1) = f(a) ;
i =1 ;
temp = false;
while temp == false 
    x(i+1) = f(x(i));
    if x(i+1) == x(i)
        temp = true;
        root = x(i-1); 
        iteration = i-1;
        return;
    end
    i = i +1 ;
end
x
end

[Sindar]

I don't understand your code. It looks like you are finding locations where f(x)=x, not roots

what about this case:

f(x) = 2x+1

a=0

x(1) = f(a) = 1

x(2) = f(x(1)) = f(1) = 3

...

x(n) = 2n-1

infinite loop

[Alan Stevens]

@ Sindar In the following

% let the equation whose root we want to find be x^3-5*x-7 = 0
% simplified eqation example:-  f = @(x) (5x+7)^(1/3)

the function is a rearranged form of x^3 - 5x - 7 = 0 so x = (5x + 7)^(1/3).

@ Milind This one does work. However, it's not always so easy to find a rearrangement that does work. Fixed point iterations often diverge. You might have to try different rearrangements of your equation to get the procedure to work (even then, there is no guarantee).

[Milind Amga]

I tried different arrangements for below equation :-

x-exp(-x)=0

But it runs into an infinite loop. To the best of my knowledge, the value of derivative of the simplified form of above equation at initial approximation should be less than 1 for better convergence. So i simplified the above equation to meet this condition but for this case code runs into an infinite loop.

Any idea what could cause the problem?

I appreciate your effort.

[John D'Errico]

Yes. There is a guarantee. Let me do some writing.

[Alan Stevens]

@ John

Nice explanation. In practice, with a complicated function, it can be extremely difficult to ensure an arrangement such that abs(derivative near root)<1. Is it theoretically possible always to find such an arrangement?

[John D'Errico]

Admittedly, if it was easy to do in some completely general way, we would probably see a fixed point solver in the optimization toolbox, or at least in MATLAB proper as a complement to tools like fminbnd and fzero.

You need to find a transformation in a fixed point form of the function around the root, where the absolute derivative is bounded by 1. And to be able to do that, you want to know where the root of interest is, and also be able to differentiate the function, and to know the derivative is bounded. But wait! If you knew all of that, then just use a Newton method. Or, better yet, a tool like fzero, which will be more robust to some problems and still has good convergence properties.

In the end, the answer really is to just use fzero, or whatever solver is appropriate. That is what I try to preach time and again - that while learning to use methods like fixed point iteration is a good thing for a student, after you get past being a student, use the right tools and don't write your own.

But can we use fixed point on some general problem? Lets see. find a root of the quadratic function x^2-3*x+2. Yeah, I know, creative as hell on my part. :)

syms x
solve(x^2 - 3*x + 2)
ans =
 1
 2

Two roots where I would expect to find them. Will my fixed point solver converge? Which root will it work for?

F = @(x) x.^2 - 3*x + 2;

We can make a good guess from this plot:

syms x
fplot(diff(x^2 - 3*x + 2) + 1)
yline(-1,'r');
yline(1,'r');
xline(1,'g')
xline(2,'g')

I've plotted the derivative of my fixed point function, along with horizontal bands at +/- 1.

As you can see, the root at x==1 is within the band, but the root at x==2 is not.

So, I'll predict my solver will converge for the root at x==1, but NOT at x==2.

>> [xfinal,fval,ferr,itercount] = myfp(F,1.1,0.001,10,1);
Current x, Fval, Ferr
                      1.01       -0.0899999999999999        0.0899999999999999
                    1.0001      -0.00990000000000002       0.00990000000000002
                1.00000001     -9.99900000002718e-05      9.99900000002718e-05

>> [xfinal,fval,ferr,itercount] = myfp(F,2.1,0.001,10,1);
Current x, Fval, Ferr
                      2.21         0.109999999999999         0.109999999999999
                    2.4641         0.254099999999998         0.254099999999998
          3.14358880999999         0.679488809999997         0.679488809999997
           5.5949729863572           2.4513841763572           2.4513841763572
          22.1137767453524          16.5188037589952          16.5188037589952
          446.791568452582           424.67779170723           424.67779170723
          198731.122503413          198284.330934961          198284.330934961
          39493661591.2216          39493462860.0991          39493462860.0991
      1.55974930580295e+21      1.55974930576345e+21      1.55974930576345e+21
      2.43281789695278e+42      2.43281789695278e+42      2.43281789695278e+42

[Milind Amga]

@John

Well, I am new to matlab. Started learning matlab a few months back. My Math teacher asked me to write a code for fixed point iteration method. So i reached out for help here.

Anyways, thanks again for explanation.

[John D'Errico]

@Milind Amga - you misunderstood my point, which was to Alan's question.

You would rarely want to use fixed point iteration in the real world. (I have, but only on rare occasions.) Yes, you are solving this as requested because it was homework. It teaches you both about MATLAB as well as teaching you one method for solving a problem, even though that method is not really a very good one in practice, because it has some serious flaws as outlined. (In the end, I hope your teacher spends the time for you to understand why the method works, as well as the limitations on the method and why you probably want to use better methods in the end. It is the latter parts that often get skipped over.)

I made the statement about fzero because I see far too many scientists and engineers using crude methods they learned in school to solve a problem, not knowing that good code already exists to do the same thing. For example, long ago, I was asked to help a scientist to do a computation using determinants to solve a problem. They were solving 2x2 systems of equations, and they wanted to change it so they could solve 5x5 or 6x6 systems. But they did not know how to form the necessary higher order determinants.

My answer was no, I would not help them to do what they asked me to do, but that I would show them how to solve the problem using better methods.

Answer 1

[John D'Errico]

To understand fixed point iteration, we need to know why and when it will diverge. When does it converge? For some problems, the iterations might never be absolutely convergent, at least not without getting creative.

You comment about "better convergence". It is a question of divergence. Consider the trivial problem,

f = @(x) 1.1*x - 2.2;

Yes, we know the solution is x == 2. But does MATLAB?

function [xfinal,fval,ferr,itercount] = myfp(F,x0,ftol,maxit,verbosity)
% [xfinal,fval,itercount] = myfp(F,x0,ftol,maxit)
% 
%  myfp: solve for a root of the function F(x) == 0 using fixed point
%  iteration. This will be achieved using the iteration as
%    x = F(x) + x
%
% Arguments: (input)
% F - function of one variable
% x0 - initial value
% ftol - tolerance on the function value being zero
% maxit - maximum number of iterations allowed
% verbosity - (optional) flag to define output behavior, 0 = none, 1 = stuff on screen
% 
% xfinal - final estimate of the root
% fval - function value at that point
% itercount - the number of iterations taken
 
if (nargin < 5) || isempty(verbosity)
  % the default is to be comatose
  verbosity = 1;
end
 
% initialize for iteration
itercount = 0;
xfinal = x0;
fval = F(x0);
ferr = abs(fval);
 
if verbosity
  disp('Current x, Fval, Ferr')
end
 
% just a while loop that goes until it runs out of time, or it gets lucky
% and "converges"
while (ferr > ftol) && (itercount < maxit)
  itercount = itercount + 1;
  
  xprev = xfinal;
  
  % the iteration itself
  fval = F(xfinal);
  xfinal = fval + xfinal;
  
  % how far from 0 is the current iteration?
  ferr = abs(fval);
  
  if verbosity
    disp([xfinal,fval,ferr])
  end
  
end

That should work. As I set it up, I am solving the problem as

F(x) - x + x == 0

and therefore we have the iterations as

x = F(x) + x

Now, let me try it in MATLAB.

[xfinal,fval,ferr,itercount] = myfp(F,1.5,0.001,5,1);
Current x, Fval, Ferr
                      0.95                     -0.55                      0.55
                    -0.205                    -1.155                     1.155
                   -2.6305                   -2.4255                    2.4255
                  -7.72405                  -5.09355                   5.09355
                -18.420505                -10.696455                 10.696455

So it clearly diverges. The longer I allow it to run, the further it will go, even if I start quite close to the solution (unless of course it is within the tolerance at the start point.)

So let me fix it by just dividing by 2.

F = @(x) (1.1*x - 2.2)/2;
[xfinal,fval,ferr,itercount] = myfp(F,1.5,0.001,5,1);
Current x, Fval, Ferr
                     1.225                    -0.275                     0.275
                   0.79875                  -0.42625                   0.42625
                 0.1380625                -0.6606875                 0.6606875
              -0.886003125              -1.024065625               1.024065625
            -2.47330484375            -1.58730171875             1.58730171875

Hmm. Still diverges. Why?

Because the derivative of the right hand side as I created the fixed point iteration must have the property

abs(diff(F(x) + x)) < 1

in the vicinity of the root. If that derivative is larger than 1, divergence occurs. The scheme will wander away.

But mow let me change it to be:

F = @(x) (1.1*x - 2.2)/(-2);
>> [xfinal,fval,ferr,itercount] = myfp(F,1.5,0.001,5,1);
Current x, Fval, Ferr
                     1.775                     0.275                     0.275
                   1.89875                   0.12375                   0.12375
                 1.9544375        0.0556875000000001        0.0556875000000001
               1.979496875        0.0250593749999999        0.0250593749999999
             1.99077359375             0.01127671875             0.01127671875
             
[xfinal,fval,ferr,itercount] = myfp(F,1.5,0.001,100,0)
xfinal =
          1.99962165967871
fval =
      0.000462415948242256
ferr =
      0.000462415948242256
itercount =
      9

We can see it is clearly converging now, and in fact, required only 9 iterations.

Now, suppose we wish to find a root of the function f(x) == exp(x) - 2. My code does that by formulating the iteration as

x = exp(x) - 2 + x

And of course, since the derivative of the right hand side is just

exp(x) + 1

then fixed point iteratiion must always diverge. The starting value will not matter, unless it is EXACTLY at log(2). and even then, even the tiniest difference in the least significant bits will start to push it away from the root. The value of ftol would save you there though.

F = @(x) exp(x) - 2;

[xfinal,fval,ferr,itercount] = myfp(F,1.5,0.001,5,1);
Current x, Fval, Ferr
          3.98168907033806          2.48168907033806          2.48168907033806
          55.5891937168657          51.6075046465276          51.6075046465276
      1.38701157272672e+24      1.38701157272672e+24      1.38701157272672e+24
   Inf   Inf   Inf
   Inf   Inf   Inf

So things got bad very fast. Modify the function simply as

F2 = @(x) -0.5*F(x)

[xfinal,fval,ferr,itercount] = myfp(F2,1.5,0.001,5,1);
Current x, Fval, Ferr
         0.259155464830968         -1.24084453516903          1.24084453516903
         0.611237841842661         0.352082377011693         0.352082377011693
          0.68988235566251        0.0786445138198487        0.0786445138198487
         0.693141856814415       0.00325950115190499       0.00325950115190499
         0.693147180545774      5.32373135941899e-06      5.32373135941899e-06

As you can see, that final iteration will have converged.

But now, try to use the same trick to solve the related problem

F = @(x) exp(x) - 10;
>> F2 = @(x) -0.5*F(x);
>> [xfinal,fval,ferr,itercount] = myfp(F2,1.5,0.001,10,1);
Current x, Fval, Ferr
          4.25915546483097          2.75915546483097          2.75915546483097
         -26.1159481242028         -30.3751035890338          30.3751035890338
         -21.1159481242051          4.99999999999773          4.99999999999773
         -16.1159481245427          4.99999999966238          4.99999999966238
         -11.1159481746502          4.99999994989251          4.99999994989251
         -6.11595561126095          4.99999256338923          4.99999256338923
         -1.11705929394996            4.998896317311            4.998896317311
          3.71932035616684           4.8363796501168           4.8363796501168
          -11.898858916735         -15.6181792729018          15.6181792729018
         -6.89886231581379          4.99999660092118          4.99999660092118

And you should see it bounces around. It is easily fixed of course, as long as I choose a proper transformation of f.

F2 = @(x) -0.1*F(x);
[xfinal,fval,ferr,itercount] = myfp(F2,1.5,0.001,10,1);
Current x, Fval, Ferr
          2.05183109296619         0.551831092966194         0.551831092966194
          2.27361730438218         0.221786211415983         0.221786211415983
          2.30216954873883        0.0285522443566585        0.0285522443566585
          2.30258500666749      0.000415457928655094      0.000415457928655094
          
log(10)
ans =
          2.30258509299405

As you see, after 4 iterations it was quite happy.

The problem with an exponential like that however, is exp(x) gets arbitrarily large, so any constant multiplier will fail for some variation of the function. We could get tricky though, if we can find some transformation of F that will always have a small derivative of the fixed point iterant. since the derivative of exp(x) will always be large for some values of x, we need to be tricky.

[Milind Amga]

Thanks for the explanation. It is very helpful. I understood what you are trying to convey.

Answer 2

[Alan Stevens]

I defer to John on the guarantee. Certainly it's possible with x - exp(-x) = 0:

% x - exp(-x) = 0;  rearrange as x = exp(-x)
x0 = 0.5; % Initial guess
f = @(x) exp(-x);

tol = 10^-10;
flag = true;
its = 0;

while flag == true && its<100
    
    x = f(x0);
    if abs(x-x0)<tol
       flag = false;
    end
    x0 = x;
    its = its+1;
end

disp([x, x-exp(-x), its])

[Milind Amga]

Thanks. This code works for logarithmic and exponential function. Really appreciate your effort.

Answer 3

[Akash G M]

Find the root of the equation cos(x) - 1.3x = 0, taking initial approximation as 0.2 by fixed point iteration method

[Torsten]

% cos(x) - 1.3 x  = 0;  rearrange as x = cos(x)/1.3
x0 = 0.2; % Initial guess
f = @(x) cos(x)/1.3;
tol = 10^-10;
flag = true;
its = 0;
while flag == true && its<100    
    x = f(x0);
    if abs(x-x0)<tol
        flag = false;
    end
    x0 = x;
    its = its+1;
end
disp([x, x-f(x), its])

    0.6242    0.0000   29.0000