I need to change this code into one that can solve complex roots

6 views (last 30 days)
I have this Newton-Raphson method code, but I need to change it into a code that can actually solve polynomials and give as an answer complex roots.
Code:
function [tabla, raiz]=newtonraphsonMN(f,xa,errorD,imax)
tabla=[];
error=Inf;
i=0;
xr=NaN;
df=diff(f);
while error>errorD && i<=imax
i=i+1;
fxa=double(subs(f,xa));
fpxa=double(subs(df,xa));
xr=xa-fxa/fpxa;
error=100*abs((xr-xa)/xr);
tabla=[tabla; [xa fxa fpxa xr error]];
xa=xr;
end
%raiz
raiz=xr;

Accepted Answer

John D'Errico
John D'Errico on 12 Sep 2022
Edited: John D'Errico on 13 Sep 2022
Interesting. It does not work? :) Gosh, You could have fooled me. I'll try an example. (I've attached the code you gave, so it will be used.)
syms X
f = X^2 + X + 1;
Does f has complex roots?
solve(f)
ans = 
vpa(ans)
ans = 
Of course. I'd not have used an example that lacks complex roots, since that is your question.
[tabla, raiz]=newtonraphsonMN(f,1 + i,1e-12,100);
Strange.
tabla
tabla =
1.0e+02 * 0.0100 + 0.0100i 0.0200 + 0.0300i 0.0300 + 0.0200i 0.0008 + 0.0062i 1.6125 + 0.0000i 0.0008 + 0.0062i 0.0070 + 0.0071i 0.0115 + 0.0123i -0.0052 + 0.0063i 0.7267 + 0.0000i -0.0052 + 0.0063i 0.0035 - 0.0002i -0.0003 + 0.0126i -0.0049 + 0.0091i 0.2687 + 0.0000i -0.0049 + 0.0091i -0.0008 + 0.0001i 0.0001 + 0.0182i -0.0050 + 0.0087i 0.0424 + 0.0000i -0.0050 + 0.0087i -0.0000 + 0.0000i 0.0000 + 0.0173i -0.0050 + 0.0087i 0.0010 + 0.0000i -0.0050 + 0.0087i -0.0000 + 0.0000i 0.0000 + 0.0173i -0.0050 + 0.0087i 0.0000 + 0.0000i -0.0050 + 0.0087i -0.0000 + 0.0000i 0.0000 + 0.0173i -0.0050 + 0.0087i 0.0000 + 0.0000i -0.0050 + 0.0087i 0.0000 + 0.0000i 0.0000 + 0.0173i -0.0050 + 0.0087i 0.0000 + 0.0000i
format long
raiz
raiz =
-0.500000000000000 + 0.866025403784439i
To me, it seems to have worked. But then, what do I know? :)
(Hint: Do you see what I did different?)
  2 Comments
Torsten
Torsten on 13 Sep 2022
To solve for polynomial roots, use MATLAB's "roots" command.
Newton's method will only give you one solution per call. Which one you get depends on the initial guess.
Example:
To determine all roots of p(x) = x^2+x+1:
p = [1 1 1];
sol = roots(p)
sol =
-0.5000 + 0.8660i -0.5000 - 0.8660i
polyval(p,sol)
ans = 2×1
1.0e-15 * 0.3331 0.3331

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!