# Multiple roots formula with Matlab

21 Ansichten (letzte 30 Tage)
uruc am 24 Okt. 2020
Beantwortet: Duncan Carlsmith am 6 Mär. 2023
f(x)= x^4- 6x^3 + 12x^2 - 10x +3
Starting point x0=0
Find x1,x2,x3
I need the find these points with the multiple roots formula = fi+1 = fi - (f(xi)*f '(xi)) / ([f '(xi)]^2 - f(xi)* f ''(xi))
Can anyone understand this and help us out, thanks a lot.
##### 0 Kommentare-1 ältere Kommentare anzeigen-1 ältere Kommentare ausblenden

Melden Sie sich an, um zu kommentieren.

### Akzeptierte Antwort

Alan Stevens am 24 Okt. 2020
I think your formula shoud be
First define your function and its first two derivatives in Matlab
f = @(x) x.^4- 6*x.^3 + 12*x.^2 - 10*x +3; % function
fp = @(x) ...; % first derivative here
fpp = @(x) ...; % second derivative here
Set an initial guess for x and a tolerance. Initialise an error measure and an iteration number. e.g:
x = 5; %initial guess
tol = 10^-6;
err = 1;
its = 0;
Then use a while loop
while err>tol && its<100
xold = x;
x = ... put your iteration formua here
err = abs(x-xold);
its=its+1;
end
disp(x)
Be careful not to use an initial guess that makes fp equal zero, or the routine will fail!
##### 3 Kommentare2 ältere Kommentare anzeigen2 ältere Kommentare ausblenden
John D'Errico am 25 Okt. 2020
Bearbeitet: John D'Errico am 25 Okt. 2020
A big point of this homework assignment may be to understand what happens near that multiple root.
Thus even if we look at the output of roots, we see three approximate roots at 1, but none seem to nail the root at 1.
>> roots([1 -6 12 -10 3])
ans =
2.99999999999999 + 0i
1.00000947488643 + 0i
0.999995262556785 + 8.20550203732316e-06i
0.999995262556785 - 8.20550203732316e-06i
You should see that typically, for a triple root, do not expect accuracty in the root as found to be better than roughly the cube root of eps. For a root r of multiplicity n, the resulting accuracy will typically be nthroot(eps( r ),n).
We can see why there is a problem. Here is your function:
>> fun = @(x) x.^4 - 6*x.^3 + 12*x.^2 - 10*x + 3;
>> fun(1 + [-1e-6 0 eps 1e-6])
ans =
-8.88178419700125e-16 0 8.88178419700125e-16 0
As you can see, if we try to evaluate fun at any point in the rough interval [1-1e-6 , 1+1e-6], we always get something as close to zero as MATLAB can resolve. We need to go quite a bit further out before it starts to show a different result than zero.
>> fun(1 + [-1e-4 1e-4])
ans =
2.00106597958438e-12 -1.99928962274498e-12
Next, you need to understand that a numerical rootfinder, once it finds that root around x == 1, need not understand the root is a multiple root. As far is it is concerned, a root is just a root.

Melden Sie sich an, um zu kommentieren.

### Weitere Antworten (2)

uruc am 26 Okt. 2020
f = @(x) x.^4- 6*x.^3 + 12*x.^2 - 10*x +3;
fp = @(x) 4*x^3 - 18*x^2 + 24*x - 10;
fpp = @(x) 12*x^2 - 36*x + 24;
x = 0;
tol = 10^-6;
err = 1;
its = 0;
while err>tol && its<10
xold = x;
x = x-f(x)*fp(x)/( fp(x)^2 - f(x)*fpp(x))
err = abs(x-xold);
its=its+1;
end
f = @(x) x.^4- 6*x.^3 + 12*x.^2 - 10*x +3;
fp = @(x) 4*x^3 - 18*x^2 + 24*x - 10;
fpp = @(x) 12*x^2 - 36*x + 24;
x = 0;
tol = 10^-6;
err = 1;
its = 0;
while err>tol && its<25
xold = x;
x = x-f(x)/fp(x)
err = abs(x-xold);
its=its+1;
end
We managed to calculte the results with these codes thanks to Alan Stevens.
Also thank you John D'Errico for the effort aswell. It was from another perspective which we appretiated.
Have a nice day
##### 1 KommentarKeine anzeigenKeine ausblenden
Alan Stevens am 26 Okt. 2020
Well done!

Melden Sie sich an, um zu kommentieren.

Duncan Carlsmith am 6 Mär. 2023
Not likely the point of the exercise but this works:
syms x; solve(x^4- 6*x^3 + 12*x^2 - 10*x +3==0,'ReturnConditions',true)
ans = struct with fields:
x: [4×1 sym] parameters: [1×0 sym] conditions: [4×1 sym]
ans.x
ans =
##### 0 Kommentare-1 ältere Kommentare anzeigen-1 ältere Kommentare ausblenden

Melden Sie sich an, um zu kommentieren.

### Kategorien

Find more on Matrix Indexing in Help Center and File Exchange

R2020b

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by