Bairstow Method to find polynomial roots matlab code problem
11 views (last 30 days)
Show older comments
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 Comments
Accepted Answer
Wayne King
on 10 Oct 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 Comment
Wayne King
on 10 Oct 2011
I mean should contain the same roots, I don't know how they are ordered since I'm just going by the help.
More Answers (3)
Walter Roberson
on 10 Oct 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
on 27 Oct 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 Comment
Steven Lord
on 27 Oct 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.
See Also
Categories
Find more on Polynomials in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!