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?
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.
