unit2.pas

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

unit Unit2;

interface
uses
  Windows, Messages, SysUtils, Classes, Graphics, unit1,Controls, Forms, Dialogs;
Procedure LINMIN(var P,XI:array of real; N:integer;var FRET:real);
Function F1DIM(X:real):real;

implementation
Function F1DIM(X:real):real;
var
    XT:array[0..50] of real;
    J:integer;
begin
    For J:=1 To NCOM do
        XT[J]:=PCOM[J] + X * XICOM[J];
    F1DIM:=FUNC2(XT, NCOM);
end;

Procedure MNBRAK(var AX, BX, CX, FA, FB, FC:real);
const
    GOLD = 1.618034;  GLIMIT = 100;  TINY = 1E-20;
var 
    DUM,R,Q,U,ULIM,FU:real;   DONE:boolean;
begin
    FA:=FUNC(AX);
    FB:=FUNC(BX);
    If FB > FA Then
    begin
        DUM:=AX;
        AX:=BX;
        BX:=DUM;
        DUM:=FB;
        FB:=FA;
        FA:=DUM;
    end;
    CX:=BX + GOLD * (BX - AX);
    FC:=FUNC(CX);
    repeat
        If FB < FC Then break;
        DONE:=true;
        R:=(BX - AX) * (FB - FC);
        Q:=(BX - CX) * (FB - FA);
        DUM:=Q - R;
        If Abs(DUM) < TINY Then DUM:=TINY;
        U:=BX - ((BX - CX) * Q - (BX - AX) * R) / (2 * DUM);
        ULIM:=BX + GLIMIT * (CX - BX);
        If (BX - U) * (U - CX) > 0  Then
        begin
            FU:=FUNC(U);
            If FU < FC Then
            begin
                AX:=BX;
                FA:=FB;
                BX:=U;
                FB:=FU;
                Exit;
            end
            Else If FU > FB Then
            begin
                CX:=U;
                FC:=FU;
                Exit;
            end;
            U:=CX + GOLD * (CX - BX);
            FU:=FUNC(U);
        end
        Else If (CX - U) * (U - ULIM) > 0  Then
        begin
            FU:=FUNC(U);
            If FU < FC Then
            begin
                BX:=CX;
                CX:=U;
                U:=CX + GOLD * (CX - BX);
                FB:=FC;
                FC:=FU;
                FU:=FUNC(U);
            end;
        end
        Else If (U - ULIM) * (ULIM - CX) >= 0  Then
        begin
            U:=ULIM;
            FU:=FUNC(U);
        end
        Else
        begin
            U:=CX + GOLD * (CX - BX);
            FU:=FUNC(U);
        end;
        If DONE Then
        begin
            AX:=BX;
            BX:=CX;
            CX:=U;
            FA:=FB;
            FB:=FC;
            FC:=FU;
        end
        Else
            DONE:=false;
    until not DONE
end;

Function BRENT(AX, BX, CX, TOL:real;var XMIN:real):real;
label 1,2,3;
const
    ITMAX = 100;  CGOLD = 0.381966;  ZEPS = 0.1e-9;
var
    A,B,D,V,X,W,E,FV1,XM,TOL1,TOL2,FX,FW,P,Q,R,ETEMP,DUM,U,ZZ,FU:real;
    ITER,I:integer;
begin
    A:=AX;
    If CX < AX Then A:=CX;
    B:=AX;
    If CX > AX Then B:=CX;
    V:=BX;
    W:=V;
    X:=V;
    E:=0;
    FX:=FUNC(X);
    FV1:=FX;
    FW:=FX;
    For ITER:=1 To ITMAX do
    begin
        XM:=0.5 * (A + B);
        TOL1:=TOL * Abs(X) + ZEPS;
        TOL2:=2 * TOL1;
        If Abs(X - XM) <= TOL2 - 0.5 * (B - A) Then goto 3;
        If Abs(E) > TOL1 Then
        begin
            R:=(X - W) * (FX - FV1);
            Q:=(X - V) * (FX - FW);
            P:=(X - V) * Q - (X - W) * R;
            Q:=2  * (Q - R);
            If Q > 0  Then P:=-P;
            Q:=Abs(Q);
            ETEMP:=E;
            E:=D;
            DUM:=Abs(0.5 * Q * ETEMP);
            If (Abs(P)>=DUM) or (P<=Q*(B - X)) or (P>=Q*(B - X)) Then
                GOTO 1;
            D:=P / Q;
            U:=X + D;
            If (U - A < TOL2) Or (B - U < TOL2) Then
            begin
              if XM-X>=0 then
                ZZ:=1
              else
                ZZ:=-1;
              D:=Abs(TOL1) * ZZ;
            end;
            goto 2;
        end;
1:      If X >= XM Then
          E:=A - X
        Else
          E:=B - X;
        D:=CGOLD * E;
2:      If Abs(D) >= TOL1 Then
            U:=X + D
        Else 
        begin
            if D >= 0 THEN
                ZZ:=1
            else
                ZZ:=-1;
            U:=X + Abs(TOL1) * ZZ;
        end;
        FU:=FUNC(U);
        If FU <= FX Then
        begin
            If U >= X Then
                A:=X
            Else
                B:=X;
            V:=W;
            FV1:=FW;
            W:=X;
            FW:=FX;
            X:=U;
            FX:=FU;
        end
        else
        begin
            If U < X Then
                A:=U
            else
                B:=U;
            If (FU <= FW) Or (W = X) Then
            begin
                V:=W;
                FV1:=FW;
                W:=U;
                FW:=FU;
            end
            Else If (FU <= FV1) Or (V = X) Or (V = W) Then
            begin
                V:=U;
                FV1:=FU;
            end;
        end;
    end; 
    If ITER > ITMAX Then ShowMessage('Brent exceed maximum iterations.');
3:  XMIN:=X;
    BRENT:=FX;
end;

Procedure LINMIN(var P,XI:array of real; N:integer;var FRET:real);
var
    TOL,AX,BX,XX,FA,FB,FX,XMIN:real; J:integer;
begin
    TOL:=0.0001;
    NCOM:=N;
    For J:=1 To N do
    begin
        PCOM[J]:=P[J];
        XICOM[J]:=XI[J];
    end;
    AX:=0 ;
    XX:=1;
    MNBRAK(AX, XX, BX, FA, FX, FB);
    FRET:=BRENT(AX, XX, BX, TOL, XMIN);
    For J:=1 To N do
    begin
        XI[J]:=XMIN * XI[J];
        P[J]:=P[J] + XI[J];
    end;
end;

end.