Function roots. Fixed-point iteration

Hello,

I can't figure out how to fix my fixed-point iteration method function funtions:

I have a function:

q= @(x) 0.0008*x.^7-0.0332*x.^6+0.5501*x.^5-4.7539*x.^4+23.5423*x.^3-68.9035*x.^2+110.8455*x-65.6061;
n=10;
F(q,5,10)

My function functions is:

function root = F(q, x0, n)
%UNTITLED2 Summary of this function goes here
%   Detailed explanation goes here
maxIter = 100;
epsilon = 5*10^-(n+1);
x=x0;
xold=x0;
for i=maxIter
    x=q(x);
    err = abs(x-xold);
    xold = x;
    if(err<epsilon)
        break
    end
end
root=x
end

Answer 1

[John D'Errico]

First, what are the roots? Are you trying to solve the problem q(x) == 0? Or are you trying to solve the problem q(x)==x?

A fixed point iteration as you have done it, implies that you want to solve the problem q(x) == x. So note that in the symbolic solve I use below, I subtracted off x from what you had as q(x).

syms x
fun = (0.0008*x.^7-0.0332*x.^6+0.5501*x.^5-4.7539*x.^4+23.5423*x.^3-68.9035*x.^2+110.8455*x-65.6061) - x;
format long g
double(solve(fun))
ans =
           1.25178388553229 +                     0i
           2.48825030999686 -      2.86450820415501i
           2.48825030999686 +      2.86450820415501i
           4.26318986807972 +                     0i
           8.67433389206779 -      1.70199422256665i
           8.67433389206779 +      1.70199422256665i
           13.6598578422587 +                     0i

So there are real roots near x=1.25, 4.26, and 13.66.

Now, what do we know about fixed point iteration? When will it be convergent?

https://en.wikipedia.org/wiki/Fixed-point_iteration

https://en.wikipedia.org/wiki/Fixed-point_theorem

http://wwwf.imperial.ac.uk/metric/metric_public/numerical_methods/iteration/fixed_point_iteration.html

A good rule for fixed point iteration is that near the root, the derivative should be less than 1 in absolute value. Does that hold near the roots?

q = 0.0008*x.^7-0.0332*x.^6+0.5501*x.^5-4.7539*x.^4+23.5423*x.^3-68.9035*x.^2+110.8455*x-65.6061;
double(subs(diff(q),x,[1.25,4.26,13.66]))
ans =
          17.9299775390625         -4.74152326450803          345.381167197486

I'd suggest it is highly unlikely that a fixed point iteration will converge. Now, how can you change the problem so convergence is likely? Rewrite it, by moving some terms around. Here, I just picked a term with a moderately large coefficient and a high power.

23.5423*x.^3 = x - (0.0008*x.^7-0.0332*x.^6+0.5501*x.^5-4.7539*x.^4+23.5423*x.^3-68.9035*x.^2+110.8455*x-65.6061)

Divide by the coefficient, then take the cube root. Now we have a fixed point iteration that looks like this:

x = nthroot((x - (0.0008*x.^7-0.0332*x.^6+0.5501*x.^5-4.7539*x.^4-68.9035*x.^2+110.8455*x-65.6061))/23.5423,3)

See if it works now.

x0 = [1 4 13];
FPI = @(x) nthroot((x - (0.0008*x.^7-0.0332*x.^6+0.5501*x.^5-4.7539*x.^4-68.9035*x.^2+110.8455*x-65.6061))/23.5423,3);
for i = 1:1000
  x0 = FPI(x0);
end
x0

x0 =
   1.25178388553228          1.25178388553229          13.6598578422554

So it looks like when we start near the root at 4.26, this variation still does not converge. But we manage to find the roots around 1.25 and 13.66.

The point is, fixed point iteration need not converge always. We would conclude that from a plot.

syms x
q = nthroot((x - (0.0008*x.^7-0.0332*x.^6+0.5501*x.^5-4.7539*x.^4-68.9035*x.^2+110.8455*x-65.6061))/23.5423,3);
fplot(diff(q),[0,15])
yline(1);
xline(1.25);
xline(4.26);
xline(13.66);

double(subs(diff(q),x,[1.25,4.26,13.66]))
ans =
         0.846215679769845          1.00448632683656         0.973868818386675

Fixed point iteration can be finicky. Sometimes you need to be creative about how you build an iteration so as to be convergent.

[ASHA RANI]

If fun involve a constant such as a in its expression who takes values like a=[1:1:10];

then what should be command instead of double(solve(fun))

syms x
fun = (a*0.0008*x.^7-0.0332*x.^6+0.5501*x.^5-4.7539*x.^4+23.5423*x.^3-68.9035*x.^2+110.8455*x-65.6061) - x;
format long g
double(solve(fun))
ans =
           1.25178388553229 +                     0i
           2.48825030999686 -      2.86450820415501i
           2.48825030999686 +      2.86450820415501i
           4.26318986807972 +                     0i
           8.67433389206779 -      1.70199422256665i
           8.67433389206779 +      1.70199422256665i
           13.6598578422587 +                     0i

Answer 2

[pragiedruliai]

Hello,

Thank you for your answer.

I need to find q(x) == 0.

Unfortunately my university has MATLAB Version: 8.6.0.267246 (R2015b), so I can't use syms x.

[Steven Lord]

Symbolic Math Toolbox is available for release R2015b, and the syms function has been part of that toolbox for a very long time (it was "Introduced before R2006a".)