unit2.pas

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

unit Unit2;

interface
uses
  Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms,unit1, Dialogs;
procedure POLINT(XA, YA:array of real; N:integer; X:real; var Y, DY:real);
procedure TRAPZD(A, B:real;var S:real; N:integer);
procedure QROMB(A, B:real;var SS:real);

implementation
procedure POLINT(XA, YA:array of real; N:integer; X:real; var Y, DY:real);
var
    C,D:array[1..10] of real;
    DIF,DIFT,HO,HP,W,DEN,te:real;
    NS,I,M:integer;
begin
    NS:=1;
    DIF:=Abs(X - XA[1]);
    For I:=1 To N do
      begin
        DIFT:=Abs(X - XA[I]);
        If DIFT < DIF Then
          begin
            NS:=I;
            DIF:=DIFT;
          end;
        C[I]:=YA[I];
        D[I]:=YA[I];
      end;
    Y:=YA[NS];
    NS:=NS - 1;
    For M:=1 To N - 1 do
      begin
        For I:=1 To N - M do
          begin
            HO:=XA[I] - X;
            HP:=XA[I + M] - X;
            W:=C[I + 1] - D[I];
            DEN:=HO - HP;
            If DEN = 0 Then
              begin
                ShowMessage('PAUSE');
                Exit;
              end;
            DEN:=W / DEN;
            D[I]:=HP * DEN;
            C[I]:=HO * DEN;
          end;
        If 2 * NS < N - M Then
            DY:=C[NS + 1]
        Else
          begin
            DY:=D[NS];
            NS:=NS - 1;
          end;
        Y:=Y + DY;
      end;
end;

procedure TRAPZD(A, B:real;var S:real; N:integer);
var
    J,IT:integer;
    TNM,DEL,X,SUM:real;
begin
    If N = 1 Then
    begin
        S:=0.5 * (B - A) * (FUNC(A) + FUNC(B));
        IT:=1;
    end
    Else
    begin
        IT:=Trunc(exp((N - 2)*Ln(2)));
        TNM:=IT;
        DEL:=(B - A) / TNM;
        X:=A + 0.5 * DEL;
        Sum:=0;
        For J:=1 To IT do
        begin
            Sum:=Sum + FUNC(X);
            X:=X + DEL;
        end;
        S:=0.5 * (S + (B - A) * Sum / TNM);
    end;
end;

procedure QROMB(A, B:real; var SS:real);
CONST
    EPS = 0.000001;    JMAX = 20;
VAR
    S,H:ARRAY[0..21] of real;
    JMAXP,J,K,KM:integer;  DSS:real;
begin
    JMAXP:=JMAX + 1;
    K:=5;
    KM:=K - 1;
    H[1]:=1;
    For J:=1 To JMAX do
    begin
        TRAPZD(A, B, S[J], J);
        If J >= K Then
        begin
            POLINT(H, S, K, 0, SS, DSS);
            If (AbS(DSS) < (EPS * Abs(SS))) Then Exit;
        end;
        S[J + 1]:=S[J];
        H[J + 1]:=0.25 * H[J];
    end;
    ShowMessage('Too many steps.');
end;
end.