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.

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

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.

Which MatLab function are you calling?

greencurry

05-23-2008, 09:33 AM

Pate,

I used bessely in Matlab and BesselY in Mathematica.

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

Powered by vBulletin® Version 4.2.3 Copyright © 2020 vBulletin Solutions, Inc. All rights reserved.