Bessel functions
10-30-2017, 10:21 PM
Post: #1 salvomic Senior Member Posts: 1,394 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
10-31-2017, 01:12 AM (This post was last modified: 11-01-2017 11:48 AM by toml_12953.)
Post: #2
 toml_12953 Senior Member Posts: 1,949 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?
10-31-2017, 08:37 AM
Post: #3 salvomic Senior Member Posts: 1,394 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
10-31-2017, 09:14 PM
Post: #4
 toml_12953 Senior Member Posts: 1,949 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?
10-31-2017, 11:56 PM
Post: #5 salvomic Senior Member Posts: 1,394 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
11-01-2017, 12:22 AM
Post: #6
 toml_12953 Senior Member Posts: 1,949 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?
11-01-2017, 12:36 AM
Post: #7 salvomic Senior Member Posts: 1,394 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
11-01-2017, 04:35 AM (This post was last modified: 11-01-2017 04:41 AM by toml_12953.)
Post: #8
 toml_12953 Senior Member Posts: 1,949 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?
11-01-2017, 09:11 AM
Post: #9 salvomic Senior Member Posts: 1,394 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
 « Next Oldest | Next Newest »

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