unit2.pas

数值分析常用算法编程

unit Unit2;

interface
uses
  Windows, Messages, SysUtils, Classes, Graphics, UNIT1,Controls, Forms, Dialogs;
Function BESSJ0(X:real):real;
Procedure ZBRAK(X1, X2:real; N:integer;
                        var XB1, XB2:array of real;var NB:integer);
Function ZBRENT(X1, X2, TOL:real):real;

implementation
Function BESSJ0(X:real):real;
var
   AAA,BBB,CCC,Y,AX,Z,DDD,EEE,XX:real;
const
  P1=1;                    P2=-0.001098628627;
  P3=0.2734510407e-4;      P4=-0.2073370639e-5;
  P5=2.093887211E-07;
  Q1=-0.1562499995e-1;     Q2=0.1430488765e-3;
  Q3=-0.6911147651e-5;     Q4=7.621095161E-07;
  Q5=-9.34945152E-08;
  R1=57568490574;          R2=-13362590354;
  R3=651619640.7;          R4=-11214424.18;
  R5=77392.33017;          R6=-184.9052456;
  S1=57568490411;          S2=1029532985;
  S3=9494680.718;          S4=59272.64853;
  S5=267.8532712;          S6=1;
begin
  If Abs(X) < 8 Then
    begin
      Y:=X * X;
      BBB:=Y* (R4+ Y* (R5+ Y* R6));
      AAA:=R1+ Y* (R2+ Y* (R3+ BBB));
      CCC:=Y* (S3+ Y* (S4+ Y* (S5+ Y* S6)));
      BESSJ0:= AAA / (S1+ Y* (S2+ CCC));
    end
  Else
    begin
      AX:=Abs(X);
      Z:=8/ AX;
      Y:=Z* Z;
      XX:= AX- 0.785398164;
      CCC:=Y* (P3+ Y* (P4+ Y* P5));
      AAA:=P1+ Y* (P2+ CCC);
      DDD:=Y* (Q3+ Y* (Q4+ Y* Q5));
      EEE:=Z* Sin(XX) * (Q1+ Y* (Q2+ DDD));
      BESSJ0:= Sqrt(0.636619772 / AX) * (Cos(XX) * AAA- EEE);
    End;
End;

Procedure ZBRAK(X1, X2:real; N:integer;
                        var XB1, XB2:array of real;var NB:integer);
var
    NBB,I:integer;   X,DX,FP,FC:real;
begin
    NBB:=NB;
    NB:=0;
    X:=X1;
    DX:=(X2 - X1) / N;
    FP:=FUN(X);
    For I:=1 To N do
    begin
        X:=X + DX;
        FC:=FUN(X);
        If FC * FP < 0  Then
        begin
            NB:=NB + 1;
            XB1[NB]:=X - DX;
            XB2[NB]:=X;
        end;
        FP:=FC;
        If NBB = NB Then Exit;
    end; 
end;

Function ZBRENT(X1, X2, TOL:real):real;
var
    ITMAX,ITER:integer;  ZZ,S,AAA,EPS,FA,FB,FC,C,TOL1,XM,P,Q,R,D,E,A,B:real;
begin
    ITMAX:=100;
    EPS:=0.00000003;
    A:=X1;
    B:=X2;
    FA:=FUN(A);
    FB:=FUN(B);
    If FB * FA > 0  Then
        ShowMessage('Root must be bracketed for ZBRENT');
    FC:=FB;
    For ITER:=1 To ITMAX do
    begin
        If FB * FC > 0 Then
        begin
            C:=A;
            FC:=FA;
            D:=B - A;
            E:=D;
        end;
        If Abs(FC) < Abs(FB) Then
        begin
            A:=B;
            B:=C;
            C:=A;
            FA:=FB;
            FB:=FC;
            FC:=FA;
        end;
        TOL1:=2 * EPS * Abs(B) + 0.5 * TOL;
        XM:=0.5 * (C - B);
        If (Abs(XM) <= TOL1) Or (FB = 0) Then
        begin
            ZBRENT:=B;
            Exit;
        end;
        If (Abs(E) >= TOL1) And (Abs(FA) > Abs(FB)) Then
        begin
            S:=FB / FA;
            If A = C Then
            begin
                P:=2 * XM * S;
                Q:=1 - S;
            end
            Else
            begin
                Q:=FA / FC;
                R:=FB / FC;
                P:=S * (2 * XM * Q * (Q - R) - (B - A) * (R - 1 ));
                Q:=(Q - 1 ) * (R - 1 ) * (S - 1 );
            end;
            If P > 0  Then Q:=-Q;
            P:=Abs(P);
            If 3 * XM * Q - Abs(TOL1*Q) < Abs(E*Q) Then
                AAA:=3  * XM * Q - Abs(TOL1 * Q)
            Else
                AAA:=Abs(E * Q);
            If 2 * P < AAA Then
            begin
                E:=D;
                D:=P / Q;
            end
            Else
            begin
                D:=XM;
                E:=D;
            end;
        end
        Else
        begin
            D:=XM;
            E:=D;
        end;
        A:=B;
        FA:=FB;
        If Abs(D) > TOL1 Then
            B:=B + D
        Else
        Begin
            If XM>=0 then
                ZZ:=1
            Else
                ZZ:=-1;
            B:=B + Abs(TOL1) * ZZ;
        end;
        FB:=FUN(B);
    end; 
    ShowMessage('ZBRENT exceeding maximum iterations.');
    ZBRENT:=B;
end;
end.