HP Forums / HP Calculators (and very old HP Computers) / HP Prime v / Bessel functions

Post Reply 
Bessel functions
10-30-2017, 10:21 PM
Post: #1
salvomic Offline
Senior Member
 		Posts: 1,396
Joined: Jan 2015
Bessel functions
hi everybody,
there is a good post in the HP Prime Software part of the Forum (here) about the Bessel functions of 1st kind, by Eddie W. Shore and roadrunner.
I need some help to implementate in a simple way also the 2nd kind functions, starting or not by the functions proposed by Eddie and roadrunner.
Any help will be appreciated, thank you.

Salvo
∫aL√0mic (IT9CLU) :: HP Prime 50g 41CX 71b 42s 39s 35s 12C 15C - DM42, DM41X - WP34s Prime Soft. Lib
Visit this user's website Find all posts by this user
Quote this message in a reply
10-31-2017, 01:12 AM (This post was last modified: 11-01-2017 11:48 AM by toml_12953.)
Post: #2
toml_12953 Offline
Senior Member
 		Posts: 2,200
Joined: Dec 2013
RE: Bessel functions
(10-30-2017 10:21 PM)salvomic Wrote:  hi everybody,
there is a good post in the HP Prime Software part of the Forum (here) about the Bessel functions of 1st kind, by Eddie W. Shore and roadrunner.
I need some help to implementate in a simple way also the 2nd kind functions, starting or not by the functions proposed by Eddie and roadrunner.
Any help will be appreciated, thank you.

Salvo

Here's one I converted from a C program I found at Bessel Second Kind (the Prime comes up with a different answer. The answer in the code is right - I don't know why the Prime is so far off.):

N.B. The code below changes LOG to LN which fixes the problem in the original code.

Code:
//************************************************************************
//*                                                                      *
//*    Program to calculate the second kind Bessel function of integer   *
//*    order N, for any REAL X, using the function BESSY(N,X).           *
//*                                                                      *
//* -------------------------------------------------------------------- *
//*                                                                      *
//*    SAMPLE RUN:                                                       *
//*                                                                      *
//*    (Calculate Bessel function for N:=2, X:=0.75).                    *
//*                                                                      *
//*    Second kind Bessel function of order  2 for X =  0.7500:          *
//*                                                                      *
//*         Y := -2.62974604                                             *
//*                                                                      *
//* -------------------------------------------------------------------- *
//*   Reference: From Numath Library By Tuan Dang Trong in Fortran 77.   *
//*                                                                      *
//*                               C++ Release 1.0 By J-P Moreau, Paris.  *
//*                                        (www.jpmoreau.fr)             *
//*                                  Converted to HPPL by Tom Lake       *
//***********************************************************************/

BESSJ0 (X)
BEGIN
//***********************************************************************
//      This subroutine calculates the First Kind Bessel Function of
//      order 0, for any real number X. The polynomial approximation by
//      series of Chebyshev polynomials is used for 0<X<8 and 0<8/X<1.
//      REFERENCES:
//      M.ABRAMOWITZ,I.A.STEGUN, HANDBOOK OF MATHEMATICAL FUNCTIONS, 1965.
//      C.W.CLENSHAW, NATIONAL PHYSICAL LABORATORY MATHEMATICAL TABLES,
//      VOL.5, 1962.
//************************************************************************/

    LOCAL P1:=1.0, P2:=-0.1098628627E-2, P3:=0.2734510407E-4;
    LOCAL P4:=-0.2073370639E-5, P5:= 0.2093887211E-6;
    LOCAL Q1:=-0.1562499995E-1, Q2:= 0.1430488765E-3, Q3:=-0.6911147651E-5;
    LOCAL Q4:= 0.7621095161E-6, Q5:=-0.9349451520E-7;
    LOCAL R1:= 57568490574.0, R2:=-13362590354.0, R3:=651619640.7;
    LOCAL R4:=-11214424.18, R5:= 77392.33017, R6:=-184.9052456;
    LOCAL S1:= 57568490411.0, S2:=1029532985.0, S3:=9494680.718;
    LOCAL S4:= 59272.64853, S5:=267.8532712, S6:=1.0;
    LOCAL AX,FR,FS,Z,FP1,FQ,XX,Y;
    LOCAL TMP;

    IF (X==0.0) THEN RETURN 1.0; END;
    AX := ABS(X);
    IF (AX < 8.0) THEN
        Y := X*X;
        FR := R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6))));
        FS := S1+Y*(S2+Y*(S3+Y*(S4+Y*(S5+Y*S6))));
        TMP := FR/FS;     
    ELSE 
        Z := 8./AX;
        Y := Z*Z;
        XX := AX-0.785398164;
        FP1 := P1+Y*(P2+Y*(P3+Y*(P4+Y*P5)));
        FQ := Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*Q5)));
        TMP := sqrt(0.636619772/AX)*(FP1*COS(XX)-Z*FQ*SIN(XX));
    END;
    RETURN TMP;
END;

Signu (X,Y)
BEGIN
    IF (Y<0.0) THEN 
        RETURN (-ABS(X));
    ELSE 
        RETURN (ABS(X));
    END;
END;

BESSJ1 (X)
BEGIN  
//**********************************************************************
//      This subroutine calculates the First Kind Bessel Function of
//     order 1, for any real number X. The polynomial approximation by
//      series of Chebyshev polynomials is used for 0<X<8 and 0<8/X<1.
//      REFERENCES:
//      M.ABRAMOWITZ,I.A.STEGUN, HANDBOOK OF MATHEMATICAL FUNCTIONS, 1965.
//      C.W.CLENSHAW, NATIONAL PHYSICAL LABORATORY MATHEMATICAL TABLES,
//      VOL.5, 1962.
//***********************************************************************/

    LOCAL P1:=1.0, P2:=0.183105E-2, P3:=-0.3516396496E-4, P4:=0.2457520174E-5;
    LOCAL P5:=-0.240337019E-6,  P6:=0.636619772;
    LOCAL Q1:= 0.04687499995, Q2:=-0.2002690873E-3, Q3:=0.8449199096E-5;
    LOCAL Q4:=-0.88228987E-6, Q5:= 0.105787412E-6;
    LOCAL R1:= 72362614232.0, R2:=-7895059235.0, R3:=242396853.1;
    LOCAL R4:=-2972611.439,   R5:=15704.48260,  R6:=-30.16036606;
    LOCAL S1:=144725228442.0, S2:=2300535178.0, S3:=18583304.74;
    LOCAL S4:=99447.43394,    S5:=376.9991397,  S6:=1.0;

    LOCAL AX,FR,FS,Y,Z,FP1,FQ,XX;
    LOCAL TMP;

    AX := ABS(X);
    IF (AX < 8.0) THEN
        Y := X*X;
        FR := R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6))));
        FS := S1+Y*(S2+Y*(S3+Y*(S4+Y*(S5+Y*S6))));
        TMP := X*(FR/FS);      
    ELSE
        Z := 8.0/AX;
        Y := Z*Z;
        XX := AX-2.35619491;
        FP1 := P1+Y*(P2+Y*(P3+Y*(P4+Y*P5)));
        FQ := Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*Q5)));
        TMP := sqrt(P6/AX)*(COS(XX)*FP1-Z*SIN(XX)*FQ)*Signu(S6,X);
    END;
    RETURN TMP;
END;


BESSY0(X)
BEGIN
//* --------------------------------------------------------------------
//      This subroutine calculates the Second Kind Bessel Function of
//      order 0, for any real number X. The polynomial approximation by
//      series of Chebyshev polynomials is used for 0<X<8 and 0<8/X<1.
//      REFERENCES:
//      M.ABRAMOWITZ,I.A.STEGUN, HANDBOOK OF MATHEMATICAL FUNCTIONS, 1965.
//      C.W.CLENSHAW, NATIONAL PHYSICAL LABORATORY MATHEMATICAL TABLES,
//      VOL.5, 1962.
//  --------------------------------------------------------------------- */

    LOCAL P1:= 1.0, P2:=-0.1098628627E-2, P3:=0.2734510407E-4;
    LOCAL P4:=-0.2073370639E-5, P5:= 0.2093887211E-6;
    LOCAL Q1:=-0.1562499995E-1, Q2:= 0.1430488765E-3, Q3:=-0.6911147651E-5;
    LOCAL Q4:= 0.7621095161E-6, Q5:=-0.9349451520E-7;
    LOCAL R1:=-2957821389.0, R2:=7062834065.0, R3:=-512359803.6;
    LOCAL R4:= 10879881.29,  R5:=-86327.92757, R6:=228.4622733;
    LOCAL S1:= 40076544269.0, S2:=745249964.8, S3:=7189466.438;
    LOCAL S4:= 47447.26470,   S5:=226.1030244, S6:=1.0;
    LOCAL FS,FR,Z,FP1,FQ,XX,Y;
    IF (X == 0.0) THEN RETURN -1e30; END;
    IF (X < 8.0) THEN
        Y := X*X;
        FR := R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6))));
        FS := S1+Y*(S2+Y*(S3+Y*(S4+Y*(S5+Y*S6))));
        RETURN (FR/FS+0.636619772*BESSJ0(X)*LN(X));
    ELSE 
        Z := 8.0/X;
        Y := Z*Z;
        XX := X-0.785398164;
        FP1 := P1+Y*(P2+Y*(P3+Y*(P4+Y*P5)));
        FQ := Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*Q5)));
        RETURN (sqrt(0.636619772/X)*(FP1*SIN(XX)+Z*FQ*COS(XX)));
    END;
END;

BESSY1(X)
BEGIN
//* ---------------------------------------------------------------------
//      This subroutine calculates the Second Kind Bessel Function of
//      order 1, for any real number X. The polynomial approximation by
//      series of Chebyshev polynomials is used for 0<X<8 and 0<8/X<1.
//      REFERENCES:
//      M.ABRAMOWITZ,I.A.STEGUN, HANDBOOK OF MATHEMATICAL FUNCTIONS, 1965.
//      C.W.CLENSHAW, NATIONAL PHYSICAL LABORATORY MATHEMATICAL TABLES,
//      VOL.5, 1962.
//  ---------------------------------------------------------------------- */

    LOCAL P1:= 1.0, P2:=0.183105E-2, P3:=-0.3516396496E-4;
    LOCAL P4:= 0.2457520174E-5, P5:=-0.240337019E-6;
    LOCAL Q1:= 0.04687499995, Q2:=-0.2002690873E-3, Q3:=0.8449199096E-5;
    LOCAL Q4:=-0.88228987E-6, Q5:= 0.105787412E-6;
    LOCAL R1:=-0.4900604943E13, R2:= 0.1275274390E13, R3:=-0.5153438139E11;
    LOCAL R4:= 0.7349264551E9,  R5:=-0.4237922726E7,  R6:= 0.8511937935E4;
    LOCAL S1:= 0.2499580570E14, S2:= 0.4244419664E12, S3:= 0.3733650367E10;
    LOCAL S4:= 0.2245904002E8,  S5:= 0.1020426050E6,  S6:= 0.3549632885E3, S7:=1.0;
    LOCAL FR,FS,Z,FP1,FQ,XX, Y;
    IF (X == 0.0) THEN RETURN -1e30; END;
    IF (X < 8.0) THEN
        Y := X*X;
        FR := R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6))));
        FS := S1+Y*(S2+Y*(S3+Y*(S4+Y*(S5+Y*(S6+Y*S7)))));
        RETURN (X*(FR/FS)+0.636619772*(BESSJ1(X)*LN(X)-1.0/X));      
    ELSE 
        Z := 8./X;
        Y := Z*Z;
        XX := X-2.356194491;
        FP1 := P1+Y*(P2+Y*(P3+Y*(P4+Y*P5)));
        FQ := Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*Q5)));
        RETURN (sqrt(0.636619772/X)*(SIN(XX)*FP1+Z*COS(XX)*FQ));
    END;
END;

EXPORT BESSY(N,X)
  BEGIN
//* -----------------------------------------------------------------
//      This subroutine calculates the second kind Bessel Function of
//      integer order N, for any real X. We use here the classical
//      recursive formula. 
//  ------------------------------------------------------------------ */
    LOCAL TOX,BY,BYM,BYP,J;
    IF (N == 0) THEN RETURN BESSY0(X); END;
    IF (N == 1) THEN RETURN BESSY1(X); END;
    IF (X == 0.0) THEN RETURN -1e30;  END;
    TOX := 2.0/X;
    BY  := BESSY1(X);
    BYM := BESSY0(X);
    FOR J:=1 TO N-1 DO
        BYP := J*TOX*BY-BYM;
        BYM := BY;
        BY  := BYP;
    END;
    RETURN BY;
  END;
// --------------------------------------------------------------------------
BESSYP(N,X)
BEGIN
    IF (N == 0) THEN
        RETURN (-BESSY(1,X));
    ELSE 
        IF (X == 0.0) THEN
            RETURN 1e-30;
        ELSE
            RETURN (BESSY(N-1,X)-(1.0*N/X)*BESSY(N,X));
        END;
    END;
END;

EXPORT Tbessy()
BEGIN

  LOCAL X,Y,N;

  N:=2;
  X:=0.75;

  Y := BESSY(N,X);

  PRINT("\n Second Kind Bessel Function of order " + N + ", for X= " + X + "\n\n");
  PRINT("      Y = "+CAS(format(Y,"s9"))+"\n\n");

END;

// End of file Tbessy.hpprgm
Tom L
Cui bono?
Find all posts by this user
Quote this message in a reply
10-31-2017, 08:37 AM
Post: #3
salvomic Offline
Senior Member
 		Posts: 1,396
Joined: Jan 2015
RE: Bessel functions
(10-31-2017 01:12 AM)toml_12953 Wrote:  Here's one I converted from a C program I found at Bessel Second Kind (the Prime comes up with a different answer. The answer in the code is right - I don't know why the Prime is so far off.):

...

Thank you.
In fact is very strange: VC returns -2.622... instead of -2.629...
∫aL√0mic (IT9CLU) :: HP Prime 50g 41CX 71b 42s 39s 35s 12C 15C - DM42, DM41X - WP34s Prime Soft. Lib
Visit this user's website Find all posts by this user
Quote this message in a reply
10-31-2017, 09:14 PM
Post: #4
toml_12953 Offline
Senior Member
 		Posts: 2,200
Joined: Dec 2013
RE: Bessel functions
(10-31-2017 08:37 AM)salvomic Wrote:  
(10-31-2017 01:12 AM)toml_12953 Wrote:  Here's one I converted from a C program I found at Bessel Second Kind (the Prime comes up with a different answer. The answer in the code is right - I don't know why the Prime is so far off.):

...

Thank you.
In fact is very strange: VC returns -2.622... instead of -2.629...

Problem solved! The HPPL function LOG computes the common log while in C (and BASIC) LOG returns the natural log. In the HPPL program change all the LOG to LN and the program works as expected.
Tom L
Cui bono?
Find all posts by this user
Quote this message in a reply
10-31-2017, 11:56 PM
Post: #5
salvomic Offline
Senior Member
 		Posts: 1,396
Joined: Jan 2015
RE: Bessel functions
(10-31-2017 09:14 PM)toml_12953 Wrote:  
(10-31-2017 08:37 AM)salvomic Wrote:  Thank you.
In fact is very strange: VC returns -2.622... instead of -2.629...

Problem solved! The HPPL function LOG computes the common log while in C (and BASIC) LOG returns the natural log. In the HPPL program change all the LOG to LN and the program works as expected.

well, Tom!
please, could you put a correct version in the Prime Software section of the Forum under something like "Bessel 2nd kind..."? Or, if you want could do it for you...

Salvo
∫aL√0mic (IT9CLU) :: HP Prime 50g 41CX 71b 42s 39s 35s 12C 15C - DM42, DM41X - WP34s Prime Soft. Lib
Visit this user's website Find all posts by this user
Quote this message in a reply
11-01-2017, 12:22 AM
Post: #6
toml_12953 Offline
Senior Member
 		Posts: 2,200
Joined: Dec 2013
RE: Bessel functions
(10-31-2017 11:56 PM)salvomic Wrote:  
(10-31-2017 09:14 PM)toml_12953 Wrote:  Problem solved! The HPPL function LOG computes the common log while in C (and BASIC) LOG returns the natural log. In the HPPL program change all the LOG to LN and the program works as expected.

well, Tom!
please, could you put a correct version in the Prime Software section of the Forum under something like "Bessel 2nd kind..."? Or, if you want could do it for you...

Salvo

I did it. It's called "Bessel Function Second Kind" (I don't know where I come up with these catchy names!)
Tom L
Cui bono?
Find all posts by this user
Quote this message in a reply
11-01-2017, 12:36 AM
Post: #7
salvomic Offline
Senior Member
 		Posts: 1,396
Joined: Jan 2015
RE: Bessel functions
(11-01-2017 12:22 AM)toml_12953 Wrote:  I did it. It's called "Bessel Function Second Kind" (I don't know where I come up with these catchy names!)

thank you, it works.
However I suggest to make the program "more friendly", not using Tbessly() to get the correct value for n=2, but with a function like bessel2(n,x) where the user could choose the number and the value to calc for, like the program for bessel1...
∫aL√0mic (IT9CLU) :: HP Prime 50g 41CX 71b 42s 39s 35s 12C 15C - DM42, DM41X - WP34s Prime Soft. Lib
Visit this user's website Find all posts by this user
Quote this message in a reply
11-01-2017, 04:35 AM (This post was last modified: 11-01-2017 04:41 AM by toml_12953.)
Post: #8
toml_12953 Offline
Senior Member
 		Posts: 2,200
Joined: Dec 2013
RE: Bessel functions
(11-01-2017 12:36 AM)salvomic Wrote:  
(11-01-2017 12:22 AM)toml_12953 Wrote:  I did it. It's called "Bessel Function Second Kind" (I don't know where I come up with these catchy names!)

thank you, it works.
However I suggest to make the program "more friendly", not using Tbessly() to get the correct value for n=2, but with a function like bessel2(n,x) where the user could choose the number and the value to calc for, like the program for bessel1...

It already has that function. It's called BESSY(N,X). The wrapper (Tbessy) calls it. Just insert the word EXPORT at the beginning of the BESSY routine and you can call BESSY from any other program or the home screen. I made the change to the program in the library. Thanks for the suggestion!
Tom L
Cui bono?
Find all posts by this user
Quote this message in a reply
11-01-2017, 09:11 AM
Post: #9
salvomic Offline
Senior Member
 		Posts: 1,396
Joined: Jan 2015
RE: Bessel functions
Thank you Tom!

I've seen it in the Library section, I tested it with some values in the literature and the function is very precise.
So, with this one and that of Eddie's we have both Bessel 1st and 2nd kind in the Prime!

Salvo Micciché
∫aL√0mic (IT9CLU) :: HP Prime 50g 41CX 71b 42s 39s 35s 12C 15C - DM42, DM41X - WP34s Prime Soft. Lib
Visit this user's website Find all posts by this user
Quote this message in a reply
Post Reply 




User(s) browsing this thread: 2 Guest(s)

Contact Us | The Museum of HP Calculators | Return to Top | Return to Content | Lite (Archive) Mode | RSS Syndication