Bairstow Method
- To: mathgroup at smc.vnet.net
- Subject: [mg71220] Bairstow Method
- From: "ms z" <ms-z- at hotmail.com>
- Date: Fri, 10 Nov 2006 06:38:46 -0500 (EST)
I'm sorry for mailing again. Left out something in the previous email.
I posted a question on how to write an automated program to solve the
polynomial f(x)=(x-4)(x+2)(x-1)(x+5)(x-7) (without using NSolve or Solve)
earlier on.
I suggested this program: (though not a good one)
solution := {Plot[{(x - 4)(x + 2)(x - 1)(x + 5)(x - 7)}, {x, -10, 10},
AxesLabel -> TraditionalForm /@ {x, y}]}
I've tried to write another program using Bairstow Method. But it doesn't
not seem to work. Could I have some help?
The program is as follows:
Off[General::spell];
bairstow[rl_, sl_, esl_] :=
Block[{r, x, rnew, snew, es, f,
f1, f2, cl1, cl2, dr, ds, eqn1, eqn2, ear, eas},
f = Expand[(x - 4) (x + 2) (x - 1) (x + 5) (x - 7)];
ea = 100; es = esl;
For[i = 1, ear && eas > es, i++,
(b5 = 1;
b4 = -5 + r*b5;
b3 = -35 + r*b4 + s*b5;
b2 = 125 + r*b3 + s*b4;
b1 = 194 + r*b2 + s*b3;
b0 = -280 + r*b1 + s*b2;
c5 = b5;
c4 = b4 + r*c5;
c3 = b3 + r*c4 + s*c5;
c2 = b2 + r*c3 + s*c4;
c1 = b1 + r*c2 + s*c3;
eqn1 = c2*dr + c3*ds == -b1;
eqn2 = c1*dr + c2*ds == -b0;
dr = (-b1*c2 - (-b0)*c3)/(c2*c2 - c1*c3);
ds = (-b1*c1 - (-b0)*c2)/(c3*c1 - c2*c2);
If[i == 1, (r = rl; s = sl), (r = rnew; s = snew)];
dr; ds;
rnew = r + dr;
snew = s + ds;
ear = Abs[dr/rnew]*100;
eas = Abs[ds/snew]*100;
Clear[r, s, dr, ds];)];
x1 = (rnew + Sqrt[(rnew^2 + 4*snew)])/2;
x2 = (rnew - Sqrt[(rnew^2 + 4*snew)])/2;
Print["Value of x1 is ", x1];
Print["Value of x2 is ", x2];
Clear[r, s, dr, ds];
f = x^5 - 5 x^4 - 35 x^3 + 125 x^2 + 194 x - 280;
f1 = PolynomialQuotient[f, x^2 - rnew*x - snew, x];
cl1 = CoefficientList[f1, x];
Clear[rnew, snew];
ea = 100; es = esl;
For[i = 1, ear && eas > es, i++,
(b3 = cl1[[4]];
b2 = cl1[[3]] + r*b3;
b1 = cl1[[2]] + r*b2 + s*b3;
b0 = cl1[[1]] + r*b1 + s*b2;
c3 = b3;
c2 = b2 + r*c3;
c1 = b1 + r*c2 + s*c3;
eqn1 = c2*dr + c3*ds == -b1;
eqn2 = c1*dr + c2*ds == -b0;
dr = (-b1*c2 - (-b0)*c3)/(c2*c2 - c1*c3);
ds = (-b1*c1 - (-b0)*c2)/(c3*c1 - c2*c2);
If[i == 1, (r = rl; s = sl), (r = rnew; s = snew)];
dr; ds;
rnew = r + dr;
snew = s + ds;
ear = Abs[dr/rnew]*100;
eas = Abs[ds/snew]*100;
Clear[r, s, dr, ds];)];
x3 = (rnew + Sqrt[(rnew^2 + 4*snew)])/2;
x4 = (rnew - Sqrt[(rnew^2 + 4*snew)])/2;
Print["Value of x3 is ", x3];
Print["Value of x4 is ", x4];
Clear[r, s, dr, ds];
f2 = PolynomialQuotient[f1, x^2 - rnew*x - snew, x];
cl2 = CoefficientList[f2, x];
x5 = cl2[[1]]/cl2[[2]];
Print["Value of x5 is ", x5];]
_________________________________________________________________
Find just what you are after with the more precise, more powerful new MSN
Search. http://search.msn.com.sg/ Try it now.