Metoda Bairstowa — wybór początkowych wartości, zapewnienie zbieżności

1 
uses crt;
type TElem = double;
     TArray = array of TElem;
const EPS = 1e-12;

{
   Procedura znajdujaca czynnik kwadratowy
   Znajduje takze wspolczynniki wielomianu po deflacji znalezionym czynnikiem kwadratowym
   Korzysta z metody stycznych Newtona aby rozwiazac uklad rownan nieliniowych
   powstaly po przyrownaniu do zera reszty z dzielenia wielomianu przez trojmian kwadratowy

}

procedure Bairstow(deg:integer;a:TArray;var b:TArray;r0,s0:TElem;var r:TElem;var s:TElem);
var i:integer;
    d,dr,ds:TElem;
    c:TArray;
begin
   SetLength(c,deg+1);
   repeat
     b[deg] := a[deg];
     b[deg-1] := a[deg-1]+r0*b[deg];
     c[deg] := 0;
     c[deg-1] := b[deg];
     for i := deg - 2 downto 0 do
     begin
       b[i] := a[i] + r0 * b[i + 1] + s0 * b[i + 2];
       c[i] := b[i + 1] + r0 * c[i + 1] + s0 * c[i + 2];
     end;
     d := (c[0] * c[2] - sqr(c[1]));
     dr := (c[1] * b[1] - c[2] * b[0])/d;
     ds := (c[1] * b[0] - c[0] * b[1])/d;
     r0 := r0 + dr;
     s0 := s0 + ds;
   until ((abs(dr)<EPS)and(abs(ds) < EPS));
   r := r0;
   s := s0;
end;

{
    Procedura rozwiazujaca numerycznie rownanie wielomianowe
    wykorzystujaca czynnik kwadratowy i wspolczynniki wielomianu
    po deflacji znalezionym czynnikiem kwadratowym
    Powyzsze dane otrzymujemy w wyniku dzialania procedury Bairstow
    Procedura znajduje wszystkie pierwiastki wielomianu takze te zespolone
    unikajac arytmetyki zespolonej przy zalozeniu ze wspolczynniki wielomianu sa rzeczywiste
}

procedure PolyRoots(deg:integer;a:TArray;p0,q0:TElem;var xr:TArray;var xi:TArray);
var i:integer;
        b,c:TArray;
        p,q,d:TElem;
begin
  SetLength(b,deg+1);
  SetLength(c,deg+1);
  for i := 0 to deg do
      c[i] := a[i]/a[deg];
  for i := 0 to deg-1 do
  begin
    xr[i] := 0;
    xi[i] := 0;
  end;
  while (deg >= 2) do
  begin
    Bairstow(deg,c,b,p0,q0,p,q);
    d := sqr(p) + 4*q;
    if (d < 0) then
    begin
      xr[deg-1] := 0.5 * p;
      xr[deg-2] := 0.5 * p;
      xi[deg-1] := -0.5 * sqrt(-d);
      xi[deg-2] := 0.5 * sqrt(-d);
    end
    else
    begin
      xr[deg-1] := 0.5 * (p - sqrt(d));
      xr[deg-2] := 0.5 * (p + sqrt(d));
    end;
    deg := deg - 2;
    for i := 0 to deg do
        c[i] := b[i+2];
  end;
  if (deg = 1) then
     xr[0] := -c[0]/c[1];
end;

{
  Glowny blok programu demonstrujacy dzialanie procedur
}

var
    i,deg:integer;
    a,xr,xi:TArray;
    p0,q0:TElem;
    esc : char;

begin
  ClrScr;
  WriteLn('Przyblizone rozwiazywanie rownan wielomianowych metoda Bairstowa');
  repeat
  WriteLn('Podaj stopien wielomianu');
  ReadLn(deg);
  SetLength(a,deg+1);
  SetLength(xr,deg);
  SetLength(xi,deg);
  WriteLn('Podaj poczatkowe wartosci wspolczynnikow trojmianu');
  ReadLn(p0,q0);
  WriteLn('Podaj wspolczynniki wielomianu');
  for i := deg downto 0 do
  begin
    Write('P[',i,']=');
    ReadLn(a[i]);
  end;
  PolyRoots(deg,a,p0,q0,xr,xi);
  WriteLn('Przyblizone wartosci pierwiastkow wielomianu to: ');
  for i := 0 to deg - 1 do
    if xi[i] < 0 then
        WriteLn('x[',i,']=',xr[i]:1:12,xi[i]:1:12,'*i')
    else
       WriteLn('x[',i,']=',xr[i]:1:12,'+',xi[i]:1:12,'*i');
    esc := ReadKey;
  until esc = #27;
end.

Co do mozliwosci modyfikacji to:

  • Czy istnieje jakis sposob wyboru wartosci poczatkowych?
  • Czy mozna zmodyfikowac badz zastapic metode Newtona tak aby zagwarantowac zbieznosc procedurze Bairstow?
0

Ostatnio zmodyfikowałem powyższy kod programu aby metody Bairstowa używał tylko dla równania wielomianowego parzystego stopnia
ale nadal problemem jest metoda Newtona i powyższe dwa pytania pozostają otwarte

uses crt;
type TElem = double;
     TArray = array of TElem;
const EPS = 1e-12;

procedure Horner(deg:integer;a:TArray;z:TElem;var b:TArray;var alpha:TElem;var beta:TElem);
var k:integer;
begin
  alpha := a[deg];
  beta := 0;
  for k := deg-1 downto 0 do
  begin
    b[k] := alpha;
    beta := alpha + z * beta;
    alpha := a[k] + z * alpha;
  end;
end;

procedure Bairstow(deg:integer;a:TArray;var b:TArray;u0,v0:TElem;var u:TElem;var v:TElem);
var
    k:integer;
    d,du,dv:TElem;
    c:TArray;
begin
  SetLength(c,deg+1);
  repeat
     b[deg] := a[deg];
     c[deg] := 0;
     c[deg-1] := a[deg];
     b[deg-1] := a[deg-1]+u0*b[deg];
     for k := deg - 2 downto 0 do
     begin
       b[k] := a[k] + u0 * b[k + 1] + v0 * b[k + 2];
       c[k] := b[k + 1] + u0 * c[k + 1] + v0 * c[k + 2];
     end;
     d := c[0] * c[2] - sqr(c[1]);
     du := (c[1] * b[1] - c[2] * b[0])/d;
     dv := (c[1] * b[0] - c[0] * b[1])/d;
     u0 := u0 + du;
     v0 := v0 + dv;
  until ((abs(du) < EPS)and(abs(dv) < EPS));
  u := u0;
  v := v0;
end;


procedure PolyRoots(deg:integer;a:TArray;p0,q0:TElem;var xr:TArray;var xi:TArray);
var i:integer;
    b,c:TArray;
    p,q,d:TElem;
    stop:boolean;
begin
  SetLength(b,deg+1);
  SetLength(c,deg+1);
  for i := 0 to deg do
     c[i] := a[i]/a[deg];
  for i := 0 to deg - 1 do
  begin
    xr[i] := 0;
    xi[i] := 0;
  end;
  if odd(deg) then
  begin
    d := 0;
    for i := 0 to deg - 1 do
       if(abs(c[i]) > d) then
          d := abs(c[i]);
    d := d + 1;
    repeat
       Horner(deg,c,d,b,p,q);
       xr[deg-1] := d - p/q;
       stop := (abs(xr[deg-1] - d) < EPS);
       d := xr[deg-1];
    until stop;
    deg := deg - 1;
    for i := 0 to deg do
       c[i] := b[i];
  end;
  while (deg >= 2) do
  begin
    Bairstow(deg,c,b,p0,q0,p,q);
    d := sqr(p)+4*q;
    if (d < 0) then
    begin
      xr[deg - 1] := 0.5 * p;
      xr[deg - 2] := 0.5 * p;
      xi[deg - 1] := -0.5 * sqrt(-d);
      xi[deg - 2] := 0.5 * sqrt(-d);
    end
    else
    begin
      xr[deg - 1] := 0.5 * (p - sqrt(d));
      xr[deg - 2] := 0.5 * (p + sqrt(d));
    end;
    deg := deg - 2;
    for i := 0 to deg do
       c[i] := b[i + 2];
  end;
end;


var
    i,deg:integer;
    a,xr,xi:TArray;
    p0,q0:TElem;
    esc : char;

begin
  ClrScr;
  WriteLn('Przyblizone rozwiazywanie rownan wielomianowych metoda Bairstowa');
  repeat
    WriteLn('Podaj stopien wielomianu');
    ReadLn(deg);
    SetLength(a,deg+1);
    SetLength(xr,deg);
    SetLength(xi,deg);
    WriteLn('Podaj poczatkowe wartosci wspolczynnikow trojmianu kwadratowego');
    ReadLn(p0,q0);
    WriteLn('Podaj wspolczynniki wielomianu');
    for i :=deg downto 0 do
    begin
      Write('P[',i,']=');
      ReadLn(a[i]);
    end;
    PolyRoots(deg,a,p0,q0,xr,xi);
    WriteLn('Przyblizone wartosci pierwiastkow wielomianu to: ');
    for i := 0 to deg - 1 do
        if xi[i] < 0 then
            WriteLn('x[',i,']=',xr[i]:1:12,xi[i]:1:12,'*i')
        else
            WriteLn('x[',i,']=',xr[i]:1:12,'+',xi[i]:1:12,'*i');
    esc := ReadKey;
    until esc = #27;
end.

1 użytkowników online, w tym zalogowanych: 0, gości: 1