MATLAB Answers

Multiple roots formula with Matlab

7 views (last 30 days)
uruc
uruc on 24 Oct 2020
Commented: Alan Stevens on 26 Oct 2020
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 Comments

Sign in to comment.

Accepted Answer

Alan Stevens
Alan Stevens on 24 Oct 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
Display your result
disp(x)
Be careful not to use an initial guess that makes fp equal zero, or the routine will fail!

  3 Comments

uruc
uruc on 25 Oct 2020
Alan Stevens,
Thank you for your answer, but we had a hard time writing the iteration formula in matlab.
Even though we found some answers, we couldn't be sure if we were on the right track or not.
Could you please show us how to put the formula and see what are the first 3 or 4 roots of the problem?
Thanks in advance.
Alan Stevens
Alan Stevens on 25 Oct 2020
It's often a good idea to plot the function first to get an idea of where the roots are. Here the plot shows:
You have a fourth order polynomial, so you should have 4 roots (places where the function hits zero). This function looks like it has a roots at or near 1 and 3. At least one of these is a multiple root. That is, it is repeated more than once. That is why your iteration routine is a multi-root routine. You need to supply an initial guess to get the actual vaues of the roots. Choose different initial guesses to get the two roots. Even though there are technically four roots, because of the repetition, the routine will only spit out one (two in total, depending on the initial guess).
The iteration expression should look like
x(i+1) = x(i) - f(x(i))*fp(x(i))/( fp(x(i))^2 - f(x(i))*fpp(x(i)) ) ;
Now let's see some coding from you.
John D'Errico
John D'Errico on 25 Oct 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.

Sign in to comment.

More Answers (1)

uruc
uruc on 26 Oct 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

Products


Release

R2020b

Community Treasure Hunt

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

Start Hunting!

Translated by