回 帖 发 新 帖 刷新版面

主题:请教fortran中如何计算第二类修正贝塞尔函数?


我现在要计算Modified Bessel functions of the second kind,也就叫第二类虚宗量贝塞尔函数,或者第二类修正贝塞尔函数。

请问fortran自带的函数库中可以直接调用吗?应该use 什么?(用什么模板)
不能直接调用的话,需要安装什么?以及怎么进行安装和后面的操作?

回复列表 (共11个回复)

沙发

intel编译器带bessel函数,放在模块里,要use一下。看他的pdf说明书就知道了,该怎么办说的清清楚楚。

板凳

但是好像没有我说的那种的,它只有计算第一类和第二类贝塞尔函数的语句。

3 楼

还有,“pdf说明书”是?是Fortran 90 Handbook吗?里面好像没有。

4 楼

你指的是不是这个函数,在fortran的imsl库函数中有:

  
Evaluate a sequence of exponentially scaled modified Bessel functions of the third kind of fractional order.
Usage
CALL BSKES (XNU, X, NIN, BKE)
Arguments
XNU ?Fractional order of the function.   (Input) 
XNU must be less than 1.0 in absolute value.
X ?Argument for which the sequence of Bessel functions is to be evaluated.   (Input)
NIN ?Number of elements in the sequence.   (Input)
BKE ?Vector of length NIN containing the values of the function through the series.   (Output)
Comments
1.    If NIN is positive, BKE(1) contains EXP(X) times the value of the function of order XNU, BKE(2) contains EXP(X) times the value of the function of order XNU + 1,  . . . , and BKE(NIN) contains EXP(X) times the value of the function of order XNU + NIN - 1.
2.    If NIN is negative, BKE(1) contains EXP(X) times the value of the function of order XNU, BKE(2) contains EXP(X) times the value of the function of order XNU - 1,  . . . , and BKE(ABS(NIN)) contains EXP(X) times the value of the function of order XNU + NIN + 1.
Algorithm
Function BSKES evaluates exKn + k - 1(x), for k = 1,  . . . , n. For the definition of Kn(x), see BSKS. 
Currently, n is restricted to be less than 1 in absolute value. A total of |n| values is stored in the array BKE. For n positive, BKE(1) contains exKn(x), BKE(2) contains exKn + 1(x),  . . . , and BKE(N) contains exKn + n - 1(x). For n negative, BKE(1) contains exKn(x), BKE(2) contains exKn - 1(x),  . . . , and BKE(|n|) contains exKn + n + 1(x). This routine is particularly useful for calculating sequences for large x provided n ?x. (Overflow becomes a problem if n << x.) n must not be zero, and x must not be greater than zero. Moreover, |n| must be less than 1. Also, when |n| is large compared with x, |n + n| must not be so large that exKn+n(x) ?exG(|n + n|)/[2(x/2)|n + n|] overflows. 
BSKES is based on the work of Cody (1983).
Example
In this example, Kn - 1/2(2.0), n = 1,  . . . , 6 is computed and printed.
  
C                                 Declare variables
      INTEGER    NIN
      PARAMETER  (NIN=6)
C
      INTEGER    K, NOUT
      REAL       BKE(NIN), X, XNU
      EXTERNAL   BSKES, UMACH
C                                 Compute
      XNU = 0.5
      X   = 2.0
      CALL BSKES (XNU, X, NIN, BKE)
C                                 Print the results
      CALL UMACH (2, NOUT)
      DO 10  K=1, NIN
         WRITE (NOUT,99999) X, XNU+K-1, X, BKE(K)
   10 CONTINUE
99999 FORMAT (' exp(', F6.3, ') * K sub ', F6.3,
     &     ' (', F6.3, ') = ', F8.3)
      END
  
Output
  
exp( 2.000) * K sub  0.500 ( 2.000) =    0.886
exp( 2.000) * K sub  1.500 ( 2.000) =    1.329
exp( 2.000) * K sub  2.500 ( 2.000) =    2.880
exp( 2.000) * K sub  3.500 ( 2.000) =    8.530
exp( 2.000) * K sub  4.500 ( 2.000) =   32.735
exp( 2.000) * K sub  5.500 ( 2.000) =  155.837
  

 

5 楼

若是只是想用用,则调用现成的库函数就可以了,上面的几位网友给你推荐的都不错;若是想学习如何写代码, 推荐你参考 computation of special functions, Shan-Jie Zhang, 好像是南京大学的老师。我这几天好好的研究了一下所有 Bessel 函数的代码;还是蛮有长进的。

6 楼

现成的库函数在哪里怎么调用呀?
我要计算的是第二类虚宗量贝塞尔函数的值

7 楼

我不知道你是什么原因不理解上面各位网友的帮助,是因为你的编译器没有提供 IMSL 库,还是因为你太“菜鸟”。我这里有一个本人昔日写的小程序,调用 IMSL 中的 Bessel 函数, CVF6.6。

! this module is to the date type Operating System independent
module TypeKind
  implicit none
  integer(kind = kind(1)), parameter:: IP = kind(1)
  integer(kind = IP), parameter:: LP = kind(.true.)
  integer(kind = IP), parameter:: SP = kind(1.0) 
  integer(kind = IP), parameter:: DP = kind(1.d0)
end module TypeKind



program main
  implicit none

!  call Test_BSK0
  call Test_CBJS

  stop
end program main



subroutine Test_BSK0
  use TypeKind
  use imslf90
  implicit none

  real(kind = DP), external:: dbsk0
  real(kind = DP):: x = 0.5d0
  real(kind = DP):: y

  y = dbsk0(x)
  write(*, *) y

  return
end subroutine


subroutine Test_CBJS
  use TypeKind
  use imslf90
  implicit none
  ! Number of elements in the sequence 
  integer(kind = IP), parameter:: N = 4
  ! complex argument for which the sequence of Bessel functions
  ! is to be evaluated
  complex(kind = DP):: z 
  ! real argument which is the lowest order desired
  real(kind = DP):: XNU
  ! Vector of length N containing the values of the function
  ! through the series
  complex(kind = DP):: CBS(N)

  ! external subroutine / function
  external:: DCBJS, DCBYS, DCBIS, DCBKS
  integer(kind = IP):: i   ! for loop

  XNU = 0.30d0
  z = (1.2d0, 0.5d0)
  
  write(*, *) "call CBJS:"
  call DCBJS(XNU, z, N, CBS)
  do i = 1, N, 1
    write(*, *) CBS(i)
  end do

  write(*, *) "call CBYS:"
  call DCBYS(XNU, z, N, CBS)
  do i = 1, N, 1
    write(*, *) CBS(i)
  end do

  write(*, *) "call CBIS:"
  call DCBIS(XNU, z, N, CBS)
  do i = 1, N, 1
    write(*, *) CBS(i)
  end do

  write(*, *) "call CBKS:"
  call DCBKS(XNU, z, N, CBS)
  do i = 1, N, 1
    write(*, *) CBS(i)
  end do


  return
end subroutine

8 楼


我肯定是你说的第二种了,很菜!
我只会用“F1”,可是,在帮助里连dbsk0都查不到,所以,就是菜到看不懂大家在说什么。[em8][em8]需要自我反省!
你的程序我拿来运行了一下,跟徐士良书中计算的差不多,前8位有效数字都一样,他的是双精度的,dbsk0也是双精度吗?但是后面的数字不同。

9 楼

是怎么知道用dbsk0的呀?是通过1楼说的什么pdf说明书吗?那是啥呀?

10 楼

开始——〉Visual Fortran 6.0——〉IMSL Fortran 90 MP Library 3.0 Help 

[em8]

我来回复

您尚未登录,请登录后再回复。点此登录或注册