Weierstrass Elliptic Functions for the HP-41
Overview
1°) A Laurent Series ( real arguments )
2°) Duplication formula ( real arguments
)
3°) A Laurent Series ( complex arguments
)
4°) Duplication formula ( complex arguments
)
5°) Weierstrass & Jacobian Elliptic
Functions - Half Periods
1°) A Laurent Series ( real arguments
)
-The Weierstrass Elliptic Function P(x;g2;g3)
satisfies the differential equation: (P')2
= 4.P3 - g2.P - g3
-This program calculates P(x;g2;g3)
by a Laurent series:
P(x;g2;g3) = x-2 + c2.x2 + c3.x4 + ...... + cn.x2n-2 + ....
where c2 = g2/20 ; c3 = g3/28 and cn = 3 ( c2. cn-2 + c3. cn-3 + ....... + cn-2. c2 ) / (( 2n+1 )( n-3 )) ( n > 3 )
-The successive cn are computed and stored into registers R05 R06 ......
Data Registers: R00 to R04 are used for temporary
data storage and R05 = c2 ; R06 = c3 ; R07
= c4 ; ......
Flags: /
Subroutines: /
| 01 LBL "WEF"
02 X^2 03 STO 02 04 CLX 05 20 06 / 07 STO 05 08 X<>Y 09 28 10 / 11 STO 06 12 RCL 02 13 * 14 + |
15 STO 01
16 5 17 STO 00 18 LBL 01 19 2 20 STO 04 21 LBL 02 22 RCL 00 23 .1 24 % 25 5 26 + 27 STO 03 28 CLX |
29 LBL 03
30 RCL IND Y 31 RCL IND 03 32 * 33 + 34 DSE Y 35 ISG 03 36 GTO 03 37 3 38 * 39 RCL 00 40 ENTER^ 41 ST+ Y 42 DSE Y |
43 4
44 - 45 * 46 / 47 ISG 03 48 CLX 49 STO IND 03 50 RCL 02 51 RCL 00 52 3 53 - 54 Y^X 55 * 56 RCL 01 |
57 +
58 STO 01 59 ISG 00 60 CLX 61 LASTX 62 X#Y? 63 GTO 01 64 DSE 04 65 GTO 02 66 RCL 02 67 ST* Y 68 1/X 69 + 70 END |
( 95 bytes / SIZE ? )
| STACK | INPUTS | OUTPUTS |
| Z | g3 | / |
| Y | g2 | / |
| X | x | P(x;g2;g3) |
Example: Calculate
P(x;g2;g3) for x = 0.6
g2 = 0.9 g3 = 1.4
1.4 ENTER^
0.9 ENTER^
0.6 XEQ "WEF" >>>> 2.800500784
( in 32 seconds )
-The SIZE must be ( at least ) 016 in this example.
-If you get the message "NON EXISTENT" ( line 49 ) , increase the SIZE
and R/S
-The Laurent series may diverge if x is too great. In this case, one
may use the duplication formula below.
2°) Duplication formula ( real arguments
)
-We have P(2x) = -2 P(x) + ( 6 P2(x)
- g2/2 )2 / ( 4 ( 4 P3(x) - g2
P(x) - g3 ) )
-This formula is used recursively ( n times ) to obtain P(2nx)
if we know P(x) ( n = 1 ; 2 ; 3 ;
..... )
Data Registers: • R00 = n ( n > 0 ) ( Register R00 is to be initialized before executing "DUPW" )
R01: temp R02 = g2 R03 = g3
Flags: /
Subroutines: /
| 01 LBL "DUPW"
02 RDN 03 STO 02 04 X<>Y 05 STO 03 06 R^ 07 LBL 01 |
08 STO 01
09 X^2 10 6 11 * 12 RCL 02 13 2 14 / |
15 -
16 X^2 17 RCL 01 18 ST+ X 19 X^2 20 RCL 02 21 - |
22 RCL 01
23 * 24 RCL 03 25 - 26 4 27 * 28 / |
29 RCL 01
30 ST+ X 31 - 32 DSE 00 33 GTO 01 34 END |
( 47 bytes / SIZE 004 )
| STACK | INPUTS | OUTPUTS |
| Z | g3 | / |
| Y | g2 | / |
| X | P(x;g2;g3) | P(2nx;g2;g3) |
Example: Calculate P(x;g2;g3)
for x = 4.8 g2 = 0.9 g3
= 1.4
-The Laurent series diverges for these arguments but 4.8 = 23* 0.6 and we have already found that P(0.6;0.9;1.4) = 2.800500784 , thus:
3 STO 00
1.4 ENTER^
0.9 ENTER^
2.800500784 XEQ "DUPW" >>>> P(4.8;0.9;1.4)
= 1.954704160 ( in 3seconds )
3°) A Laurent Series ( complex arguments
)
-The same Laurent series is now applied to z = x + i.y
( but g2 and g3 real )
Data Registers: R00 to R06 are used for temporary
data storage and R07 = c2 ; R08 = c3 ; R09
= c4 ; ......
Flags: /
Subroutines: /
| 01 LBL "WEFZ"
02 R-P 03 X^2 04 STO 03 05 RDN 06 ST+ X 07 STO 04 08 CLX 09 20 10 / 11 STO 07 12 CLX 13 28 14 / 15 STO 08 16 RCL 03 17 * 18 RCL 04 19 X<>Y 20 P-R 21 RCL 07 22 + |
23 STO 01
24 X<>Y 25 STO 02 26 7 27 STO 00 28 LBL 01 29 2 30 STO 06 31 LBL 02 32 RCL 00 33 .1 34 % 35 7 36 + 37 STO 05 38 CLX 39 LBL 03 40 RCL IND Y 41 RCL IND 05 42 * 43 + 44 DSE Y |
45 ISG 05
46 GTO 03 47 3 48 * 49 RCL 00 50 ENTER^ 51 ST+ Y 52 DSE X 53 5 54 ST- Z 55 - 56 * 57 / 58 ISG 05 59 CLX 60 STO IND 05 61 RCL 04 62 RCL 03 63 RCL 00 64 5 65 - 66 ST* Z |
67 Y^X
68 RCL IND 05 69 * 70 P-R 71 RCL 01 72 + 73 STO 01 74 LASTX 75 - 76 X<>Y 77 RCL 02 78 + 79 STO 02 80 LASTX 81 - 82 R-P 83 ISG 00 84 CLX 85 X#0? 86 GTO 01 87 RCL 02 88 RCL 01 |
89 R-P
90 RCL 03 91 * 92 X<>Y 93 RCL 04 94 + 95 X<>Y 96 P-R 97 RCL 04 98 CHS 99 RCL 03 100 1/X 101 P-R 102 RDN 103 ST+ Z 104 X<> T 105 + 106 END |
( 135 bytes / SIZE? )
| STACK | INPUTS | OUTPUTS |
| T | g3 | / |
| Z | g2 | / |
| Y | y | B |
| X | x | A |
with P(x+iy;g2;g3) = A + i.B
Example: Calculate P(z;g2;g3) for z = 0.6 + 0.4 i g2 = 0.9 g3 = 1.4
1.4 ENTER^
0.9 ENTER^
0.4 ENTER^
0.6 XEQ "WEFZ" >>>> 0.7390436447
X<>Y -1.744031391 ( in 48seconds
, SIZE 018 ( at least ) )
Whence P(0.6+0.4i ;0.9;1.4)
= 0.7390436447 -1.744031391 i
4°) Duplication formula ( complex arguments
)
-The same duplication formula holds for complex arguments.
Data Registers: • R00 = n ( n > 0 ) ( Register R00 is to be initialized before executing "DUPWZ" )
R01: temp R02 = g2 R03 = g3
R04 to R06: temp
Flags: /
Subroutines: /
| 01 LBL "DUPWZ"
02 R^ 03 STO 03 04 R^ 05 STO 02 06 R^ 07 R^ 08 LBL 01 09 STO 01 10 X<>Y 11 STO 04 12 X<>Y 13 R-P 14 3 |
15 ST* Z
16 Y^X 17 4 18 * 19 P-R 20 RCL 04 21 RCL 02 22 * 23 ST- Z 24 CLX 25 RCL 01 26 LASTX 27 * 28 - |
29 RCL 03
30 - 31 R-P 32 STO 05 33 X<>Y 34 STO 06 35 RCL 04 36 RCL 01 37 R-P 38 2 39 ST* Z 40 Y^X 41 6 42 * |
43 P-R
44 RCL 02 45 2 46 / 47 - 48 R-P 49 2 50 ST/ Y 51 ST* Z 52 Y^X 53 RCL 06 54 ST- Z 55 CLX 56 RCL 05 |
57 /
58 P-R 59 RCL 04 60 ST+ X 61 ST- Z 62 CLX 63 RCL 01 64 ST+ X 65 - 66 DSE 00 67 GTO 01 68 END |
( 89 bytes / SIZE 007 )
| STACK | INPUTS | OUTPUTS |
| T | g3 | / |
| Z | g2 | / |
| Y | B | B' |
| X | A | A' |
with P(z;g2;g3) = A + i.B and P(2nz;g2;g3) = A' + i.B'
Example: Deduce P(4.8+3.2 i ;0.9;1.4) from P(0.6+0.4i ;0.9;1.4) = 0.7390436447 -1.744031391 i
3 STO 00
1.4 ENTER^
0.9 ENTER^
0.4 ENTER^
0.6 XEQ "DUPWZ" >>>> 0.152165025
X<>Y -1.254979593 ( in 22 seconds
)
whence P(4.8+3.2 i ;0.9;1.4)
= 0.152165025 - 1.254979593 i
5°) Weierstrass & Jacobian Elliptic Functions
- Half Periods
-The Weierstrass Elliptic Functions may also be computed by mean of
the Jacobian Elliptic Functions.
-This program also gives the first derivative and the half-periods
of P(x;g2;g3)
-We have P(x+2k*Omega+2ik'*Omega') = P(x)
for all x ( except the poles ) if k
and k' are integers & Omega and i.Omega' are the half-periods.
Formulae: Let r1
; r2 ; r3 the 3 roots of
the polynomial 4x3- g2.x -g3
and m the parameter of the complete elliptic integrals
of the 1st kind: K(m)
IF THE 3 ROOTS ARE REAL ( R1 > R2 > R3 ) IF ONLY ONE ROOT IS REAL ( SAY R1 )
m = ( r2 - r3 )/( r1 - r3
)
m = 1/2 - 3 r1/(4H)
Omega = K(m)/( r1 - r3 )1/2
Omega2 = K(m)/H1/2
Omega' = K(1-m)/( r1 - r3 )1/2
Omega'2 = K(1-m)/H1/2
P(x;g2;g3)
= r3 + ( r1 - r3 )/sn2( y |
m )
P(x;g2;g3) = r1 + H ( 1 + cn(y'|m)
) / ( 1 - cn(y'|m) )
P'(x;g2;g3)
= -2 ( r1 - r3 )3/2 cn (y|m) dn(y|m) /
sn3(y|m)
P'(x;g2;g3) = -4 H3/2 sn(y'|m)
dn(y'|m) / ( 1 - cn(y'|m) )2
with y = ( r1 - r3 )1/2.x
with y' = 2.x.H1/2 and
H = ( 2.r12 + r2.r3 )1/2
Omega2 = Omega + i. Omega'
Omega'2 = Omega' + i. Omega
Data Registers: R00 to R07 are used for temporary data storage by the Jacobian Elliptic Function and Complete Elliptic Integral programs
and when the program stops:
R08 = x
R09 = Omega if flag F01 is clear , R09 = Omega2
if flag F01 is set , R11: temp
R10 = Omega2 if flag F01 is clear , R10 = Omega'2 if
flag F01 is set , R12: temp
Flags: F01 - F02 - F05
Subroutines:
"CEI" ( Complete Ellptic Integrals ) see "Jacobian
Elliptic Functions for the HP-41"
"JEF" ( Jacobian Elliptic Functions ) see
"Jacobian Elliptic Functions for the HP-41"
"P3" ( Cubic Equations ) see "Polynomials
for the HP-41" or another program which returns the larger real solution
of a cubic in X-register.
Note:
-If you have used registers R08 & R09 instead of synthetic registers
M & N in the "JEF" program,
replace R08 and R09 by R13 and R14 in the following listing.
| 01 LBL "WEF2"
02 STO 08 03 RDN 04 CHS 05 0 06 X<> Z 07 CHS 08 4 09 RDN 10 XEQ "P3" 11 FS? 01 12 GTO 02 13 RCL Z 14 STO 11 15 ST- Z 16 - 17 STO 12 18 XEQ 01 19 ST/ Y 20 X^2 21 STO 01 22 / 23 * |
24 RCL 01
25 1/X 26 GTO 03 27 LBL 01 28 CF 05 29 / 30 STO 10 31 1 32 X<>Y 33 X=Y? 34 SF 05 35 FC? 05 36 XEQ "CEI" 37 X<>Y 38 FS?C 05 39 E90 40 STO 09 41 1 42 STO Y 43 RCL 10 44 - 45 X=Y? 46 SF 05 |
47 FC? 05
48 XEQ "CEI" 49 X<>Y 50 FS?C 05 51 E90 52 X<> 10 53 RCL 12 54 SQRT 55 ST/ 09 56 ST/ 10 57 RCL 08 58 FS? 01 59 ST+ X 60 * 61 XEQ "JEF" 62 RTN 63 LBL 02 64 STO 11 65 X^2 66 ST+ X 67 RDN 68 R-P 69 X^2 |
70 R^
71 + 72 SQRT 73 STO 12 74 RCL 11 75 3 76 * 77 X<>Y 78 / 79 2 80 - 81 CHS 82 4 83 XEQ 01 84 SF 01 85 X<>Y 86 STO 01 87 1 88 - 89 STO 02 90 X^2 91 / 92 * |
93 ST+ X
94 1 95 RCL 01 96 + 97 RCL 02 98 CHS 99 / 100 LBL 03 101 RCL 12 102 SQRT 103 ST+ X 104 CHS 105 ST* Z 106 X<> L 107 ST* Z 108 * 109 RCL 11 110 + 111 END |
( 170 bytes / SIZE 013 )
| STACK | INPUTS | OUTPUTS |
| Z | g3 | / |
| Y | g2 | P'(x;g2;g3) |
| X | x | P(x;g2;g3) |
Example1: Calculate P(x;g2;g3)
& P'(x;g2;g3) for x
= 2 g2 = 4 g3 = 1
1 ENTER^
4 ENTER^
2 XEQ "WEF2" >>>> P(2;4;1)
= 4.950267724 X<>Y
P'(2;4;1)
= 21.55057197 ( in 28 seconds )
-We have R09 = 1.225694692 & R10 = 1.496729323 ( Omega & Omega' because F01 is clear )
Therefore the primitive half-periods
are: 1.225694692 & 1.496729323
i
Example2: Calculate P(x;g2;g3) & P'(x;g2;g3) for x = 1 g2 = 2 g3 = 3
3 ENTER^
2 ENTER^
1 XEQ "WEF2" >>>> P(1;2;3)
= 1.214433709 X<>Y
P'(1;2;3) = -1.317406193 ( in 27
seconds )
-We have R09 = 1.197220889 & R10 = 2.350281226 ( Omega2 & Omega'2 because F01 is set )
Whence:
Omega = 0.598610445 - 1.175140613 i &
Omega' = 0.598610445 + 1.175140613 i
Notes:
-I've called Omega' & Omega'2 what Abramowitz and Stegun call
Omega'/i & (Omega'2)/i
-This program doesn't work if g2 = g3 =
0 but in this case, P(x) = 1/x2
-Lines 28-32 -33-34-35-38-39-42-45-46-47-50-51 are useful only
if 2 roots of the polynomial 4x3- g2.x -g3
are equal:
in this case, one of the half-periods is infinite and a number
near E90 is stored in register R09 or R10.
-Without these lines, there would be a DATA ERROR at line 38 of the "CEI" program ( K(1) = infinite )
For instance P(1;48;-64) = 2.181597704 ;
P'(1;48;-64)
= -0.903006185 ; R09 = 4.082 1089 = infinite =
Omega ; R10 = 0.641274915 = Omega' ( CF01 )
P(1;12;8)
= 2.079381536 ; P'(1;12;8) = 1.735214810 ; R09 = 0.906899682
= Omega ; R10 = 5.7735 1089 = infinite = Omega'
( CF01 )
-Numerous other relations with Jacobian Elliptic Functions and Theta
Functions may be found in the reference below
Reference:
[1] Abramowitz and Stegun , "Handbook of Mathematical Functions" -
Dover Publications - ISBN 0-486-61272-4