主题:请教fortran中如何计算第二类修正贝塞尔函数?
qirenqiqi [专家分:0] 发布于 2007-06-14 21:45:00
我现在要计算Modified Bessel functions of the second kind,也就叫第二类虚宗量贝塞尔函数,或者第二类修正贝塞尔函数。
请问fortran自带的函数库中可以直接调用吗?应该use 什么?(用什么模板)
不能直接调用的话,需要安装什么?以及怎么进行安装和后面的操作?
回复列表 (共11个回复)
沙发
f2003 [专家分:7960] 发布于 2007-06-15 01:02:00
intel编译器带bessel函数,放在模块里,要use一下。看他的pdf说明书就知道了,该怎么办说的清清楚楚。
板凳
qirenqiqi [专家分:0] 发布于 2007-06-15 08:26:00
但是好像没有我说的那种的,它只有计算第一类和第二类贝塞尔函数的语句。
3 楼
qirenqiqi [专家分:0] 发布于 2007-06-15 08:33:00
还有,“pdf说明书”是?是Fortran 90 Handbook吗?里面好像没有。
4 楼
zhangjie9u [专家分:20] 发布于 2007-06-15 11:42:00
你指的是不是这个函数,在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 楼
asymptotic [专家分:16630] 发布于 2007-06-16 21:26:00
若是只是想用用,则调用现成的库函数就可以了,上面的几位网友给你推荐的都不错;若是想学习如何写代码, 推荐你参考 computation of special functions, Shan-Jie Zhang, 好像是南京大学的老师。我这几天好好的研究了一下所有 Bessel 函数的代码;还是蛮有长进的。
6 楼
qirenqiqi [专家分:0] 发布于 2007-06-20 15:29:00
现成的库函数在哪里怎么调用呀?
我要计算的是第二类虚宗量贝塞尔函数的值
7 楼
asymptotic [专家分:16630] 发布于 2007-06-21 21:57:00
我不知道你是什么原因不理解上面各位网友的帮助,是因为你的编译器没有提供 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 楼
qirenqiqi [专家分:0] 发布于 2007-06-30 16:05:00
我肯定是你说的第二种了,很菜!
我只会用“F1”,可是,在帮助里连dbsk0都查不到,所以,就是菜到看不懂大家在说什么。[em8][em8]需要自我反省!
你的程序我拿来运行了一下,跟徐士良书中计算的差不多,前8位有效数字都一样,他的是双精度的,dbsk0也是双精度吗?但是后面的数字不同。
9 楼
qirenqiqi [专家分:0] 发布于 2007-06-30 16:07:00
是怎么知道用dbsk0的呀?是通过1楼说的什么pdf说明书吗?那是啥呀?
10 楼
qirenqiqi [专家分:0] 发布于 2007-06-30 16:46:00
开始——〉Visual Fortran 6.0——〉IMSL Fortran 90 MP Library 3.0 Help
[em8]
我来回复