[其他行业算法] 蚂蚁的问题

!--------------------------------------------------------------------------------
PROGRAM SOLVE
    IMPLICIT NONE
    INTEGER, PARAMETER::DP=KIND(1.0_16)
    INTEGER,PARAMETER::M=8
    REAL(KIND=DP),PARAMETER::EPS=5.Q-34
    REAL(KIND=DP)::A,B,ALPHA,PP,INTGVAL,C,P,D,X
    REAL(KIND=DP),EXTERNAL::YOUR_FUNC
    INTEGER::I
    P=0.345_DP
    PP=1.0_DP-P
    A=0.0_DP
    B=100.0
    ALPHA=1.5Q0
!这一部分确定积分上限的大体范围
    DO I=1000,1,-1
        D=YOUR_FUNC(B,ALPHA)
        IF(D>1.E-1)THEN
            EXIT
        ELSE
            B=B-1.0
        END IF
    END DO

    WRITE(*,FMT="(A,E41.34)")"B = ", B

    X=B
!这一部分用牛顿法精细求解积分区间上限
    DO I=1,200000
        CALL PHI_GAMMA(A,X,M,ALPHA,INTGVAL)
        C=INTGVAL-PP
        D=YOUR_FUNC(X,ALPHA)

        IF(ABS(C)>D)THEN
            IF(C<0.0_DP)THEN
                X=X+C
            ELSE
                X=X-C
            END IF
        ELSE
            X=X-C/D
        ENDIF

        WRITE(*,FMT="(A,I2.2,2X,A,E41.34,2X,A,E41.34,2X,A,E41.34)")"I = ",I,"C = ",C,"X = ",X,"FUNVAL = ",D

        IF(ABS(C)<EPS.OR.X<A)EXIT
        !READ(*,*)
    END DO
    !The precision of X is 4--5digits
END PROGRAM
!--------------------------------------------------------------------------------

SUBROUTINE PHI_GAMMA(A,B,M,ALPHA,INTGVAL)
    IMPLICIT NONE
    INTEGER, PARAMETER::DP=KIND(1.0_16)

    REAL(KIND=DP),INTENT(IN)::A,B

    REAL(KIND=DP),INTENT(IN)::ALPHA
    REAL(KIND=DP),INTENT(OUT)::INTGVAL
    INTEGER,INTENT(IN)::M

    REAL(KIND=DP),EXTERNAL::YOUR_FUNC
    REAL(KIND=DP),EXTERNAL::GUALEGINTG12


    ! A1   B1/A2 B2/A3  B3/A4...BM-1/AM  BM
    ! A|-----|-----|------|.......|-----|B

    REAL(KIND=DP)::AX(M)
    REAL(KIND=DP)::BX(M)
    REAL(KIND=DP)::DX

    INTEGER::I

    DX=(B-A)/REAL(M,KIND=DP)
    INTGVAL=0.0_DP
    DO I=1,M
        AX(I)=A+REAL(I-1,KIND=DP)*DX
        BX(I)=AX(I)+DX
        INTGVAL=INTGVAL+GUALEGINTG12(YOUR_FUNC,ALPHA,AX(I),BX(I))
    ENDDO
    RETURN
END SUBROUTINE PHI_GAMMA
!--------------------------------------------------------------------------------
FUNCTION YOUR_FUNC(X,ALPHA)
    IMPLICIT NONE
    INTEGER, PARAMETER::DP=KIND(1.0_16)
    REAL(KIND=DP)::YOUR_FUNC
    REAL(KIND=DP),INTENT(IN)::X,ALPHA
    REAL(KIND=DP)::ALF1,TEMP
    ALF1=ALPHA-1.0_DP
    TEMP=GAMMA(ALPHA)
    !................
    !FOR EXAMPLE
    !YOUR_FUNC=SIN(X)+ALPHA !OR YOUR_FUNC=X*X+ALPHA
    YOUR_FUNC=((X**ALF1)*EXP(-X))/TEMP
    RETURN
END FUNCTION YOUR_FUNC
!--------------------------------------------------------------------------------

FUNCTION GUALEGINTG12(YOUR_FUNC,ALPHA,A,B)
    IMPLICIT NONE
    INTEGER, PARAMETER::DP=KIND(1.0_16)
    REAL(KIND=DP)::GUALEGINTG12
    REAL(KIND=DP),INTENT(IN)::A,B
    REAL(KIND=DP)::YOUR_FUNC,ALPHA
    EXTERNAL::YOUR_FUNC
    REAL(KIND=DP)::RSLT
    REAL(KIND=DP)::X,APB2,ASB2
    INTEGER::I
!--------------------------------------------------------------------------------
!12-POINT GAUSS-LEGENDRE WEIGHTS AND ABSCISSAS
!--------------------------------------------------------------------------------
    !GAUSS WEIGHTS
    REAL(KIND=DP)::WEIGHTS(12)=(/&
     0.4717533638651182719461596148501741D-01,&
     0.1069393259953184309602547181939961D+00,&
     0.1600783285433462263346525295433591D+00,&
     0.2031674267230659217490644558097984D+00,&
     0.2334925365383548087608498989248781D+00,&
     0.2491470458134027850005624360429511D+00,&
     0.2491470458134027850005624360429511D+00,&
     0.2334925365383548087608498989248781D+00,&
     0.2031674267230659217490644558097984D+00,&
     0.1600783285433462263346525295433591D+00,&
     0.1069393259953184309602547181939961D+00,&
     0.4717533638651182719461596148501741D-01 &
    /)
    !GAUSS NODES IN INTERVAL[-1.0, 1.0]
    REAL(KIND=DP)::NODES(12)=(/&
    -0.9815606342467192506905490901492809D+00,&
    -0.9041172563704748566784658661190962D+00,&
    -0.7699026741943046870368938332128180D+00,&
    -0.5873179542866174472967024189405343D+00,&
    -0.3678314989981801937526915366437176D+00,&
    -0.1252334085114689154724413694638531D+00,&
     0.1252334085114689154724413694638531D+00,&
     0.3678314989981801937526915366437176D+00,&
     0.5873179542866174472967024189405343D+00,&
     0.7699026741943046870368938332128180D+00,&
     0.9041172563704748566784658661190962D+00,&
     0.9815606342467192506905490901492809D+00 &
    /)
!--------------------------------------------------------------------------------
    RSLT=0.0_DP
    ASB2=0.5Q0*(B-A)
    APB2=0.5Q0*(B+A)
    DO I=1,12
        X=ASB2*NODES(I)+APB2
        RSLT=RSLT+WEIGHTS(I)*YOUR_FUNC(X,ALPHA)
    ENDDO
    RSLT=RSLT*ASB2
    GUALEGINTG12=RSLT
    RETURN
END FUNCTION GUALEGINTG12
!--------------------------------------------------------------------------------