The Museum of HP Calculators
The Museum of HP Calculators
This program is Copyright © 2007 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.
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"1 = f1(x;y1,.....,yn)
, ....... , y"n = fn(x;y1,.....,yn)
( X-Functions Module required )
2°) A formula of order 7
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"1 = f1(x;y1,.....,yn)
, ....... , y"n = fn(x;y1,.....,yn)
( X-Functions Module required )
-Numerov's method has been devised to solve second order differential
equations where the first derivative y' does not appear explicitly.
-The formula is: yn+1 = 2.yn - yn-1
+ (h2/12).( fn+1 + 10.fn + fn-1
) where h = xn - xn-1 is the stepsize
and fk = f(xk,yk) ;
yn-j = y(xn-j.h)
-It requires 2 starting values y-1 and y0
instead of y0 and y'0 .
-The method duplicates the Taylor series up to the term in h5
-Since it is an implicit formula ( yn+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 = x0
• R02 = y0 • R03 = y-1 =
y(x0 - 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
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
this line is only useful to display the successive approximations, but
it's not necessary for the calculations.
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(x0+N.h) |
| X | / | x0+N.h |
Example: y" = ( x2 - 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 ( -x2/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'0 .
-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 h2 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 = x0
• R02 = y0 • R04 = h = stepsize
• R06 = y-1 = y(x0 - h )
• R03 = z0 • R05 = N = number of steps
• R07 = z-1 = z(x0 - 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(x0+N.h) |
| Y | / | y(x0+N.h) |
| X | / | x0+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"1 = f1(x;y1,.....,yn)
, ....... , y"n = fn(x;y1,.....,yn)
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 = number of equations
• R10 = x0
• R11 = y1(x0)
• R11+n = y1(x0-h)
• R12 = y2(x0)
• R12+n = y2(x0-h)
R11 thru R10+2n contain
.....................
..........................
the 2 starting values
• R10+n = yn(x0)
• R10+2n = yn(x0-h)
-During the calculations, R11+n to R10+2n = fi(1)
the n values of the functions for x = x0+h
R11+2n to R10+3n = fi(0)
------------------------------- x = x0
R11+3n to R10+4n = fi(-1)
------------------------------- x = x0-h
R11+4n to R10+5n = yi(0)
--------------- yi(x0)
R11+5n to R10+6n = yi(-1)
--------------- yi(x0-h)
Flags: /
Subroutine: A program which computes
and
stores f1(x;y1,.....,yn)
, ....... , fn(x;y1,.....,yn) into
R11+n , .......... , R10+2n respectively
with x ; y1,.....,yn in R10 ; R11 , .........
, R10+n respectively
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
Lines 01 to 56 are only executed the first time in order to initialize
some registers.
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
the loop stops when 2 consecutive approximations are equal
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 | / | y3(x0+N.h) |
| Z | / | y2(x0+N.h) |
| Y | / | y1(x0+N.h) |
| X | / | x0+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 y1 , y2
, y3 are denoted x , y , z ( the coordinates
are expressed in Astronomical Units )
x" = -k2 x/(x2+y2+z2)3/2
-The equations are y" = -k2
y/(x2+y2+z2)3/2
where k = 0.01720209895 is Gauss' constant
z" = -k2 z/(x2+y2+z2)3/2
-At t = 0 x = 0.092
at t = -1 x = 0.070
Find the position of the planet at t = 2 days
y = -0.445
y = -0.451
z = -0.045
z = -0.043
-We store -k2 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
---- h = 0.5 day, -------------------------- 1.6
10 -6 AU -------------
-When h is divided by 2, the errors are approximately divided by 16
= 24
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 h7
yn+1 = yn + yn-2 - yn-3 + (h2/240).( 17. fn+1 + 232.fn + 222.fn-1 + 232.fn-2 + 17.fn-3 ) [ h = the stepsize , fk = f(xk,yk) , yn-j = y(xn-j.h) ]
-4 values are needed to initialize the algorithm: y0
, 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 = x0
• R02 = y0 • R03 = y-1 =
y(x0 - h ) • R04 =
h = stepsize • R05 = N = number of steps
• R06 = y-2 = y(x0 - 2.h )
• R07 = y-3 = y(x0 - 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(x0+N.h) |
| X | / | x0+N.h |
Example: y" = ( x2 - 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 = x0
• R02 = y0 • R04 = h = stepsize
• R06 = y-1 = y(x0 - h )
• R09 = z-1 = z(x0 - h )
• R03 = z0 • R05 = N = number of steps
• R07 = y-2 = y(x0 - 2h )
• R10 = z-2 = z(x0 - 2h )
R12 to R24: temp
• R08 = y-3 = y(x0 - 3h )
• R11 = z-3 = z(x0 - 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
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
a three-byte GTO
136 RCL 03
137 RCL 02
138 RCL 01
139 END
( 207 bytes / SIZE 025 )
| STACK | INPUTS | OUTPUTS |
| Z | / | z(x0+N.h) |
| Y | / | y(x0+N.h) |
| X | / | x0+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"1 = f1(x;y1,.....,yn)
, ....... , y"n = fn(x;y1,.....,yn)
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 = number of equations
• R14 = x0
• R15 = y1(x0)
• R15+n = y1(x0-h)
• R15+2n = y1(x0-2h)
• R15+3n = y1(x0-3h)
• R16 = y2(x0)
• R16+n = y2(x0-h)
• R16+2n = y2(x0-2h)
• R16+3n = y2(x0-3h)
R15 thru R14+4n
.....................
..........................
..............................
contain
• R14+n = yn(x0)
• R14+2n = yn(x0-h)
• R14+3n = yn(x0-2h)
• R14+4n = yn(x0-3h)
the 4 starting values
-During the calculations, R15+n to R14+2n = fi(1)
the n values of the functions for x = x0+h
R15+2n to R14+3n = fi(0)
------------------------------- x = x0
R15+3n to R14+4n = fi(-1)
------------------------------- x = x0-h
R15+4n to R14+5n = fi(-2)
------------------------------- x = x0-2h
R15+5n to R14+6n = fi(-3)
------------------------------- x = x0-3h
R15+6n to R14+7n = yi(0)
--------------- yi(x0)
R15+7n to R14+8n = yi(-1)
--------------- yi(x0-h)
R15+8n to R14+9n = yi(-2)
--------------- yi(x0-2h)
R15+9n to R14+10n = yi(-3) ---------------
yi(x0-3h)
Flags: /
Subroutine: A program which computes
and
stores f1(x;y1,.....,yn)
, ....... , fn(x;y1,.....,yn) into
R15+n , .......... , R14+2n respectively
with x ; y1,.....,yn in R14 ; R15 , .........
, R14+n respectively
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
Lines 01 to 62 are only executed the first time in order to initialize
some registers.
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
a three-byte GTO
167 RCL 17
168 RCL 16
169 RCL 15
170 RCL 14
171 RTN
172 GTO 01
a three-byte GTO
173 END
( 246 bytes / SIZE 10n+15 )
| STACK | INPUTS | OUTPUTS |
| T | / | y3(x0+N.h) |
| Z | / | y2(x0+N.h) |
| Y | / | y1(x0+N.h) |
| X | / | x0+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 y1 , y2
, y3 are denoted x , y , z ( the coordinates
are expressed in Astronomical Units )
x" = -k2 x/(x2+y2+z2)3/2
-The equations are y" = -k2
y/(x2+y2+z2)3/2
where k = 0.01720209895 is Gauss' constant
z" = -k2 z/(x2+y2+z2)3/2
-At t = 0 x = 0.293510249
at t = -1 x = 0.301200207
at t = -2 x = 0.305864609
at t = -3 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 -k2 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
---- h = 0.5 day, -------------------------- 5.8
10 -9 AU -------------
-When h is divided by 2, the errors are approximately divided by 64 = 26
-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
Go back to the HP-41 software library
Go back to the general software library
Go
back to the main exhibit hall