Numerov

Numerov's Method for the HP-41

Overview

1°) Numerov's method

a) 1 differential equation      y" = f(x,y)
b) 2 differential equations    y" = f(x,y,z) ; z" = g(x,y,z)
c) n differential equations    y"₁ = f₁(x;y₁,.....,y_n) , ....... , y"_n = f_n(x;y₁,.....,y_n)     ( X-Functions Module required )

2°) A formula of order 7

-Numerov's method has been devised to solve second order differential equations where the first derivative y' does not appear explicitly.
-The formula is: y_n+1 = 2.y_n - y_n-1 + (h²/12).( f_n+1 + 10.f_n + f_n-1 ) where h = x_n - x_n-1 is the stepsize and f_k = f(x_k,y_k) ; y_n-j = y(x_n-j.h)

-It requires 2 starting values y_-1 and y₀ instead of y₀ and y'₀ .
-The method duplicates the Taylor series up to the term in h⁵

-Since it is an implicit formula ( y_n+1 appears on both sides ) , we use an iterative method at each step.

1°) Numerov's Method

a) 1 differential equation y" = f(x,y)

Data Registers: • R00 = f name ( Registers R00 thru R05 are to be initialized before executing "NUM1" )

                                      • R01 = x₀
                                      • R02 = y₀   • R03 = y_-1 = y(x₀ - h )    • R04 = h = stepsize    • R05 = N = number of steps      R06 to R10: temp
Flags: /
Subroutine:   A program which computes f(x,y) assuming x is in X-register and y is in Y-register upon entry

-Line 40 is only useful to display the successive approximations, but it's not necessary for the calculations.

01 LBL "NUM1"
02 RCL 05
03 STO 06
04 RCL 03
05 RCL 01
06 RCL 04
07 -
08 XEQ IND 00
09 STO 07
10 RCL 02
11 STO 10

12 RCL 01
13 XEQ IND 00
14 STO 08
15 LBL 01
16 RCL 10
17 RCL 01
18 RCL 04
19 +
20 XEQ IND 00
21 STO 09
22 RCL 08

23 10
24 *
25 +
26 RCL 07
27 +
28 RCL 04
29 X^2
30 *
31 12
32 /
33 RCL 02

34 ST+ X
35 RCL 03
36 -
37 +
38 ENTER^
39 X<> 10
40 VIEW X
41 X#Y?
42 GTO 01
43 X<> 02
44 STO 03

45 RCL 04
46 ST+ 01
47 RCL 09
48 X<> 08
49 STO 07
50 DSE 06
51 GTO 01
52 RCL 02
53 RCL 01
54 CLD
55 END

( 78 bytes / SIZE 011 )

STACK	INPUTS	OUTPUTS
Y	/	y(x₀+N.h)
X	/	x₀+N.h

Example: y" = ( x² - 1 ).y ; y(0) = 1 & y(-0.1) = 0.995012479

-Let's evaluate y(1)

-First, we load the required subroutine in main memory, for instance:

LBL "T" any global label, maximum 6 characters
X^2
1
-
*
RTN

"T" ASTO 00

     0   STO 01
     1   STO 02      0.995012479   STO 03
   0.1 STO 04                       10   STO 05     ( here, h = 0.1 and N = 10 )

XEQ "NUM1" gives x = 1 ( in 53 seconds ( or 48 seconds if you delete line 40 ) )
X<>Y y(1) = 0.606528753

-To continue with the same N-value, simply press R/S , it yields y(2) = 0.135332761

-The solution was actually y(x) = exp ( -x²/2 ) so y(1) = 0.606530660 and y(2) = 0.135335283 the error is of the order of 3 E-6

Notes:

-In a trajectory problem , we can use tabulated values to begin the calculations.
-If we don't have 2 starting values, the second may be found by a Runge-Kutta method, provided you know y'₀ .

-The loop stops when 2 consecutive approximations are equal ( line 41 ).
-This might lead to an infinite loop but practically, the risk is tiny because of the factor h² in the formula - provided h remains "reasonably small" however.

b) 2 differential equations y" = f(x,y,z) ; z" = g(x,y,z)

-Of course, the method may be generalized to systems of differential equations:

Data Registers: • R00 = subroutine name ( Registers R00 thru R07 are to be initialized before executing "NUM2" )

                                      • R01 = x₀
                                      • R02 = y₀   • R04 = h = stepsize                 • R06 = y_-1 = y(x₀ - h )
                                      • R03 = z₀   • R05 = N = number of steps    • R07 = z_-1 = z(x₀ - h )              R08 to R16: temp
Flags: /
Subroutine:   A program which computes f(x,y,z) in Y-register and g(x,y,z) in X-register
                       assuming x is in X-register, y is in Y-register and z is in Z-register upon entry.

01 LBL "NUM2"
02 RCL 05
03 STO 16
04 RCL 07
05 RCL 06
06 RCL 01
07 RCL 04
08 -
09 XEQ IND 00
10 STO 13
11 X<>Y
12 STO 10
13 RCL 03
14 STO 09
15 RCL 02
16 STO 08
17 RCL 01
18 XEQ IND 00
19 STO 14

20 X<>Y
21 STO 11
22 LBL 01
23 RCL 09
24 RCL 08
25 RCL 01
26 RCL 04
27 +
28 XEQ IND 00
29 STO 15
30 RCL 14
31 10
32 *
33 +
34 RCL 13
35 +
36 X<>Y
37 STO 12
38 RCL 11

39 10
40 *
41 +
42 RCL 10
43 +
44 RCL 04
45 X^2
46 12
47 /
48 ST* Z
49 *
50 RCL 06
51 -
52 RCL 02
53 ST+ X
54 +
55 ENTER^
56 X<> 08
57 -

58 ABS
59 X<>Y
60 RCL07
61 -
62 RCL 03
63 ST+ X
64 +
65 ENTER^
66 X<> 09
67 -
68 ABS
69 +
70 X#0?
71 GTO 01
72 RCL 08
73 X<> 02
74 STO 06
75 RCL 09
76 X<> 03

77 STO 07
78 RCL 12
79 X<> 11
80 STO 10
81 RCL 15
82 X<> 14
83 STO 13
84 RCL 04
85 ST+ 01
86 DSE 16
87 GTO 01
88 RCL 03
89 RCL 02
90 RCL 01
91 END

( 120 bytes / SIZE 017 )

STACK	INPUTS	OUTPUTS
Z	/	z(x₀+N.h)
Y	/	y(x₀+N.h)
X	/	x₀+N.h

Example:

y" = ( x - 2 ).z y(1) = 0.367879441 y(0.9) = 0.365912694
z" = y/x z(1) = 0.367879441 z(0.9) = 0.406569660

-Evaluate y(2) & z(2)

LBL "T"
ST/ Y
2
-
ST* Z
RDN
RTN

"T" ASTO 00
1 STO 01 0.367879441 STO 02 STO 03 h = 0.1 STO 04 , N = 10 STO 05

0.365912694 STO 06
0.406569660 STO 07

   XEQ "NUM2" >>>>       x = 2                        ( in 81 seconds )
                            RDN y(2) = 0.270670254
                            RDN z(2) = 0.135335322

-In fact, y = x exp(-x) and z = exp(-x) , so the exact results are y(2) = 0.270670566 & z(2) = 0.135335283

-Press R/S to continue the calculations ( after changing N in R05 if needed )

c) n differential equations y"₁ = f₁(x;y₁,.....,y_n) , ....... , y"_n = f_n(x;y₁,.....,y_n)

Data Registers: • R00 = subroutine name ( Registers R00 thru R02 and R09 thru R2n+10 are to be initialized before executing "NUM" )

• R01 = N = number of steps • R02 = h = stepsize R03 to R08: temp

                                      • R09 = n = numberof equations
                                      • R10 = x₀
                                      • R11 = y₁(x₀)               • R11+n = y₁(x₀-h)
                                      • R12 = y₂(x₀)               • R12+n = y₂(x₀-h)                         R11 thru R10+2n contain
                                        .....................                 ..........................                            the 2 starting values

• R10+n = y_n(x₀) • R10+2n = y_n(x₀-h)

-During the calculations,

                R11+n to R10+2n = f_i⁽¹⁾          the n values of the functions for x = x₀+h
                R11+2n to R10+3n = f_i⁽⁰⁾        the n values of the functions for x = x₀
                R11+3n to R10+4n = f_i^(-1)      the n values of the functions for x = x₀-h
                R11+4n to R10+5n = y_i⁽⁰⁾       the n values of y_i(x₀)
                R11+5n to R10+6n = y_i^(-1)      the n values of y_i(x₀-h)

Flags: /
Subroutine: A program which computes and stores f₁(x;y₁,.....,y_n) , ....... , f_n(x;y₁,.....,y_n) into R11+n , .......... , R10+2n respectively
with x ; y₁,.....,y_n in R10 ; R11 , ......... , R10+n respectively

-Lines 01 to 56 are only executed the first time in order to initialize some registers.

01 LBL "NUM"
02 RCL 09
03 4
04 *
05 11.011
06 STO 07
07 INT
08 +
09 STO 03
10 E3
11 /
12 RCL 07
13 INT
14 +
15 RCL 09
16 E6
17 /
18 STO 04
19 ST+ X
20 +
21 REGMOVE
22 XEQ IND 00
23 RCL 09
24 .1
25 %
26 ST+ X
27 +
28 RCL 07
29 +
30 RCL 04

31 +
32 STO 05
33 REGMOVE
34 RCL 03
35 RCL 07
36 FRC
37 +
38 RCL 04
39 +
40 STO 06
41 RCL 09
42 +
43 REGMOVE
44 RCL 02
45 ST- 10
46 XEQ IND 00
47 RCL 05
48 RCL 09
49 E3
50 /
51 +
52 REGMOVE
53 RCL 06
54 REGMOVE
55 RCL 02
56 ST+ 10
57 LBL 01
58 RCL 01
59 STO 08
60 LBL 02

61 RCL 02
62 ST+ 10
63 LBL 03
64 XEQ IND 00
65 RCL 09
66 RCL 09
67 RCL 09
68 10.01
69 +
70 STO 07
71 INT
72 +
73 STO 06
74 +
75 STO 05
76 +
77 STO 04
78 +
79 STO 03
80 CLX
81 LBL 04
82 RCL IND 04
83 RCL IND 05
84 10
85 *
86 +
87 RCL IND 06
88 +
89 RCL 02
90 X^2

91 *
92 12
93 /
94 RCL IND 03
95 ST+ X
96 +
97 RCL 03
98 RCL 09
99 +
100 RDN
101 RCL IND T
102 -
103 ENTER^
104 X<> IND 07
105 -
106 ABS
107 +
108 DSE 03
109 DSE 04
110 DSE 05
111 DSE 06
112 DSE 07
113 GTO 04
114 X#0?
115 GTO 03
116 RCL 05
117 1
118 +
119 .1
120 %

121 +
122 RCL 09
123 ST- Y
124 E6
125 /
126 STO 07
127 4
128 *
129 +
130 REGMOVE
131 RCL 03
132 1
133 +
134 E3
135 /
136 RCL 07
137 +
138 11
139 +
140 REGMOVE
141 DSE 08
142 GTO 02
143 RCL 13
144 RCL 12
145 RCL 11
146 RCL 10
147 RTN
148 GTO 01
149 END

( 210 bytes / SIZE 6n+11 )

STACK	INPUTS	OUTPUTS
T	/	y₃(x₀+N.h)
Z	/	y₂(x₀+N.h)
Y	/	y₁(x₀+N.h)
X	/	x₀+N.h

Example:

-We want to study the motion of a planet around a point sun ( the mass of the planet is neglected ).
-Here, the variable is t = time ( in days ), and y₁ , y₂ , y₃ are denoted x , y , z ( the coordinates are expressed in Astronomical Units )

      x" = -k² x/(x²+y²+z²)^3/2
      y" = -k² y/(x²+y²+z²)^3/2          where k = 0.01720209895 is Gauss' constant
      z" = -k² z/(x²+y²+z²)^3/2

-At t = 0

       x = 0.092
       y = -0.445
       z = -0.045

-At t = -1

       x = 0.070                               Find the position of the planet at t = 2 days
       y = -0.451
       z = -0.043

-We store -k² in register R30 17.20209895 E-3 X^2 CHS STO 30
-we load the following subroutine in main memory:

LBL "T"      X^2          RCL 13        SQRT        STO 14      RCL 29       RCL 30      STO 16
RCL 30      RCL 12    X^2              *                RCL 12       /                  *                RTN
RCL 11      X^2          +                  STO 29      RCL 30      STO 15       RCL 29
ST* Y        +               ENTER^      /                 *                 RCL 13       /

"T" ASTO 00 number of steps = 2 STO 01 stepsize h = 1 STO 02 there are 3 equations: 3 STO 09 then the initial values:

          0       STO 10
       0.092   STO 11      0.070    STO 14
     -0.445   STO 12     -0.451    STO 15
     -0.045   STO 13     -0.043    STO 16

     XEQ "NUM" >>>>   t =   2                                      ( execution time = 56 seconds )
                            RDN   x =   0.135070 = R11
                            RDN   y = -0.428856 = R12
                            RDN   z = -0.048573 = R13

N.B:

-To continue, press R/S or XEQ 01 ( after changing N in register R01 if needed )
but not XEQ "NUM" because the second new initial values are not in registers R14 thru R16 but in R26-R27-R28

-With the same N-value ( N = 2 ), you should find t = 4 , x = 0.176408 , y = -0.407227 , z = -0.051524

-For a planet like Mercury

with h = 1 day, the error is of the order of 2.7 10^-5 AU after 88 days
with h = 0.5 day, the error is of the order of 1.6 10^-6 AU after 88 days

-When h is divided by 2, the errors are approximately divided by 16 = 2⁴

2°) A formula of order 7

-We can use more than 2 starting values in order to reach a better accuracy.
-The following formula duplicates the Taylor's series up to the term in h⁷

y_n+1 = y_n + y_n-2 - y_n-3 + (h²/240).( 17. f_n+1 + 232.f_n + 222.f_n-1 + 232.f_n-2 + 17.f_n-3 ) [ h = the stepsize , f_k = f(x_k,y_k) , y_n-j = y(x_n-j.h) ]

-4 values are needed to initialize the algorithm: y₀ , y_-1 , y_-2 , y_-3

a) 1 differential equation y" = f(x,y)

Data Registers: • R00 = f name ( Registers R00 thru R07 are to be initialized before executing "7NUM1" )

• R01 = x₀
• R02 = y₀ • R03 = y_-1 = y(x₀ - h ) • R04 = h = stepsize • R05 = N = number of steps

                                                          • R06 = y_-2 = y(x₀ - 2.h )
                                                          • R07 = y_-3 = y(x₀ - 3.h )                                            R08 to R14: temp
Flags: /
Subroutine:   A program which computes f(x,y) assuming x is in X-register and y is in Y-register upon entry

01 LBL "7NUM1"
02 RCL 05
03 STO 08
04 RCL 07
05 RCL 01
06 RCL 04
07 3
08 *
09 -
10 XEQ IND 00
11 STO 09
12 RCL 06
13 RCL 01
14 RCL 04
15 ST+ X
16 -
17 XEQ IND 00

18 STO 10
19 RCL 03
20 RCL 01
21 RCL 04
22 -
23 XEQ IND 00
24 STO 11
25 RCL 02
26 STO 14
27 RCL 01
28 XEQ IND 00
29 STO 12
30 LBL 01
31 RCL 14
32 RCL 01
33 RCL 04
34 +

35 XEQ IND 00
36 STO 13
37 RCL 09
38 +
39 17
40 *
41 RCL 10
42 RCL 12
43 +
44 232
45 *
46 +
47 RCL 11
48 222
49 *
50 +
51 240

52 /
53 RCL 04
54 X^2
55 *
56 RCL 07
57 -
58 RCL 06
59 +
60 RCL 02
61 +
62 ENTER^
63 X<> 14
64 X#Y?
65 GTO 01
66 X<> 02
67 X<> 03
68 X<> 06

69 STO 07
70 RCL 13
71 X<> 12
72 X<> 11
73 X<> 10
74 STO 09
75 RCL 04
76 ST+ 01
77 DSE 08
78 GTO 01
79 RCL 02
80 RCL 01
81 END

( 115 bytes / SIZE 015 )

STACK	INPUTS	OUTPUTS
Y	/	y(x₀+N.h)
X	/	x₀+N.h

Example: y" = ( x² - 1 ).y Knowing y(0) = 1 ; y(-0.1) = 0.995012479 ; y(-0.2) = 0.980198673 ; y(-0.3) = 0.955997482

-Evaluate y(1)

   LBL "T"
   X^2
   1
   -
   *
   RTN

"T" ASTO 00

     0   STO 01
     1   STO 02       0.995012479   STO 03           h = 0.1 STO 04 , N = 10   STO 05
                             0.980198673   STO 06
                             0.955997482   STO 07

XEQ "7NUM1" gives x = 1 ( in 70 seconds )
X<>Y y(1) = 0.606530689

-To continue with the same N-value, simply press R/S , it yields y(2) = 0.135335319

-The errors are of the order of 3 E-8

b) 2 differential equations y" = f(x,y,z) ; z" = g(x,y,z)

Data Registers: • R00 = subroutine name ( Registers R00 thru R11 are to be initialized before executing "7NUM2" )

                                      • R01 = x₀
                                      • R02 = y₀   • R04 = h = stepsize                 • R06 = y_-1 = y(x₀ - h )          • R09 = z_-1 = z(x₀ - h )
                                      • R03 = z₀   • R05 = N = number of steps    • R07 = y_-2 = y(x₀ - 2h )        • R10 = z_-2 = z(x₀ - 2h )           R12 to R24: temp
                                                                                                            • R08 = y_-3 = y(x₀ - 3h )        • R11 = z_-3 = z(x₀ - 3h )
Flags: /
Subroutine:   A program which computes f(x,y,z) in Y-register and g(x,y,z) in X-register
                           assuming x is in X-register, y is in Y-register and z is in Z-register upon entry.

-Line 135 is a three-byte GTO 01

01 LBL "7NUM2"
02 RCL 05
03 STO 24
04 RCL 11
05 RCL 08
06 RCL 04
07 3
08 *
09 RCL 01
10 X<>Y
11 -
12 XEQ IND 00
13 STO 17
14 X<>Y
15 STO 12
16 RCL 10
17 RCL 07
18 RCL 01
19 RCL 04
20 ST+ X
21 -
22 XEQ IND 00
23 STO 18
24 X<>Y
25 STO 13
26 RCL 09
27 RCL 06
28 RCL 01

29 RCL 04
30 -
31 XEQ IND 00
32 STO 19
33 X<>Y
34 STO 14
35 RCL 03
36 STO 23
37 RCL 02
38 STO 22
39 RCL 01
40 XEQ IND 00
41 STO 20
42 X<>Y
43 STO 15
44 LBL 01
45 RCL 23
46 RCL 22
47 RCL 01
48 RCL 04
49 +
50 XEQ IND 00
51 STO 21
52 RCL 17
53 +
54 17
55 *
56 RCL 18

57 RCL 20
58 +
59 232
60 *
61 +
62 RCL 19
63 222
64 *
65 +
66 X<>Y
67 STO 16
68 RCL 12
69 +
70 17
71 *
72 RCL 13
73 RCL 15
74 +
75 232
76 *
77 +
78 RCL 14
79 222
80 *
81 +
82 RCL 04
83 X^2
84 240

85 /
86 ST* Z
87 *
88 RCL 08
89 -
90 RCL 07
91 +
92 RCL 02
93 +
94 ENTER^
95 X<> 22
96 -
97 ABS
98 X<>Y
99 RCL 11
100 -
101 RCL 10
102 +
103 RCL 03
104 +
105 ENTER^
106 X<> 23
107 -
108 ABS
109 +
110 X#0?
111 GTO 01
112 RCL 22

113 X<> 02
114 X<> 06
115 X<> 07
116 STO 08
117 RCL 23
118 X<> 03
119 X<> 09
120 X<> 10
121 STO 11
122 RCL 16
123 X<> 15
124 X<> 14
125 X<> 13
126 STO 12
127 RCL 21
128 X<> 20
129 X<> 19
130 X<> 18
131 STO 17
132 RCL 04
133 ST+ 01
134 DSE 24
135 GTO 01
136 RCL 03
137 RCL 02
138 RCL 01
139 END

( 207 bytes / SIZE 025 )

STACK	INPUTS	OUTPUTS
Z	/	z(x₀+N.h)
Y	/	y(x₀+N.h)
X	/	x₀+N.h

Example:

y" = ( x - 2 ).z y(1) = 0.367879441 y(0.9) = 0.365912694 y(0.8) = 0.359463171 y(0.7) = 0.347609713
z" = y/x z(1) = 0.367879441 z(0.9) = 0.406569660 z(0.8) = 0.449328964 z(0.7) = 0.496585304

-Evaluate y(2) & z(2)

   LBL "T"
   ST/ Y
   2
   -
   ST* Z
   RDN
   RTN

"T" ASTO 00
1 STO 01 0.367879441 STO 02 STO 03 h = 0.1 STO 04 , N = 10 STO 05

   0.365912694   STO 06           0.406569660   STO 09
   0.359463171   STO 07           0.449328964   STO 10
   0.347609713   STO 08           0.496585304   STO 11

   XEQ "7NUM2" >>>>     x = 2                        ( in 125 seconds )
                            RDN y(2) = 0.270670563
                            RDN z(2) = 0.135335281

-The errors are of the order of 3 E-9

-Press R/S to continue the calculations ( after changing N in R05 if needed )

c) n differential equations y"₁ = f₁(x;y₁,.....,y_n) , ....... , y"_n = f_n(x;y₁,.....,y_n)

Data Registers: • R00 = subroutine name ( Registers R00 thru R02 and R13 thru R4n+14 are to be initialized before executing "7NUM" )

• R01 = N = number of steps • R02 = h = stepsize R03 to R12: temp

                                      • R13 = n = numberof equations
                                      • R14 = x₀
                                      • R15 = y₁(x₀)               • R15+n = y₁(x₀-h)               • R15+2n = y₁(x₀-2h)               • R15+3n = y₁(x₀-3h)
                                      • R16 = y₂(x₀)               • R16+n = y₂(x₀-h)               • R16+2n = y₂(x₀-2h)               • R16+3n = y₂(x₀-3h)      R15 thru R14+4n =
                                        .....................                 ..........................                   ..............................

• R14+n = y_n(x₀) • R14+2n = y_n(x₀-h) • R14+3n = y_n(x₀-2h) • R14+4n = y_n(x₀-3h) the 4 starting values

-During the calculations,

            R15+n to R14+2n = f_i⁽¹⁾          the n values of the functions for x = x₀+h
            R15+2n to R14+3n = f_i⁽⁰⁾        the n values of the functions for x = x₀
            R15+3n to R14+4n = f_i^(-1)       the n values of the functions for x = x₀-h
            R15+4n to R14+5n = f_i^(-2)       the n values of the functions for x = x₀-2h
            R15+5n to R14+6n = f_i^(-3)       the n values of the functions for x = x₀-3h
            R15+6n to R14+7n = y_i⁽⁰⁾       the n values of y_i(x₀)
            R15+7n to R14+8n = y_i^(-1)      the n values of y_i(x₀-h)
            R15+8n to R14+9n = y_i^(-2)      the n values of y_i(x₀-2h)
            R15+9n to R14+10n = y_i^(-3)    the n values of y_i(x₀-3h)

Flags: /
Subroutine: A program which computes and stores f₁(x;y₁,.....,y_n) , ....... , f_n(x;y₁,.....,y_n) into R15+n , .......... , R14+2n respectively
with x ; y₁,.....,y_n in R14 ; R15 , ......... , R14+n respectively

-Lines 01 to 62 are only executed the first time in order to initialize some registers.
-Lines 166 & 172 are three-byte GTOs

01 LBL "7NUM"
02 RCL 14
03 STO 10
04 RCL 13
05 E6
06 /
07 STO 03
08 4
09 *
10 RCL 13
11 E3
12 /
13 STO 04
14 6
15 *
16 +
17 15.015
18 STO 05
19 +
20 REGMOVE
21 LASTX
22 RCL 03
23 +
24 RCL 04
25 +
26 RCL 13
27 +
28 STO 06
29 RCL 03
30 RCL 05
31 +
32 RCL 13
33 6
34 *
35 +

36 STO 07
37 STO 08
38 4
39 STO 09
40 LBL 10
41 XEQ IND 00
42 RCL 04
43 RCL 06
44 +
45 STO 06
46 REGMOVE
47 DSE 09
48 X=0?
49 GTO 00
50 RCL 07
51 RCL 13
52 +
53 STO 07
54 REGMOVE
55 RCL 02
56 ST- 14
57 GTO 10
58 LBL 00
59 RCL 08
60 REGMOVE
61 RCL 10
62 STO 14
63 LBL 01
64 RCL 01
65 STO 12
66 LBL 02
67 RCL 02
68 ST+ 14
69 LBL 03
70 XEQ IND 00

71 RCL 13
72 RCL 13
73 RCL 13
74 14.014
75 +
76 STO 03
77 INT
78 +
79 STO 04
80 +
81 STO 05
82 +
83 STO 06
84 +
85 STO 07
86 +
87 STO 08
88 +
89 STO 09
90 +
91 +
92 STO 10
93 +
94 STO 11
95 CLX
96 LBL 04
97 RCL IND 04
98 RCL IND 08
99 +
100 17
101 *
102 RCL IND 05
103 RCL IND 07
104 +
105 232

106 *
107 +
108 RCL IND 06
109 222
110 *
111 +
112 RCL 02
113 X^2
114 *
115 240
116 /
117 RCL IND 09
118 +
119 RCL IND 10
120 +
121 RCL IND 11
122 -
123 ENTER^
124 X<> IND 03
125 -
126 ABS
127 +
128 DSE 11
129 DSE 10
130 DSE 09
131 DSE 08
132 DSE 07
133 DSE 06
134 DSE 05
135 DSE 04
136 DSE 03
137 GTO 04
138 X#0?
139 GTO 03
140 RCL 05

141 1
142 +
143 .1
144 %
145 +
146 RCL 13
147 ST- Y
148 E6
149 /
150 STO 07
151 8
152 *
153 +
154 REGMOVE
155 RCL 09
156 1
157 +
158 E3
159 /
160 RCL 07
161 +
162 15
163 +
164 REGMOVE
165 DSE 12
166 GTO 02
167 RCL 17
168 RCL 16
169 RCL 15
170 RCL 14
171 RTN
172 GTO 01
173 END

( 246 bytes / SIZE 10n+15 )

STACK	INPUTS	OUTPUTS
T	/	y₃(x₀+N.h)
Z	/	y₂(x₀+N.h)
Y	/	y₁(x₀+N.h)
X	/	x₀+N.h

Example:

-We study again the motion of a planet around a point sun ( the mass of the planet is neglected ).
-The variable is t = time ( in days ), and y₁ , y₂ , y₃ are denoted x , y , z ( the coordinates are expressed in Astronomical Units )

      x" = -k² x/(x²+y²+z²)^3/2
      y" = -k² y/(x²+y²+z²)^3/2          where k = 0.01720209895 is Gauss' constant
      z" = -k² z/(x²+y²+z²)^3/2

t = 0 t = -1 t = -2 t = -3

       x = 0.293510249          x = 0.301200207          x = 0.305864609           x = 0.307427938
       y = 0.091967806          y = 0.061830391          y = 0.031072548          y = 0
       z = 0.040946705           z = 0.027528664          z = 0.013834390          z = 0

Find the position of the planet at t = 4 days

-We store -k² in register R45 17.20209895 E-3 X^2 CHS STO 45
and we load the following subroutine in main memory:

LBL "T"      X^2          RCL 17        SQRT        STO 18      RCL 46       RCL 45      STO 20
RCL 45      RCL 16    X^2              *                RCL 16       /                  *                RTN
RCL 15      X^2          +                  STO 46      RCL 45      STO 19       RCL 46
ST* Y        +               ENTER^      /                 *                 RCL 17       /

"T" ASTO 00 number of steps = 4 STO 01 stepsize h = 1 STO 02 there are 3 equations: 3 STO 13 then the initial values:

               0            STO 14
     0.293510249   STO 15      0.301200207    STO 18    0.305864609     STO 21   0.307427938   STO 24
     0.091967806   STO 16      0.061830391    STO 19    0.031072548     STO 22             0            STO 25
     0.040946705   STO 17      0.027528664    STO 20    0.013834390     STO 23             0            STO 26

     XEQ "7NUM" >>>   t = 4                     = R14                    ( execution time = 2m32s )
                            RDN   x = 0.235500989 = R15
                            RDN   y = 0.200940664 = R16
                            RDN   z = 0.089464547 = R17

N.B:

-To continue, press R/S or XEQ 01 ( after changing N in register R01 if needed ) but not XEQ "7NUM"

-For a planet like Mercury,

with h = 1 day, the error is of the order of 3.6 10^-7 AU after 88 days
with h = 0.5 day, the error is of the order of 5.8 10^-9 AU after 88 days

-When h is divided by 2, the errors are approximately divided by 64 = 2⁶

-In fact, the error - after 88 days with h = 0.5 day - is 5.8 10^-9 AU on an HP-48,
but it reaches 1.6 10^-8 AU on an HP-41 because of roundoff errors: the HP-41 works with 10 digits only!

Reference:

[1] Francis Scheid "Numerical Analysis" McGraw-Hill ISBN 07-055197-9

hp41programs

Numerov's Method for the HP-41