Bairstow Method to find polynomial roots matlab code problem

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;

Answer 1

[Wayne King]

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()

[Wayne King]

I mean should contain the same roots, I don't know how they are ordered since I'm just going by the help.

Answer 2

[Walter Roberson]

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

[Steve]

Walter please help me with this I really need this.

Answer 3

[Shadi Srm]

What's the solution if in oneof the steps dn=0?

Answer 4

[Shumet]

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

[Steven Lord]

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.