Anger & Weber Functions for the HP-41
Overview
1°) Ascending Series
2°) Asymptotic Expansions
-Anger's function is defined by Jn(x)
= (1/PI) §0pi cos ( n.t - x.sin t ) dt
§ = integral
-Weber's function is defined by En(x)
= (1/PI) §0pi sin ( n.t - x.sin t ) dt
1°) Ascending Series
Formulae:
Jn(x) = + (x/2) sin ( 90°n
) 1F2( 1 ; (3-n)/2 , (3+n)/2 ; -x2/4
) / Gam((3-n)/2) / Gam ((3+n)/2)
+ cos ( 90°n ) 1F2( 1 ; (2-n)/2 , (2+n)/2 ; -x2/4
) / Gam((2-n)/2) / Gam ((2+n)/2)
En(x) = - (x/2) cos ( 90°n
) 1F2( 1 ; (3-n)/2 , (3+n)/2 ; -x2/4
) / Gam((3-n)/2) / Gam ((3+n)/2)
+ sin ( 90°n ) 1F2( 1 ; (2-n)/2 , (2+n)/2 ; -x2/4
) / Gam((2-n)/2) / Gam ((2+n)/2)
Data Registers: R00 to R07: temp
Flags: /
Subroutines: "1F2" ( cf "Hypergeometric
Functions" ) , "GAM" ( cf "Gamma Function" )
01 LBL "ANWEB"
02 2 03 ST/ Z 04 / 05 STO 04 06 X<>Y 07 STO 03 08 STO 05 09 CHS 10 1 11 STO 01 12 ST+ 03 |
13 +
14 STO 02 15 RCL 04 16 X^2 17 CHS 18 XEQ "1F2" 19 STO 06 20 RCL 02 21 XEQ "GAM" 22 ST/ 06 23 RCL 03 24 XEQ "GAM" |
25 ST/ 06
26 2 27 1/X 28 ST+ 02 29 ST+ 03 30 RCL 02 31 XEQ "GAM" 32 STO 07 33 RCL 03 34 XEQ "GAM" 35 ST* 07 36 RCL 00 |
37 XEQ "1F2"
38 RCL 07 39 / 40 RCL 04 41 * 42 1 43 CHS 44 ACOS 45 RCL 05 46 * 47 STO 05 48 X<>Y |
49 P-R
50 RCL 05 51 RCL 06 52 P-R 53 X<> Z 54 - 55 X<> Z 56 + 57 END |
( 100 bytes / SIZE 008 )
STACK | INPUTS | OUTPUTS |
Y | n | En(x) |
X | x | Jn(x) |
Example:
2 SQRT
PI XEQ "ANWEB" >>>> Jsqrt(2)(PI)
= 0.366086558
---Execution time = 26s---
X<>Y Esqrt(2)(PI) = -0.315594385
Notes:
-If n is an integer, Jn(x) = Jn(x)
= Bessel function of the first kind.
-However, this routine doesn't work if n is an integer, except if n
= -1 , 0 , +1
but it will work if we use the regularized hypergeometric
function 1F~2
-The variant hereunder calls "HGF+"
Data Registers: R00 to R06: temp
Flags: /
Subroutine: "HGF+" ( cf "Hypergeometric
Functions" )
-If you want to use the M-Code routine HGF+ , replace line 30 by
HGF+ and line 20 by STO 00 HGF+
01 LBL "ANWEB"
02 2 03 ST/ Z 04 / 05 STO 04 06 X<>Y 07 STO 03 08 STO 05 09 CHS 10 STO 02 |
11 1
12 STO 01 13 ST+ 02 14 ST+ 03 15 LASTX 16 CHS 17 RCL 04 18 X^2 19 CHS 20 XEQ "HGF+" |
21 STO 06
22 2 23 1/X 24 ST+ 02 25 ST+ 03 26 1 27 LASTX 28 CHS 29 RCL 00 30 XEQ "HGF+" |
31 RCL 04
32 * 33 1 34 CHS 35 ACOS 36 RCL 05 37 * 38 STO 05 39 X<>Y 40 P-R |
41 RCL 05
42 RCL 06 43 P-R 44 X<> Z 45 - 46 X<> Z 47 + 48 END |
( 75 bytes / SIZE 007 )
STACK | INPUTS | OUTPUTS |
Y | n | En(x) |
X | x | Jn(x) |
Examples:
2 SQRT
PI XEQ "ANWEB" >>>> Jsqrt(2)(PI)
= 0.366086558
---Execution time = 19s---
X<>Y Esqrt(2)(PI) = -0.315594386
5 ENTER^
PI R/S
>>>> J5(PI) = 0.052141184
---Execution time = 19s---
X<>Y E5(PI) = 0.207000292
2°) Asymptotic Expansions ( large x-values )
Formulae:
Jn(x) ~ Jn(x) + [ ( sin n.PI ) / ( x.PI ) ] ( F - n G / x )
En(x) ~ -Yn(x) - [ ( 1 + cos (n.PI) ) F - n ( 1 - cos n.PI ) G / x ] / ( x.PI )
where F = 1 + SUM k=1 to infinity
(n2-12) ....... (n2-(2k-1)2
/ x2k
and G = 1 + SUM k=1
to infinity (n2-22) ....... (n2-(2k)2
/ x2k
Data Registers: R00 = x , R01 = n
R02 = Jn(x) , R03 = Yn(x)
R04 to R06: temp
Flags: /
Subroutine: "AEJY" Asymptotic
series for Bessel Functions ( cf "Bessel Functions for the HP-41" §6°)
)
01 LBL "AEANW"
02 XEQ "AEJY" 03 STO 02 04 X<>Y 05 STO 03 06 RCL 01 07 X^2 08 1 09 STO 05 10 - 11 RCL 00 12 X^2 13 / 14 STO 06 15 RCL 05 |
16 +
17 XEQ 01 18 STO 04 19 CLX 20 STO 05 21 SIGN 22 STO 06 23 XEQ 01 24 RCL 01 25 * 26 RCL 00 27 / 28 RCL 04 29 X<>Y 30 ST+ 04 |
31 -
32 GTO 02 33 LBL 01 34 RCL 01 35 X^2 36 RCL 05 37 2 38 + 39 STO 05 40 X^2 41 - 42 RCL 06 43 * 44 RCL 00 45 X^2 |
46 /
47 STO 06 48 + 49 X#Y? 50 GTO 01 51 RTN 52 LBL 02 53 RCL 01 54 PI 55 R-D 56 * 57 X<>Y 58 RCL 00 59 PI 60 * |
61 ST/ 04
62 / 63 P-R 64 RCL 04 65 + 66 CHS 67 RCL 03 68 - 69 X<>Y 70 RCL 02 71 + 72 END |
( 95 bytes / SIZE 007 )
STACK | INPUTS | OUTPUTS |
Y | n | En(x) |
X | x | Jn(x) |
Examples:
PI ENTER^
24.4 XEQ "AEANW" >>>> JPI(24.4)
= 0.157155848
---Execution time = 26s---
X<>Y EPI(24.4) = -0.008969347
Notes:
-With n = PI , the asymptotic series diverges too early if x = 24.3
and we fall into an infinite loop.
-We can also express F & G in terms of hypergeometric functions
and we have:
F = 3F0( (1+n)/2 , (1-n)/2 , 1 ; -4/x2 ) and G = 3F0( (2+n)/2 , (2-n)/2 , 1 ; -4/x2 )
-So, the M-Code function HGF+ may be used to get faster results:
Data Registers: R00 = x , R10 = n
R04 = Jn(x) , R05 = Yn(x)
Flags: /
Subroutine: "AEJY" ( cf "Bessel
Functions for the HP-41" §6°) ) - the 2nd variant that uses
HGF+
01 LBL "AEANW"
02 X<>Y 03 STO 10 04 X<>Y 05 XEQ "AEJY" 06 STO 04 07 X<>Y 08 STO 05 09 1 10 STO 03 11 RCL 10 12 + 13 2 |
14 /
15 STO 01 16 RCL 03 17 RCL 10 18 - 19 2 20 / 21 STO 02 22 3 23 0 24 LASTX 25 RCL 00 26 / |
27 X^2
28 CHS 29 HGF+ 30 STO 06 31 .5 32 ST+ 01 33 ST+ 02 34 3 35 0 36 LASTX 37 HGF+ 38 RCL 10 39 * |
40 RCL 00
41 / 42 RCL 06 43 X<>Y 44 ST+ 06 45 - 46 RCL 10 47 PI 48 R-D 49 * 50 X<>Y 51 RCL 00 52 PI |
53 *
54 ST/ 06 55 / 56 P-R 57 RCL 06 58 + 59 CHS 60 RCL 05 61 - 62 X<>Y 63 RCL 04 64 + 65 END |
( 89 bytes / SIZE 011 )
STACK | INPUTS | OUTPUTS |
Y | n | En(x) |
X | x | Jn(x) |
Examples:
PI ENTER^
24.4 XEQ "AEANW" >>>> JPI(24.4)
= 0.157155848
---Execution time = 16s---
X<>Y EPI(24.4) = -0.008969347
References:
[1] Abramowitz and Stegun - "Handbook of Mathematical Functions"
- Dover Publications - ISBN 0-486-61272-4
[2] http://mathworld.wolfram.com/AngerFunctions.html
[3] http://mathworld.wolfram.com/WeberFunctions.html
[4] http://dlmf.nist.gov/11.11