Function roots. Fixed-point iteration
4 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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
0 commentaires
Réponses (2)
John D'Errico
le 20 Mai 2019
Modifié(e) : John D'Errico
le 20 Mai 2019
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.
1 commentaire
ASHA RANI
le 30 Mai 2021
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
pragiedruliai
le 21 Mai 2019
1 commentaire
Steven Lord
le 21 Mai 2019
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".)
Voir également
Catégories
En savoir plus sur Assumptions dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!