Strange result with Bessel functions

[Arjen_Markus]

Hello,

in response to this post in comp.lang.fortran (http://groups.google.nl/group/comp.lang.fortran/browse_frm/thread/451dd609b24f2585?hl=nl#), I tried the following program using Intel Fortran 11.0:

! bessel.f90 --

! Posting by "steve" on clf

!

program testjn

integer, parameter :: dp = kind(1.d0)

real(dp) z

z = 2.4048255576957729_dp

print *, besjn(3,z)

do i = 1,10

z = 2.1 + i/10.0

print *, z, besjn(1,z), besjn(2,z), besjn(3,z)

enddo

end program testjn

The results were astonishing (this is the result on Linux):

6.7595576E-33

2.19999980926514 0.0000000E+00 0.0000000E+00 0.0000000E+00

2.29999995231628 -9.7165255E-11 -8.8399815E-11 9.7165255E-11

2.39999985694885 5.4210109E-20 0.0000000E+00 0.0000000E+00

2.50000000000000 0.0000000E+00 0.0000000E+00 0.0000000E+00

2.59999990463257 -0.5767248 0.3528340 -0.1289432

2.69999980926514 0.0000000E+00 0.0000000E+00 0.0000000E+00

2.79999995231628 -9.7165255E-11 -8.8399815E-11 9.7165255E-11

2.89999985694885 5.4210109E-20 0.0000000E+00 0.0000000E+00

3.00000000000000 0.0000000E+00 0.0000000E+00 0.0000000E+00

3.09999990463257 -0.5767248 0.3528340 -0.1289432

As the Bessel functions are smooth, oscillating functions, it is clear

this is incorrect.

If this problem has been solved in Intel Fortran 11.1, then I apologise for the noise.

Regards,

Arjen

[IDZ_A_Intel]

Arjen,
On Windows Vista I changed your code to fix compiler errors :

program testjn
use ifport
integer, parameter :: dp = kind(1.d0)
real(dp) z
integer i

z = 2.4048255576957729_dp

print *, dbesjn(3,z)
do i = 1,10
z = 2.1 + i/10.0
print *, z, dbesjn(1,z), dbesjn(2,z), dbesjn(3,z)
enddo
end program testjn

and got the following using V11.1.060

-Infinity

2.19999980926514 0.555963076969823
0.395058649918586 0.162325439701542

2.29999995231628 0.539872541148635
0.413914583151530 0.179978922734544

2.39999985694885 0.520185298828659
0.430980017151422 0.198114772571065

2.50000000000000 0.497094102464274
0.446059058439617 0.216600391039114

2.59999990463257 0.470818293019128
0.458972840487571 0.235293795169479

2.69999980926514 0.441601437484153
0.469561484839453 0.254045255864694

2.79999995231628 0.409709262652801
0.477685492135148 0.272698594875890

2.89999985694885 0.375427532420198
0.483227044458551 0.291092561780032

3.00000000000000 0.339058958525937
0.486091260585891 0.309062722255251

3.09999990463257 0.300921170211927
0.486207015384473 0.326442739906775

Les

[Arjen_Markus]

Hi Les,

I thought besjnwas the generic name. I tried your changes with Intel Fortran 11.0 and got all the results I expected to see - including the correct value for besjn(3,2.4048...).

I have seen the -Infinity value you get with an old version of gfortran under Linux. So there is something funny going on.

Regards,

Arjen

[Steven_L_Intel1]

The Bessel functions are not generic - at least not in the current version. Please read the Intel documentation.

[Melvin_Robinson]

How can we set the compiler options or modify the program to get the above program to give the correct values? Just wanted to know because I'm just learning about Bessel functions in an emag class that I am taking.

[Steven_L_Intel1]

Reply 1 has the corrected code.

[Arjen_Markus]

Yes, they are not part of the Fortran 95 or 2003 standard, so that figures.

I have grown so used to such functions being generic, that I did not even think about that ;)

Regards,

Arjen

[Steven_L_Intel1]

Fortran 2008 has Bessel functions (with different names from what C uses.)