Bairstow Method to find polynomial roots matlab code problem
11 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Steve
el 10 de Oct. de 2011
Comentada: Steven Lord
el 27 de Oct. de 2023
Hello Experts,
I need matlab code of the Bairstow method to find polynomial roots.
I have found here on our site a guy who wrote such function - John Penny ( http://www.mathworks.com/matlabcentral/fileexchange/2305-numerical-methods-using-matlab-2e/content/edition2/na_funcs/bairstow.m). But when you use it with the simple polynomial coeff vector A = [1,-6,11,-6] (roots: 1,2,3) you don't get the right roots.
Please help me, I need it urgently to customize at work and that's why I need the correct code.
Here is the full function:
function [rts,it]=bairstow(a,n,tol)
% Bairstow's method for finding the roots of a polynomial of %degree n.
%
% Example call: [rts,it]=bairstow(a,n,tol)
% a is a row vector of REAL coefficients so that the
% polynomial is x^n+a(1)*x^(n-1)+a(2)*x^(n-2)+...+a(n).
% The accuracy to which the polynomial is satisfied is given by tol.
% The output is produced as an (n x 2) matrix rts.
% Cols 1 & 2 of rts contain the real & imag part of root respectively.
% The number of iterations taken is given by it.
%
it=1;
while n>2
%Initialise for this loop
u=1; v=1; st=1;
while st>tol
b(1)=a(1)-u; b(2)=a(2)-b(1)*u-v;
for k=3:n
b(k)=a(k)-b(k-1)*u-b(k-2)*v;
end;
c(1)=b(1)-u; c(2)=b(2)-c(1)*u-v;
for k=3:n-1
c(k)=b(k)-c(k-1)*u-c(k-2)*v;
end;
%calculate change in u and v
c1=c(n-1); b1=b(n); cb=c(n-1)*b(n-1);
c2=c(n-2)*c(n-2); bc=b(n-1)*c(n-2);
if n>3, c1=c1*c(n-3); b1=b1*c(n-3); end;
dn=c1-c2;
du=(b1-bc)/dn; dv=(cb-c(n-2)*b(n))/dn;
u=u+du; v=v+dv;
st=norm([du dv]); it=it+1;
end;
[r1,r2,im1,im2]=solveq(u,v,n,a);
rts(n,1:2)=[r1 im1]; rts(n-1,1:2)=[r2 im2];
n=n-2;
a(1:n)=b(1:n);
end;
%Solve last quadratic or linear equation
u=a(1); v=a(2);
[r1,r2,im1,im2]=solveq(u,v,n,a);
rts(n,1:2)=[r1 im1];
if n==2
rts(n-1,1:2)=[r2 im2];
end;
function [r1,r2,im1,im2]=solveq(u,v,n,a);
% Solves x^2 + ux + v = 0 (n 1) or x + a(1) = 0 (n = 1).
%
% Example call: [r1,r2,im1,im2]=solveq(u,v,n,a)
% r1, r2 are real parts of the roots,
% im1, im2 are the imaginary parts of the roots.
% Called by function bairstow.
%
if n==1
r1=-a(1);im1=0; r2=0; im2=0;
else
d=u*u-4*v;
if d<0
d=-d;
im1=sqrt(d)/2; r1=-u/2; r2=r1; im2=-im1;
elseif d>0
r1=(-u+sqrt(d))/2; im1=0; r2=(-u-sqrt(d))/2; im2=0;
else
r1=-u/2; im1=0; r2=-u/2; im2=0;
end;
end;
0 comentarios
Respuesta aceptada
Wayne King
el 10 de Oct. de 2011
I think you are most likely using the function incorrectly. Note from the help that the polynomial modeled by the function has a 1 for the highest power which is NOT included in the input vector, a.
so
A = [1,-6,11,-6];
roots(A)
is equivalent to
A1 = A(2:end);
[rts,it]=bairstow(A1,3,1e-3);
The first column of rts should be the same as the output of roots()
1 comentario
Wayne King
el 10 de Oct. de 2011
I mean should contain the same roots, I don't know how they are ordered since I'm just going by the help.
Más respuestas (3)
Walter Roberson
el 10 de Oct. de 2011
Notice that according to the documentation there is an implied leading 1 on the coefficients passed in. When you pass in A of [1,-6,11,-6] then you are generating the polynomial
x^4 + 1*x^3 - 6*x^2 + 11*x - 6
rather than your expected
1*x^3 - 6*x^2 + 11*x - 6
Shumet
el 27 de Oct. de 2023
f(x) = x^5 - 3.5x^4 + 2.75x^3 +2.125x^2 - 3.875x + 1.25 find the roots of the given polynomial using Bairstow's method with an approximate error ε restricted to 0.1%, initial guesses r = -1 and s = -1. using matlab
1 comentario
Steven Lord
el 27 de Oct. de 2023
This sounds like a homework assignment. If it is, show us the code you've written to try to solve the problem and ask a specific question about where you're having difficulty and we may be able to provide some guidance.
If you aren't sure where to start because you're not familiar with how to write MATLAB code, I suggest you start with the free MATLAB Onramp tutorial to quickly learn the essentials of MATLAB.
If you aren't sure where to start because you're not familiar with the mathematics you'll need to solve the problem, I recommend asking your professor and/or teaching assistant for help.
Rather than resurrecting a more than ten year old discussion, you might want to ask this again (with the code and the specific question) as a new question using the Ask link near the top of this page.
Ver también
Categorías
Más información sobre Polynomials en Help Center y File Exchange.
Productos
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!