PDA

View Full Version : Problem with CBYS



greencurry
05-20-2008, 06:00 PM
Hi All,

It looks like the function, CBYS, gives us the wrong imaginary part when the index is high.


Example: index=30, argument=4.93052
CBYS gives (-6.093590468815760E+018, -32007.1674435620)

Comparing to Mathematica and MathLab, they give -6.09359*10^18.

Actually, the real argument should give us the real result, as well.

Does some one know how to fix this thing?

Now, I'm using IMSL Fortran Library 5.0, working with Intel Fortran Compiler

Thank you all.

pate
05-22-2008, 12:20 PM
What is the actual call you are making, all four arguments. Can you post your code here?

greencurry
05-22-2008, 01:31 PM
Thank you pate. Below is a part of my code.







INCLUDE 'link_f90_static.h'
USE CBYS_INT

INTEGER, PARAMETER :: DBL = SELECTED_REAL_KIND(p = 15, r = 307)
COMPLEX(kind = DBL), DIMENSION(-30:30, -30:30) :: inibesselYTMOut
REAL(kind = DBL) :: r, xnu
COMPLEX(kind = DBL) :: kapTMOut, z

xnu = 0d0
r = 1d0
kapTMOut = DCMPLX(4.93052)
z = DCMPLX(kapTMOut*r)

CALL DCBYS(xnu, z, 32, inibesselYTMOut)






That I have to define kapTMOut as complex number is because I would like to make it both real and complex numbers in the future.

Thank again,

greencurry

pate
05-22-2008, 02:19 PM
I need your entire code, please post your write statements as well so that I can see if you have made an error there.

-Merv

greencurry
05-22-2008, 02:56 PM
pate, thank for this quick reply. this is my entire code.







PROGRAM besselY
USE CBYS_INT
INCLUDE 'link_f90_static.h'

INTEGER, PARAMETER :: DBL = SELECTED_REAL_KIND(p = 15, r = 307)
INTEGER, PARAMETER :: m_max = 30

INTEGER :: m, istat
REAL(kind = DBL) :: r, xnu
COMPLEX(kind = DBL) :: kapTMOut, z
COMPLEX(kind = DBL), DIMENSION(1:m_max) :: inibesselYTMOut
CHARACTER *128 :: str


xnu = 0d0
r = 1d0
kapTMOut = DCMPLX(4.93052)
z = DCMPLX(kapTMOut*r)


CALL DCBYS(xnu, z, m_max, inibesselYTMOut)


str = '.\besselY.dat'
OPEN(unit = 101, file = str, status = 'replace', iostat = istat)
DO m = 1, m_max, 1
WRITE(unit = 101, fmt = *) m-1, inibesselYTMOut(m)
END DO
CLOSE(unit = 101, status = 'keep')

END PROGRAM besselY










The following is the result.







0 (-0.297432512682031,0.000000000000000E+000)
1 (0.171156155644781,-4.503120544633727E-017)
2 (0.366859734817855,2.775557561562891E-017)
3 (0.126467405785863,1.342532207640358E-016)
4 (-0.212960262583987,1.665334536937735E-016)
5 (-0.472005418154108,1.885616988600198E-016)
6 (-0.744353393242362,2.081668171172169E-016)
7 (-1.33961700999083,1.891021036866380E-016)
8 (-3.05943158693743,1.526556658859590E-016)
9 (-8.58852542090111,-1.567027305534095E-016)
10 (-28.2949601948126,-2.273571955702103E-015)
11 (-106.186223951630,-2.036678619772546E-013)
12 (-445.508391096271,-9.039660296347729E-014)
13 (-2062.38854325973,3.118783994790436E-012)
14 (-10430.0385081093,-3.393598334253291E-012)
15 (-57168.9044665453,-1.236537634471908E-010)
16 (-337417.067665917,-1.511049465157816E-010)
17 (-2132731.15455356,2.702806180141467E-009)
18 (-14369522.2422955,-8.194781019796581E-009)
19 (-102785775.346003,-2.474954216762473E-007)
20 (-777810494.765359,-5.388279664032477E-007)
21 (-6207384232.41281,6.346285926282026E-006)
22 (-52098992501.9770,-2.275644390908344E-004)
23 (-458724437009.677,-1.216903024144121E-003)
24 (-4227636989827.99,-3.964132368855539E-003)
25 (-40698511136087.6,3.164128447070624E-002)
26 (-408492630446347.,1.01820141822921)
27 (-4.267491402852963E+015,-12.3659897975174)
28 (-4.632989055393848E+016,-54.7893678987635)
29 (-5.219394483661518E+017,277.950662191111)





Comparing to Mathematica and Matlab, the real parts are correct, but the imaginary parts should go to zero.


Thanks again,
greencurry

greencurry
05-22-2008, 04:29 PM
OK. I tried to do many experiments, and then ended up that

the errors occur for small imaginary part of the argument (less than 0.5 in my opinion) for large order.



Is there a way to fix it?

Thank you.

pate
05-23-2008, 08:03 AM
Which MatLab function are you calling?

greencurry
05-23-2008, 09:33 AM
Pate,

I used bessely in Matlab and BesselY in Mathematica.

pate
05-23-2008, 12:04 PM
I talked to a math programmer here that is familiar with the Bessel functions and even though your arguments are complex with zero as the imaginary part, he would expect the code to behave as MatLab does. There may be a problem, so I am going to file a CR on this issue, for it to be looked at, at some point by the Development team.

But for real arguments our routines are broken down separately. CBYS evaluates Bessel functions of the second kind with complex arguments, while the routine BSYS evaluates Bessel functions of the second kind with real arguments. Since you are declaring z as a complex number but only providing the real part of the number, you may be able to use BSYS as a workaround to use and just pass in the real numbers. I tried the BSYS example in our manual in Matlab and it evaluates to -3.189 as I would expect. So for real numbers as arguments I think BSYS is the correct choice of routine to use in the is situation. In the meantime, I'll log this for someone to look at.

greencurry
05-23-2008, 03:00 PM
pate,

I really thank you for the very quick reply. I would like to conclude problems I have found to you.

CBJS and CBYS make numerical errors when assigned complex arguments are pure real or pure imaginary numbers.

The real reason I wanna use CBJS and CBYS is my research has to use pure imaginary part. Actually, BSJNS works very well for me because I consider integer order and complex argument's pure imaginary part. However, there is no this kind of BSJNS for Bessel Y function.

I really thank you again,
greencurry