The Museum of HP Calculators

# Orthogonal Polynomials for the HP-41

This program is Copyright © 2006 by Jean-Marc Baillard and is used here by permission.

This program is supplied without representation or warranty of any kind. Jean-Marc Baillard and The Museum of HP Calculators therefore assume no responsibility and shall have no liability, consequential or otherwise, of any kind arising from the use of this program material or any part thereof.

Overview

1°) Legendre Polynomials
2°) Laguerre Polynomials
3°) Generalized Laguerre Polynomials
4°) Hermite Polynomials
5°) Chebyshev Polynomials & the "Connaissanse des Temps"

a) Chebyshev Polynomials of the first and second kind
b) The "Connaissance des Temps"

6°) Ultraspherical Polynomials
7°) Jacobi's Polynomials

1°) Legendre Polynomials

Formulae:      n.Pn(x) = (2n-1).x.Pn-1(x) - (n-1).Pn-2(x)  ;  P0(x) = 1  ;  P1(x) = x

Data Registers:  /
Flags: /
Subroutines: /

-Synthetic register M may be replaced by any unused data register.

01  LBL "LEG"
02  X<>Y
03   E3
04  /
05  1
06  +
07  STO M
08  LASTX
09  LBL 01
10  STO N
11  RCL Z
12  *
13  STO O
14  -
15  RCL M
16  INT
17  1/X
18  1
19  -
20  *
21  RCL O
22  +
23  RCL N
24  X<>Y
25  ISG M
26  GTO 01
27  CLA
28  END

( 46 bytes / SIZE 000 )

 STACK INPUTS OUTPUTS T / x Z / x Y n Pn-1(x) X x Pn (x)

( 0 < n < 999 )

Example:    Calculate   P7(4.9)

7    ENTER^
4.9   XEQ "LEG"  gives   P7(4.9) =  1698444.019    ( in 4 seconds )   and   P6(4.9) = 188641.3852  is in Y-register.

2°) Laguerre Polynomials

Formulae:      Ln(x) = (2n-1-x).Ln-1(x) - (n-1)2.Ln-2(x)  ;  L0(x) = 1  ;  L1(x) = 1 - x

Data Registers:  /
Flags: /
Subroutines: /

-Synthetic register M may be replaced by any unused data register.

01  LBL "LAG"
02  X<>Y
03  .1
04  %
05  1
06  +
07  STO M
08  SIGN
09  LBL 01
10  X<>Y
11  RCL M
12  INT
13  DSE X
14  ""                TEXT0   or another NOP instruction like  LBL 02  STO X   ABS  ...  etc  ...
15  X^2
16  *
17  RCL M
18  INT
19  ST+ X
20  DSE X
21  R^
22  ST- Y
23  X<> T
24  ST* Y
25  X<> Z
26  -
27  ISG M
28  GTO 01
29  CLA
30  END

( 51 bytes / SIZE 000 )

 STACK INPUTS OUTPUTS T / x Z / x Y n Ln-1(x) X x Ln (x)

( 0 < n < 1000 )

Example:    Calculate   L7 (3.14)

7       ENTER^
3.14  XEQ "LAG"  >>>>  L7 (3.14) =   -4932.43995       ( in 5 seconds )
X<>Y       >>>>   L6 (3.14) =    -189.9784735

Note:   -Some authors divide Ln (x) by  n!     In this case simply add   RCL M  1  -  INT   FACT   /    after line 28
( but Y Z and T-registers now contain the former Ln-1 (x) )
-With this definition,  L7 (3.14) =  -0.978658720

3°) Generalized Laguerre Polynomials

Formulae:      L0(a) (x) = 1  ;   L1(a) (x) = a+1-x   ;    n.Ln(a) (x) = (2.n+a-1-x).Ln-1(a) (x) - (n+a-1).Ln-2(a) (x)

Data Registers: /
Flags: /
Subroutines: /

01  LBL "LANX"
02  STO M
03  CLX
04   E3
05  /
06  1
07  ST- Z
08  +
09  STO N
10  X<>Y
11  STO O
12  CLST
13  SIGN
14  LBL 01
15  X<>Y
16  CHS
17  RCL N
18  INT
19  RCL O
20  +
21  *
22  RCL N
23  INT
24  ST+ X
25  RCL O
26  +
27  RCL M
28  -
29  R^
30  *
31  +
32  RCL N
33  INT
34  /
35  ISG N
36  GTO 01
37  CLA
38  END

( 61 bytes / SIZE 000 )

 STACK INPUTS OUTPUTS Z a Ln-1(a)(x) Y n Ln-1(a)(x) X x Ln(a)(x)

( 0 < n < 1000 ,  n integer )

Example:

1.4  ENTER^
7    ENTER^
PI   XEQ "LANX"  >>>>   L7(1.4)(Pi) = 1.688893513       ( 5 seconds )
X<>Y   L6(1.4)(Pi) = 2.271353727

Note:   Ln(0)(x) = (1/n!).Ln(x)   where   Ln(x) = Laguerre's Polynomial

4°) Hermite Polynomials

Formulae:      Hn(x) = 2x.Hn-1(x) - 2(n-1).Hn-2(x)  ;  H0(x) = 1  ;  H1(x) = 2x

Data Registers:  /
Flags: /
Subroutines: /

-Synthetic register M may be replaced by any unused data register.

01  LBL "HMT"
02  X<>Y
03  .1
04  %
05  1
06  +
07  STO M
08  SIGN
09  LBL 01
10  X<>Y
11  RCL M
12  INT
13  DSE X
14  ""                TEXT0   or another NOP instruction like  LBL 02  STO X   ABS  ...  etc  ...
15  *
16  CHS
17  X<>Y
18  ST* Z
19  X<> Z
20  +
21  ST+ X
22  ISG M
23  GTO 01
24  CLA
25  END

( 42 bytes / SIZE 000 )

 STACK INPUTS OUTPUTS T / x Z / x Y n Hn-1(x) X x Hn (x)

( 0 < n < 1000 )

Example:    Calculate    H7 (3.14)

7       ENTER^
3.14  XEQ "HMT"  >>>>  H7 (3.14) =   73726.24332       ( in 4 seconds )
X<>Y       >>>>   H6 (3.14) =   21659.28040

5°) Chebyshev Polynomials & the "Connaissance des Temps"

a) Chebyshev Polynomials of the first and second kind

Formulae:    Tn(x) = 2x.Tn-1(x) - Tn-2(x)  ;  T0(x) = 1  ;  T1(x) = x   ( first kind )  &   Un(x) = 2x.Un-1(x) - Un-2(x)  ;  U0(x) = 1  ;  U1(x) = 2x    ( second kind )

Data Registers:  /
Flag: F02               if F02 is clear, "CHB" calculates the Chebyshev polynomials of the first kind:   Tn(x)
if F02 is set,    "CHB" calculates the Chebyshev polynomials of the second kind:   Un(x)
Subroutines: /

-Synthetic register M may be replaced by any unused data register.

01  LBL "CHB"
02  STO M
03  ST+ M
04  FS? 02
05  CLX
06  X<>Y
07   E3
08  /
09  1
10  +
11  STO Z
12  SIGN
13  LBL 01
14  RCL M
15  X<>Y
16  ST* Y
17  X<> Z
18  -
19  ISG Z
20  GTO 01
21  CLA
22  END

( 40 bytes / SIZE 000 )

 STACK INPUTS OUTPUTS Y n Tn-1(x)orUn-1(x) X x Tn (x)orUn(x)

( 0 < n < 1000 )

Example:    Compute    T7 (0.314)  and  U7(0.314)

CF 02
7       ENTER^
0.314  XEQ "CHB"  >>>>   T7 (0.314) =  -0.786900700   and     T6 (0.314) =  0.338782777  in Y-register   ( in 2 seconds )

SF 02
7       ENTER^
0.314      R/S          >>>>    U7 (0.314) =  -0.582815681   and     U6 (0.314) =  0.649952293  in Y-register

b) The "Connaissance des Temps"

-The "Connaissance des Temps" is a French ephemeride which contains coefficients of Chebyshev polynomials
for very accurate coordinates of the Sun , the Moon and the planets.
-If y(t) is a coordinate of a celestial body at a given instant t ( t belonging to the interval  [ t0 ; t0 + DT ] ) ,  y is given by the formula:

y = a0 + a1.T1(x) + ....... + an.Tn(x)             where  x = -1 + 2(t-t0)/DT      ( -1 <= x <= +1 )
and    a0 ; a1 ; ....... ; an  are published in the Connaissance des Temps

Data Registers:  only the (n+1) registers used to store the coefficients   a0 ; a1 ; ....... ; an
Flags: /
Subroutines: /

-Synthetic registers M N O may be replaced by any unused data registers.

01  LBL "CdT"
02  X<>Y
03  STO M
04  SIGN
05  ST+ M
06  STO N
07  RCL Y
08  STO O
09  RCL IND L
10  LBL 01
11  R^
12  ST+ X
13  RCL N
14  ST* Y
15  X<> O
16  -
17  STO N
18  RCL IND M
19  *
20  +
21  ISG M
22  GTO 01
23  CLA
24  END

( 46 bytes )

 STACK INPUTS OUTPUTS T / x Z / x Y bbb.eee x X x y(x)

where  bbb.eee  is the control number of the registers which contain the (n+1) coefficients   a0 ; a1 ; ....... ; an

Example:    Compute the radius vector of Saturn for 2000 March 12  at  0h  TT       ( TT = TAI+32.184s )

-We find the following coefficients in the "Connaissance des Temps 2000" ( page B37 )  valid for 2000/01/00  0h   to  2000/12/33  0h     (  DT = 368 days )

a0 = 9.14765315 ;  a1 =  -0.03544281  ;  a2 = 0.00109597  ;  a3 = 0.00002140  ;  a4 = 0.00000039  ;   a5 = -0.00000083    ( in AU )

-For instance, we store these 6 numbers in registers R10 to R15                         ( bbb.eee = 10.015  )
-In this example t-t0 = 72 days therefore  x = -1 +2*72/368 = -0.6086956522

10.015          ENTER^
-0.6086956522   XEQ "CdT"  >>>>   radius vector of Saturn = 9.16896253  AU   ( in 2 seconds )

6°) Ultraspherical Polynomials

Formulae:      C0(a) (x) = 1  ;   C1(a) (x) = 2.a.x   ;    (n+1).Cn+1(a) (x) = 2.(n+a).x.Cn(a) (x) - (n+2a-1).Cn-1(a) (x)      if  a # 0

Cn(0) (x) = (2/n).Tn(x)

Data Registers:  /
Flags:                F02   if  a = 0                                          |   none  if  a # 0
Subroutines:   CHB - Chebyshev Polynomials   if  a =0   |   -------------

01  LBL "USP"
02  1
03  R^
04  X#0?
05  GTO 00
06  CF 02
07  R^
08  R^
09  XEQ "CHB"
10  ST+ X
11  R^
12  INT
13  /
14  RTN
15  LBL 00
16  -
17  ST+ X
18  STO O
19  RDN
20  STO M
21  CLX
22   E3
23  /
24  1
25  +
26  STO N
27  CLST
28  SIGN
29  LBL 01
30  X<>Y
31  RCL O
32  RCL N
33  INT
34  -
35  ST* Y
36  X<> L
37  ST+ X
38  RCL O
39  -
40  RCL M
41  *
42  R^
43  *
44  +
45  RCL N
46  INT
47  /
48  ISG N
49  GTO 01
50  CLA
51  END

( 81 bytes / SIZE 000 )

 STACK INPUTS OUTPUTS Z a / Y n Cn-1(a)(x) X x Cn(a)(x)

( 0 < n < 1000 ,  n integer )

Cn-1(a)(x)  is in Y-register only if  a # 0

Example:    1.5  ENTER^
7    PI   1/X   XEQ "USP"  >>>>   C7(1.5)(1/Pi)  =  -0.989046386
X<>Y   C6(1.5)(1/Pi)  =    1.768780932

7°) Jacobi's Polynomials

Formulae:      P0(a;b) (x) = 1  ;   P1(a;b) (x) = (a-b)/2 + x.(a+b+2)/2

2n(n+a+b)(2n+a+b-2).Pn(a;b) (x) = [ (2n+a+b-1).(a2-b2) + x.(2n+a+b-2)(2n+a+b-1)(2n+a+b) ] Pn-1(a;b) (x)
- 2(n+a-1)(n+b-1)(2n+a+b) Pn-2(a;b) (x)

Data Registers:         R00 thru R05: temp    ( R01 = x )
Flags: /
Subroutines: /

01  LBL"JCP"
02  STO 01
03  CLX
04   E3
05  /
06  1
07  STO 04
08  +
09  STO 00
10  RDN
11  STO 03
12  X<>Y
13  STO 02
14  +
15  2
16  +
17  RCL 01
18  *
19  RCL 02
20  RCL 03
21  -
22  +
23  2
24  LBL 01
25  /
26  ISG 00
27  FS? 30
28  GTO 02
29  X<> 04
30  CHS
31  RCL 00
32  INT
33  RCL 02
34  +
35  STO 05
36  1
37  -
38  *
39  RCL 00
40  INT
41  RCL 03
42  +
43  ST+ 05
44  1
45  -
46  *
47  RCL 05
48  ENTER^
49  ST* Z
50  1
51  -
52  STO 05
53  ST* Y
54  LASTX
55  -
56  *
57  RCL 01
58  *
59  RCL 02
60  X^2
61  RCL 03
62  X^2
63  -
64  RCL 05
65  *
66  +
67  RCL 04
68  *
69  2
70  /
71  +
72  RCL 05
73  1
74  -
75  RCL 02
76  RCL 03
77  +
78  RCL 00
79  INT
80  ST+ Y
81  *
82  *
83  GTO 01
84  LBL 02
85  RCL 04
86  X<>Y
87  END

( 105 bytes / SIZE 006 )

 STACK INPUTS OUTPUTS T a / Z b / Y n Pn-1(a;b)(x) X x Pn(a;b)(x)

( 0 < n < 1000 ,  n integer )

Example:

1.4  ENTER^
1.7  ENTER^
7    ENTER^
PI   1/X  XEQ "JCP"  >>>>   P7(1.4;1.7)(1/Pi) = -0.322234421     ( in 13 seconds )
X<>Y  P6(1.4;1.7)(1/Pi) =   0.538220923

Reference:         Abramowitz and Stegun , "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4