fixed point Iterative method for finding root of an equation

595 views (last 30 days)
Milind Amga on 30 Sep 2020
Answered: Akash G M on 19 Dec 2020
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?
% 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
John D'Errico on 30 Sep 2020
@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.

John D'Errico on 30 Sep 2020
Edited: John D'Errico on 30 Sep 2020
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 on 30 Sep 2020
Thanks for the explanation. It is very helpful. I understood what you are trying to convey.

Alan Stevens on 30 Sep 2020
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 on 30 Sep 2020
Thanks. This code works for logarithmic and exponential function. Really appreciate your effort.

Akash G M on 19 Dec 2020
Find the root of the equation cos(x) - 1.3x = 0, taking initial approximation as 0.2 by fixed point iteration method