unit2.pas

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

unit Unit2;

interface
uses
  Windows, Messages, SysUtils, Classes, Graphics, Controls,unit1, Forms, Dialogs;
procedure QROMO(A, B:real;var SS:real; PICK:string);

implementation
Function SQU(X, BB:real):real;
begin
    SQU:=2 * X * FUNC(BB - Sqr(X));
end;
procedure MIDSQU(AA, BB:real; var S:real; N:integer);
var
    A,B,TNM,X,DEL,DDEL,SUM:real; IT,J:integer;
begin
    B:=Sqrt(BB - AA);
    A:=0;
    If N = 1 Then
    begin
        S:=(B - A) * SQU(0.5 * (A + B), BB);
        IT:=1;
    end
    Else
    begin
        IT:=Trunc(Exp((N - 2)*Ln(3)));
        TNM:=IT;
        DEL:=(B - A) / (3 * TNM);
        DDEL:=DEL + DEL;
        X:=A + 0.5 * DEL;
        Sum:=0;
        For J:=1 To IT do
        begin
            Sum:=Sum + SQU(X, BB);
            X:=X + DDEL;
            Sum:=Sum + SQU(X, BB);
            X:=X + DEL;
        end;
        S:=(S + (B - A) * Sum / TNM) / 3;
    end;
end;

Function SQL(X, AA:real):real;
begin
    SQL:=2 * X * FUNC(AA + Sqr(X));
end;
procedure MIDSQL(AA, BB:real; var S:real; N:integer);
var
    A,B,TNM,X,DEL,DDEL,SUM:real;  IT,J:integer;
begin
    B:=Sqrt(BB - AA);
    A:=0;
    If N = 1 Then
    begin
        S:=(B - A) * SQL(0.5 * (A + B), AA);
        IT:=1;
    end
    Else
    begin
        IT:=Trunc(Exp((N - 2)*Ln(3)));
        TNM:=IT;
        DEL:=(B - A) / (3 * TNM);
        DDEL:=DEL + DEL;
        X:=A + 0.5 * DEL;
        Sum:=0;
        For J:=1 To IT do
        begin
            Sum:=Sum + SQL(X, AA);
            X:=X + DDEL;
            Sum:=Sum + SQL(X, AA);
            X:=X + DEL;
        end;
        S:=(S + (B - A) * Sum / TNM) / 3;
    end;
end;

Function INF(X:real):real;
begin
    INF:=FUNC(1 / X) / Sqr(X);
End;
procedure MIDINF(AA, BB:real; var S:real; N:integer);
var
    A,B,TNM,X,DEL,DDEL,SUM:real;  J:integer;
begin
    B:=1 / AA;
    A:=1 / BB;
    If N = 1 Then
    begin
        S:=(B - A) * INF(0.5 * (A + B));
        IT:=1;
    end
    Else
    begin
        IT:=Round(EXP((N-2)*Ln(3)));
        TNM:=IT;
        DEL:=(B - A) / (3 * TNM);
        DDEL:=DEL + DEL;
        X:=A + 0.5 * DEL;
        Sum:=0;
        For J:=1 To IT do
        begin
            Sum:=Sum + INF(X);
            X:=X + DDEL;
            Sum:=Sum + INF(X);
            X:=X + DEL;
        end;
        S:=(S + (B - A) * Sum / TNM) / 3;
        //IT:= 3 * IT;
    end;
end;

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

procedure POLINT(XA, YA:array of real; N:integer; X:real; var Y, DY:real);
var
    C,D:array[0..15] of real;
    DIF,DIFT,HO,HP,W,DEN: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 QROMO(A, B:real;var SS:real; PICK:string);
const
    EPS = 0.00003;    JMAX = 14;
    K = 7;
var
    S, H:array[0..Jmax+1] of real;
    J:integer;
    DSS:real;
begin
    H[1]:=1;
    For J:=1 To JMAX do
    begin
        If PICK = 'MIDPNT' Then  MIDPNT(A, B, S[J], J);
        If PICK = 'MIDINF' Then  MIDINF(A, B, S[J], J);
        If PICK = 'MIDSQL' Then  MIDSQL(A, B, S[J], J);
        If PICK = 'MIDSQU' Then  MIDSQU(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]:=H[J] / 9;
    end;
    ShowMessage('Too many steps.');
end;
end.