unit2.pas

《Delphi常用数值算法集》的配书源码

unit Unit2;

interface
uses
  Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms, Dialogs;
Procedure QROOT(P:array of real; N:integer; var B, C:real; EPS:real);

implementation
procedure POLDIV(U:array of real; N:integer;V:array of real; NV:integer;
                         var Q, R:array of real);
var
    J,K:integer;
begin
    For J:=1 To N do
    begin
        R[J]:=U[J];
        Q[J]:=0 ;
    end;
    For K:=N - NV DownTo 0 do
    begin
        Q[K + 1]:=R[NV + K] / V[NV];
        For J:=NV + K - 1 DownTo K + 1 do 
            R[J]:=R[J] - Q[K + 1] * V[J - K];
    end;
    R[NV]:=0 ;
end;

Procedure QROOT(P:array of real; N:integer; var B, C:real; EPS:real);
const
    ITMAX = 20;  TINY = 0.000001;
var
    Q,REM1,QQ:array[0..20] of real; D:array[0..3] of real;
    ITER,I:integer;  S,SB,RB,R,SC,RC,DIV1,DELB,DELC,DB,DC:real;
begin 
    D[3]:=1; 
    For ITER:=1 To ITMAX do
    begin
        D[2]:=B;
        D[1]:=C;
        POLDIV(P, N, D, 3, Q, REM1);
        S:=REM1[1];
        R:=REM1[2];
        POLDIV(Q, N - 1, D, 3, QQ, REM1);
        SC:=-REM1[1];
        RC:=-REM1[2];
        For I:=N - 1 DownTo 1 do
            Q[I + 1]:=Q[I];
        Q[1]:=0; 
        POLDIV(Q, N, D, 3, QQ, REM1);
        SB:=-REM1[1];
        RB:=-REM1[2];
        DIV1:=1 / (SB * RC - SC * RB);
        DELB:=(R * SC - S * RC) * DIV1;
        DELC:=(-R * SB + S * RB) * DIV1;
        B:=B + DELB;
        C:=C + DELC;
        DB:=Abs(DELB) - EPS * Abs(B);
        DC:=Abs(DELC) - EPS * Abs(C);
        If (DB<=0) Or (Abs(B)<TINY) And (DC<=0) Or (Abs(C)<TINY) Then
            Exit;
    end; 
    ShowMessage('too many iterations in QROOT');
end;

end.