蚂蚁的问题
以下的内容根据fortran code群聊天记录，尝试写的，不知道是否能对网名“蚂蚁”的群友有所帮助!--------------------------------------------------------------------------------
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/A3B3/A4...BM-1/AMBM
! 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
!-------------------------------------------------------------------------------- 帮助大大滴，非常感谢热心的kerb老师
页:
[1]