Fortran有没有可以调用半整数阶的贝塞尔函数库？

ii08842 · 发表于 2017-1-11 21:00:39

用IVF2013调用MKL库的贝塞尔函数，结果看帮助文档中阶数n只能为整型变量。有没有编译了半整数阶的库函数可以直接调用的？

li913 · 发表于 2017-1-12 19:31:20

半整数阶贝塞尔是初等函数，利用贝塞尔性质，可以自己写任意阶的代码。

zdhuang · 发表于 2017-6-25 23:01:22

JYV -- Evaluate a sequence of Bessel functions of the first and
second kinds and their derivatives with arbitrary real orders
and real arguments.
书名《COMPUTATION OF SPECIAL FUNCTIONS》
作者 Shanjie Zhang and Jianming Jin
以下Fortran77代码来自这本书，仅供参考。

[Fortran] 纯文本查看 复制代码

         SUBROUTINE JYV(V,X,VM,BJ,DJ,BY,DY)
C
C       =======================================================
C       Purpose: Compute Bessel functions Jv(x) and Yv(x)
C                and their derivatives
C       Input :  x --- Argument of Jv(x) and Yv(x)
C                v --- Order of Jv(x) and Yv(x)
C                      ( v = n+v0, 0 ף v0 < 1, n = 0,1,2,... )
C       Output:  BJ(n) --- Jn+v0(x)
C                DJ(n) --- Jn+v0'(x)
C                BY(n) --- Yn+v0(x)
C                DY(n) --- Yn+v0'(x)
C                VM --- Highest order computed
C       Routines called:
C            (1) GAMMA for computing gamma function
C            (2) MSTA1 and MSTA2 for computing the starting
C                point for backward recurrence
C       =======================================================
C
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        DIMENSION BJ(0:*),DJ(0:*),BY(0:*),DY(0:*)
        EL=.5772156649015329D0
        PI=3.141592653589793D0
        RP2=.63661977236758D0
        X2=X*X
        N=INT(V)
        V0=V-N
        IF (X.LT.1.0D-100) THEN
           DO 10 K=0,N
              BJ(K)=0.0D0
              DJ(K)=0.0D0
              BY(K)=-1.0D+300
10            DY(K)=1.0D+300
           IF (V0.EQ.0.0) THEN
              BJ(0)=1.0D0
              DJ(1)=0.5D0
           ELSE
              DJ(0)=1.0D+300
           ENDIF
           VM=V  
           RETURN
        ENDIF
        IF (X.LE.12.0) THEN
           DO 25 L=0,1
              VL=V0+L
              BJVL=1.0D0
              R=1.0D0
              DO 15 K=1,40
                 R=-0.25D0*R*X2/(K*(K+VL))
                 BJVL=BJVL+R
                 IF (DABS(R).LT.DABS(BJVL)*1.0D-15) GO TO 20
15            CONTINUE
20            VG=1.0D0+VL
              CALL GAMMA(VG,GA)
              A=(0.5D0*X)**VL/GA
              IF (L.EQ.0) BJV0=BJVL*A
              IF (L.EQ.1) BJV1=BJVL*A
25         CONTINUE
        ELSE
           K0=11
           IF (X.GE.35.0) K0=10
           IF (X.GE.50.0) K0=8
           DO 40 J=0,1
              VV=4.0D0*(J+V0)*(J+V0)
              PX=1.0D0
              RP=1.0D0
              DO 30 K=1,K0
                 RP=-0.78125D-2*RP*(VV-(4.0*K-3.0)**2.0)*(VV-
     &              (4.0*K-1.0)**2.0)/(K*(2.0*K-1.0)*X2)
30               PX=PX+RP
              QX=1.0D0
              RQ=1.0D0
              DO 35 K=1,K0
                 RQ=-0.78125D-2*RQ*(VV-(4.0*K-1.0)**2.0)*(VV-
     &              (4.0*K+1.0)**2.0)/(K*(2.0*K+1.0)*X2)
35               QX=QX+RQ
              QX=0.125D0*(VV-1.0)*QX/X
              XK=X-(0.5D0*(J+V0)+0.25D0)*PI
              A0=DSQRT(RP2/X)
              CK=DCOS(XK)
              SK=DSIN(XK)
              IF (J.EQ.0) THEN
                 BJV0=A0*(PX*CK-QX*SK)
                 BYV0=A0*(PX*SK+QX*CK)
              ELSE IF (J.EQ.1) THEN
                 BJV1=A0*(PX*CK-QX*SK)
                 BYV1=A0*(PX*SK+QX*CK)
              ENDIF
40         CONTINUE
        ENDIF
        BJ(0)=BJV0
        BJ(1)=BJV1
        DJ(0)=V0/X*BJ(0)-BJ(1)
        DJ(1)=-(1.0D0+V0)/X*BJ(1)+BJ(0)
        IF (N.GE.2.AND.N.LE.INT(0.9*X)) THEN
           F0=BJV0
           F1=BJV1
           DO 45 K=2,N
              F=2.0D0*(K+V0-1.0D0)/X*F1-F0
              BJ(K)=F
              F0=F1
45            F1=F
        ELSE IF (N.GE.2) THEN
           M=MSTA1(X,200)
           IF (M.LT.N) THEN
              N=M
           ELSE
              M=MSTA2(X,N,15)
           ENDIF
           F2=0.0D0
           F1=1.0D-100
           DO 50 K=M,0,-1
              F=2.0D0*(V0+K+1.0D0)/X*F1-F2
              IF (K.LE.N) BJ(K)=F
              F2=F1
50            F1=F
           IF (DABS(BJV0).GT.DABS(BJV1)) THEN
               CS=BJV0/F
           ELSE
               CS=BJV1/F2
           ENDIF
           DO 55 K=0,N
55            BJ(K)=CS*BJ(K)
        ENDIF
        DO 60 K=2,N
60         DJ(K)=-(K+V0)/X*BJ(K)+BJ(K-1)
        IF (X.LE.12.0D0) THEN
           IF (V0.NE.0.0) THEN
              DO 75 L=0,1
                 VL=V0+L
                 BJVL=1.0D0
                 R=1.0D0
                 DO 65 K=1,40
                    R=-0.25D0*R*X2/(K*(K-VL))
                    BJVL=BJVL+R
                    IF (DABS(R).LT.DABS(BJVL)*1.0D-15) GO TO 70
65               CONTINUE
70               VG=1.0D0-VL
                 CALL GAMMA(VG,GB)
                 B=(2.0D0/X)**VL/GB
                 IF (L.EQ.0) BJU0=BJVL*B
                 IF (L.EQ.1) BJU1=BJVL*B
75            CONTINUE
              PV0=PI*V0
              PV1=PI*(1.0D0+V0)
              BYV0=(BJV0*DCOS(PV0)-BJU0)/DSIN(PV0)
              BYV1=(BJV1*DCOS(PV1)-BJU1)/DSIN(PV1)
           ELSE
              EC=DLOG(X/2.0D0)+EL
              CS0=0.0D0
              W0=0.0D0
              R0=1.0D0
              DO 80 K=1,30
                 W0=W0+1.0D0/K
                 R0=-0.25D0*R0/(K*K)*X2
80               CS0=CS0+R0*W0
              BYV0=RP2*(EC*BJV0-CS0)
              CS1=1.0D0
              W1=0.0D0
              R1=1.0D0
              DO 85 K=1,30
                 W1=W1+1.0D0/K
                 R1=-0.25D0*R1/(K*(K+1))*X2
85               CS1=CS1+R1*(2.0D0*W1+1.0D0/(K+1.0D0))
              BYV1=RP2*(EC*BJV1-1.0D0/X-0.25D0*X*CS1)
           ENDIF
        ENDIF
        BY(0)=BYV0
        BY(1)=BYV1
        DO 90 K=2,N
           BYVK=2.0D0*(V0+K-1.0D0)/X*BYV1-BYV0
           BY(K)=BYVK
           BYV0=BYV1
90         BYV1=BYVK
        DY(0)=V0/X*BY(0)-BY(1)
        DO 95 K=1,N
95         DY(K)=-(K+V0)/X*BY(K)+BY(K-1)
        VM=N+V0
        RETURN
        END


        SUBROUTINE GAMMA(X,GA)
C
C       ==================================================
C       Purpose: Compute gamma function ג(x)
C       Input :  x  --- Argument of ג(x)
C                       ( x is not equal to 0,-1,-2,תתת)
C       Output:  GA --- ג(x)
C       ==================================================
C
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        DIMENSION G(26)
        PI=3.141592653589793D0
        IF (X.EQ.INT(X)) THEN
           IF (X.GT.0.0D0) THEN
              GA=1.0D0
              M1=X-1
              DO 10 K=2,M1
10               GA=GA*K
           ELSE
              GA=1.0D+300
           ENDIF
        ELSE
           IF (DABS(X).GT.1.0D0) THEN
              Z=DABS(X)
              M=INT(Z)
              R=1.0D0
              DO 15 K=1,M
15               R=R*(Z-K)
              Z=Z-M
           ELSE
              Z=X
           ENDIF
           DATA G/1.0D0,0.5772156649015329D0,
     &          -0.6558780715202538D0, -0.420026350340952D-1,
     &          0.1665386113822915D0,-.421977345555443D-1,
     &          -.96219715278770D-2, .72189432466630D-2,
     &          -.11651675918591D-2, -.2152416741149D-3,
     &          .1280502823882D-3, -.201348547807D-4,
     &          -.12504934821D-5, .11330272320D-5,
     &          -.2056338417D-6, .61160950D-8,
     &          .50020075D-8, -.11812746D-8,
     &          .1043427D-9, .77823D-11,
     &          -.36968D-11, .51D-12,
     &          -.206D-13, -.54D-14, .14D-14, .1D-15/
           GR=G(26)
           DO 20 K=25,1,-1
20            GR=GR*Z+G(K)
           GA=1.0D0/(GR*Z)
           IF (DABS(X).GT.1.0D0) THEN
              GA=GA*R
              IF (X.LT.0.0D0) GA=-PI/(X*GA*DSIN(PI*X))
           ENDIF
        ENDIF
        RETURN
        END


        INTEGER FUNCTION MSTA1(X,MP)
C
C       ===================================================
C       Purpose: Determine the starting point for backward  
C                recurrence such that the magnitude of    
C                Jn(x) at that point is about 10^(-MP)
C       Input :  x     --- Argument of Jn(x)
C                MP    --- Value of magnitude
C       Output:  MSTA1 --- Starting point   
C       ===================================================
C
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        A0=DABS(X)
        N0=INT(1.1*A0)+1
        F0=ENVJ(N0,A0)-MP
        N1=N0+5
        F1=ENVJ(N1,A0)-MP
        DO 10 IT=1,20             
           NN=N1-(N1-N0)/(1.0D0-F0/F1)                  
           F=ENVJ(NN,A0)-MP
           IF(ABS(NN-N1).LT.1) GO TO 20
           N0=N1
           F0=F1
           N1=NN
 10        F1=F
 20     MSTA1=NN
        RETURN
        END


        INTEGER FUNCTION MSTA2(X,N,MP)
C
C       ===================================================
C       Purpose: Determine the starting point for backward
C                recurrence such that all Jn(x) has MP
C                significant digits
C       Input :  x  --- Argument of Jn(x)
C                n  --- Order of Jn(x)
C                MP --- Significant digit
C       Output:  MSTA2 --- Starting point
C       ===================================================
C
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        A0=DABS(X)
        HMP=0.5D0*MP
        EJN=ENVJ(N,A0)
        IF (EJN.LE.HMP) THEN
           OBJ=MP
           N0=INT(1.1*A0)
        ELSE
           OBJ=HMP+EJN
           N0=N
        ENDIF
        F0=ENVJ(N0,A0)-OBJ
        N1=N0+5
        F1=ENVJ(N1,A0)-OBJ
        DO 10 IT=1,20
           NN=N1-(N1-N0)/(1.0D0-F0/F1)
           F=ENVJ(NN,A0)-OBJ
           IF (ABS(NN-N1).LT.1) GO TO 20
           N0=N1
           F0=F1
           N1=NN
10         F1=F
20      MSTA2=NN+10
        RETURN
        END

        REAL*8 FUNCTION ENVJ(N,X)
        DOUBLE PRECISION X
        ENVJ=0.5D0*DLOG10(6.28D0*N)-N*DLOG10(1.36D0*X/N)
        RETURN
        END

[特殊函数] Fortran有没有可以调用半整数阶的贝塞尔函数库？