The Museum of
HP CalculatorsThis program is 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°) An Implicit Runge-Kutta Method with order
8
2°) A General Implicit-Runge-Kutta Program
( with an example of a 12th-order method )
3°) A Slighly less general Implicit-Runge-Kutta
Program ( with an example of a 13th-order method )
-These methods are applied to solve:
a) A first order differential equation:
dy/dx = f(x,y)
with the initial value: y(x0) = y0
b) A system of 2 differential equations:
dy/dx = f(x,y,z) ; dz/dx = g(x,y,z)
-------------------s: y(x0) = y0
; z(x0) = z0
1°) An implicit Runge-Kutta Method with order 8
a) First order differential
equations: dy/dx = f (x;y)
-The formulae are called explicit Runge-Kutta methods if the
successive k1 , k2 , .... , kn may
be directly computed.
-In the following program , we have k1 = h.f ( x +
b1h ; y + a1,1 k1 + a1,2
k2 + ..... + a1,5 k5 )
.....................................................................................
where ai,j # 0 for j > i-1
k5 = h.f ( x + b5h ; y + a5,1 k1
+ a5,2 k2 + ..... + a5,5 k5
)
and ci = a5,i
-The advantage is that only 5 stages ( instead of 11 ) are required
to obtain an 8th-order method: the program is shorter and less data registers
are needed.
-But we have to solve a 5x5 non-linear system within each step! Speed
is not a characteristic of implicit methods, especially with an HP-41...
-"IRK8" uses an iterative algorithm based on the "fixed-point theorem"
-Some implicit Runge-Kutta formulae are also satisfactory for stiff
problems, but the elementary iterative method used by "IRK8" will seldom
lead to convergence
in such a case, unless h is very small: another algorithm
( like the Newton's method ) would be better to solve the 5x5 non-linear
system, if you're in the mood...
Data Registers: ( Registers R00 to R04 and R16 to R45 are to be initialized before executing "IRK8" )
• R00: f name • R01
= x0 • R03 = h = stepsize
R05 to R10: temp ; R11 = k1 ......
R15 = k5
• R02 = y0 •
R04 = N = number of steps •
R16 thru R45 = the 30 following numbers:
R16 = 1
R26 = 29/180
R36 = 29/180
R17 = 1/20
R27 = -3/140
R37 = -0.06901154103
the decimal numbers are actually
R18 = 49/180
R28 = 1/2
R38 = 0.05200216599
rational functions of sqrt(21)
R19 = 16/45
R29 = 1/20
R39 = -3/140
( all decimals are correct )
R20 = 49/180
R30 = 0.2813091833
R40 = 0
R21 = 1/20
R31 = 73/360
R41 = 1/20
R22 = 0.8273268354
R32 = -0.05283696110
R42 = -7/60
R23 = 1/20
R33 = 3/160
R43 = 2/15
R24 = 0.2702200562
R34 = 0.1726731646
R44 = -7/60
R25 = 0.3674242394
R35 = 1/20
R45 = 1/20
Flags: /
Subroutine: a program calculating f(x;y)
assuming x is in X-register and y is in Y-register upon entry.
01 LBL "IRK8"
02 RCL 04
03 STO 09
04 CLX
05 STO 11
06 STO 12
07 STO 13
08 STO 14
09 STO 15
10 LBL 01
11 11.015
12 STO 05
13 45
14 STO 07
15 CLX
16 STO 08
17 LBL 02
18 15.01
19 STO 06
20 CLX
21 LBL 03
22 RCL IND 06
23 RCL IND 07
24 *
25 +
26 DSE 07
27 DSE 06
28 GTO 03
29 RCL 02
30 +
31 STO 10
32 RCL IND 07
33 RCL 03
34 *
35 RCL 01
36 +
37 XEQ IND 00
38 RCL 03
39 *
40 ENTER^
41 X<> IND 05
42 -
43 ABS
44 ST+ 08
45 DSE 07
46 ISG 05
47 GTO 02
48 RCL 08
49 VIEW X
50 E-9
or another "small" number:
51 X<Y?
52 GTO 01
53 RCL 03
54 ST+ 01
55 RCL 10
56 STO 02
57 DSE 09
58 GTO 01
59 RCL 01
60 CLD
61 END
( 99 bytes / SIZE 046 )
| STACK | INPUTS | OUTPUTS |
| Y | / | y(x0+N.h) |
| X | / | x0+N.h |
Example: y(0) = 1 and dy/dx = -2 x.y Evaluate y(0.5)
-Initialize registers R16 thru R45.
01 LBL "T"
02 *
03 ST+ X
04 CHS
05 RTN
"T" ASTO 00
0 STO 01
1 STO 02 and if we choose
h = 0.1
0.1 STO 03
5 STO 04 XEQ "IRK8"
yields x = 0.5
( in 5mn57 )
X<>Y
y(0.5) = 0.7788007830 ( exact result
is 0.7788007831 )
-If you replace lines 50-51 with X#0? , you'll get the exact value 0.7788007831 but this severe test might sometimes lead to an infinite loop!
-The HP-41 displays | k1 - k'1 | + ............
+ | k5 - k'5 | where ki
and k'i are 2 successive k-approximations:
-At each step, this sum must tend to zero. Otherwise, reduce the stepsize
h and start over!
-Press R/S to continue with the same h and N.
-Start over with a smaller stepsize to obtain an error estimate.
b) Systems of 2 first
order differential equations: dy/dx = f (x;y;z) ;
dz/dx = g(x;y;z)
-Now, we have to solve a 10x10 non-linear system within each step.
Data Registers: ( Registers R00 to R05 and R25 to R54 are to be initialized before executing "IRK8B" )
• R00: f name • R01
= x0 • R04 = h = stepsize
R06 to R14: temp ; R15 , ..... , R24 = ki
• R02 = y0 •
R05 = N = number of steps •
R25 thru R54 = the same 30 numbers used in the previous program.
• R03 = z0
Flags: /
Subroutine: a program calculating 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 "IRK8B"
02 RCL 05
03 STO 12
04 15.024
05 CLRGX
if you don't have an HP-41CX , replace lines 04-05 with CLX
STO 15 STO 16 STO 17 STO 18 STO 19
06 LBL 01
STO 20 STO 21 STO 22 STO 23 STO 24
07 15.019
08 STO 06
09 20
10 STO 07
11 54
12 STO 10
13 CLX
14 STO 11
15 LBL 02
16 19.014
17 STO 08
18 24
19 STO 09
20 CLST
21 LBL 03
22 X<>Y
23 RCL IND 08
24 RCL IND 10
25 *
26 +
27 X<>Y
28 RCL IND 09
29 RCL IND 10
30 *
31 +
32 DSE 10
33 DSE 09
34 DSE 08
35 GTO 03
36 RCL 03
37 +
38 STO 14
39 X<>Y
40 RCL 02
41 +
42 STO 13
43 RCL IND 10
44 RCL 04
45 *
46 RCL 01
47 +
48 XEQ IND 00
49 RCL 04
50 ST* Z
51 *
52 ENTER^
53 X<> IND 07
54 -
55 ABS
56 X<>Y
57 ENTER^
58 X<> IND 06
59 -
60 ABS
61 +
62 ST+ 11
63 SIGN
64 ST+ 07
65 ST- 10
66 ISG 06
67 GTO 02
68 RCL 11
69 VIEW X
70 E-9
or another "small" number
71 X<Y?
72 GTO 01
73 RCL 04
74 ST+ 01
75 RCL 14
76 STO 03
77 RCL 13
78 STO 02
79 DSE 12
80 GTO 01 a three-byte
GTO
81 RCL 01
82 CLD
83 END
( 138 bytes / SIZE 055 )
| STACK | INPUTS | OUTPUTS |
| Z | / | z(x0+N.h) |
| Y | / | y(x0+N.h) |
| X | / | x0+N.h |
Example: y(0) = 1 ; z(0) = 0 ; dy/dx = z , dz/dx = -2 x.z -2 y ; calculate y(0.5) and z(0.5)
-Initialize registers R25 thru R54.
01 LBL "U"
02 RCL Z
03 *
04 +
05 ST+ X
06 CHS
07 RTN
"U" ASTO 00
0 STO 01
1 STO 02
0 STO 03 and if we choose
h = 0.1
0.1 STO 04
5 STO 05 XEQ "IRK8B"
yields x = 0.5
execution time = 12mn
RDN
y(0.5) = 0.7788007830
( exact result is 0.7788007831 )
RDN
z(0.5) = -0.7788007830 ( exact
result is -0.7788007831 )
-The HP-41 displays | k1 - k'1 | + ............
+ | k10 - k'10 | where ki
and k'i are 2 successive k-approximations:
-At each step, this sum must tend to zero. Otherwise, reduce the stepsize
h and start over!
-Press R/S to continue with the same h and N.
-Start over with a smaller stepsize to obtain an error estimate.
2°) A General implicit-Runge-Kutta Program
a) First order differential
equations: dy/dx = f (x;y)
-A n-stage implicit Runge-Kutta method is defined by n(n+2) coefficients ai,j ; bi ; ci
k1 = h.f ( x + b1h ; y + a1,1 k1
+ a1,2 k2 + ..... + a1,n kn
)
.....................................................................................
( actually, ai,1 + ai,2
+ .......... + ai,n = bi )
kn = h.f ( x + bnh ; y + an,1 k1
+ an,2 k2 + ..... + an,n kn
)
then: y(x+h) = y(x) + c1.k1+ ................ + cn.kn
-The following program will work for any implicit method.
Data Registers: ( Registers R00 to R05 and Rn+11 to Rn2+3n+10 are to be initialized before executing "IRK" )
• R00: f name • R01
= x0
R06 to R10: temp
• Rn+11 to Rn2+3n+10 = the ( n2 + 2n
) coefficients ai,j ; bi ; ci
• R02 = y0
R11 = k1
which are to be stored as shown below:
• R03 = h = stepsize
R12 = k2
• R04 = N = number of steps
............
• R05 = n = number of stages
Rn+10 = kn
• Rn+11 = a1,1
• Rn+12 = a1,2 ...........................
• R2n+10 = a1,n • R2n+11
= b1
• R2n+12 = a2,1
• R2n+13 = a2,2 ...........................
• R3n+11 = a2,n • R3n+12
= b2
..................................................................................................................................................
• Rn2+n+10 = an,1 • Rn2+n+11 = an,2 ....................... • Rn2+2n+9 = an,n • Rn2+2n+10 = bn
• Rn2+2n+11 = c1 • Rn2+2n+12 = c2 ....................... • Rn2+3n+10 = cn
Flags: /
Subroutine: a program calculating f(x;y)
assuming x is in X-register and y is in Y-register upon entry.
>>>>> In
other words, this subroutine has to change the stack from Y
= y
X = x into
X' = f(x;y)
01 LBL "IRK"
02 RCL 04
03 STO 06
04 RCL 05
05 10
06 +
07 .1
08 %
09 11
10 +
11 STO 07
12 CLRGX
if you don't have an HP-41CX , replace line 12 by the 5 lines:
0 LBL 00 STO IND Y
ISG Y GTO 00
13 LBL 01
14 RCL 05
15 STO 08
16 RCL 07
17 FRC
18 11.9
19 ST+ 08
20 INT
21 +
22 STO 07
23 CLX
24 STO 10
25 LBL 02
26 RCL 07
27 FRC
28 11
29 +
30 STO 09
31 CLX
32 LBL 03
33 RCL IND 08
34 RCL IND 09
35 *
36 +
37 ISG 08
38 ISG 09
39 GTO 03
40 RCL 02
41 +
42 RCL 03
43 RCL IND 08
44 *
45 RCL 01
46 +
47 XEQ IND 00
48 RCL 03
49 *
50 ENTER^
51 X<> IND 07
52 -
53 ABS
54 ST+ 10
55 ISG 08
56 ISG 07
57 GTO 02
58 RCL 10
59 VIEW X
60 E-9
or another "small" number
61 X<Y?
62 GTO 01
63 RCL 05
64 ST- 07
65 CLX
66 LBL 04
67 RCL IND 07
68 RCL IND 08
69 *
70 +
71 ISG 08
72 ISG 07
73 GTO 04
74 ST+ 02
75 RCL 03
76 ST+ 01
77 DSE 06
78 GTO 01
79 RCL 02
80 RCL 01
81 CLD
82 END
( 126 bytes / SIZE n2+3n+11 )
| STACK | INPUTS | OUTPUTS |
| Y | / | y(x0+N.h) |
| X | / | x0+N.h |
-Implicit methods are very accurate: n-stage 2n-order methods
do exist which would be impossible with explicit methods,
but, on the other hand, we have to solve a nxn non-linear system
within each step!
-For instance, here are the 48 coefficients of a 6-stage 12-order method,
based on Gauss-Legendre quadrature.
-These coefficients are given with 20 decimals in case you want to
use them with a calculator working with more than 10 digits:
R17 = a1,1
= 0.04283 11230 94792 58626
R24 = a2,1 = 0.09267 34914 30378 86319
R31 = a3,1 = 0.08224 79226 12843 87381
R18 = a1,2
= -0.01476 37259 97197 41248
R25 = a2,2 = 0.09019 03932 62034 65189
R32 = a3,2 = 0.19603 21623 33245 00606
R19 = a1,3
= 0.00932 50507 06477 75119
R26 = a2,3 = -0.02030 01022 93239 58595
R33 = a3,3 = 0.11697 84836 43172 76185
R20 = a1,4
= -0.00566 88580 49483 51191
R27 = a2,4 = 0.01036 31562 40246 42374
R34 = a3,4 = -0.02048 25277 45656 09763
R21 = a1,5
= 0.00285 44333 15099 33514
R28 = a2,5 = -0.00488 71929 28037 67146
R35 = a3,5 = 0.00798 99918 99662 33580
R22 = a1,6
= -0.00081 27801 71264 76211
R29 = a2,6 = 0.00135 55610 55485 06178
R36 = a3,6 = -0.00207 56257 84866 33419
R23 = b1
= 0.03376 52428 98423 98609
R30 = b2 = 0.16939 53067 66867 74317
R37 = b3 = 0.38069 04069 58401
54568
R38 = a4,1
= 0.08773 78719 74451 50671
R45 = a5,1 = 0.08430 66851 34100 11074
R52 = a6,1 = 0.08647 50263 60849 93463
R39 = a4,2
= 0.17239 07946 24406 96799
R46 = a5,2 = 0.18526 79794 52106 97525
R53 = a6,2 = 0.17752 63532 08969 96865
R40 = a4,3
= 0.25443 94950 32001 62132
R47 = a5,3 = 0.22359 38110 46099 09996
R54 = a6,3 = 0.23962 58253 35829 03560
R41 = a4,4
= 0.11697 84836 43172 76185
R48 = a5,4 = 0.25425 70695 79585 10965
R55 = a6,4 = 0.22463 19165 79867 77250
R42 = a4,5
= -0.01565 13758 09175 70227
R49 = a5,5 = 0.09019 03932 62034 65189
R56 = a6,5 = 0.19514 45125 21266 71626
R43 = a4,6
= 0.00341 43235 76741 29871
R50 = a5,6 = -0.00701 12452 40793 69067
R57 = a6,6 = 0.04283 11230 94792 58626
R44 = b4
= 0.61930 95930 41598 45432
R51 = b5 = 0.83060 46932 33132 25683
R58 = b6 = 0.96623 47571 01576 01391
R59 = c1
= 0.08566 22461 89585 17252
R61 = c3 = 0.23395 69672 86345 52369
R63 = c5 = 0.18038 07865 24069 30378
R60 = c2
= 0.18038 07865 24069 30378
R62 = c4 = 0.23395 69672 86345 52369
R64 = c6 = 0.08566 22461 89585 17252
Example: y(0) = 1 and dy/dx
= -2 x.y Evaluate y(1)
01 LBL "T"
02 *
03 ST+ X
04 CHS
05 RTN
-Initialize registers R17 thru R64
"T" ASTO 00
0 STO 01
1 STO 02 and if we choose
h = 0.5
0.5 STO 03
2 STO 04 ( the number of
steps )
6 STO 05 ( we are using a
6-stage method )
XEQ "IRK" yields
x = 1
( in 5mn36s )
X<>Y
y(1) = 0.3678794411 ( error = -10-10
)
( With h = N = 1 the error is only -8.4 10-9
)
-The HP-41 displays | k1 - k'1 | + ............
+ | kn - k'n | where ki
and k'i are 2 successive k-approximations:
-At each step, this sum must tend to zero. Otherwise, reduce the stepsize
h and start over!
-Press R/S to continue with the same h and N.
-Start over with a smaller stepsize to obtain an error estimate.
b) Systems of 2 first
order differential equations: dy/dx = f (x;y;z) ;
dz/dx = g(x;y;z)
-The same formulae may be generalized to a system of 2 differential
equations, and we have to solve a 2nx2n non-linear system within each step.
Data Registers: ( Registers R00 to R06 and R2n+14 to Rn2+4n+13 are to be initialized before executing "IRKB" )
• R00: f name • R01
= x0
R07 to R13: temp
• R2n+14 to Rn2+4n+13 = the ( n2 + 2n
) coefficients ai,j ; bi ; ci
• R02 = y0
R14 to R2n+13 = ki
which are to be stored as shown below:
• R03 = z0
• R04 = h = stepsize
• R05 = N = number of steps
• R06 = n = number of stages
• R2n+14 = a1,1
• R2n+15 = a1,2 ...........................
• R3n+13 = a1,n
• R3n+14 = b1
• R3n+15 = a2,1
• R3n+16 = a2,2 ...........................
• R4n+14 = a2,n
• R4n+15 = b2
.....................................................................................................................................................
• Rn2+2n+13 = an,1 • Rn2+2n+14 = an,2 ...................... • Rn2+3n+12 = an,n • Rn2+3n+13 = bn
• Rn2+3n+14 = c1 • Rn2+3n+15 = c2 ........................ • Rn2+4n+13 = cn
Flags: /
Subroutine: a program calculating 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.
Z = z
>>>>> In other words,
this subroutine has to change the stack from Y = y
into Y' = f(x;y;z)
X = x
X' = g(x;y;z)
01 LBL "IRKB"
02 RCL 05
03 STO 07
04 RCL 06
05 ST+ X
06 13
07 +
08 .1
09 %
10 14
11 +
12 STO 09
13 CLRGX
if you don't have an HP-41CX , replace line 13 by the 5 lines:
0 LBL 00 STO IND Y
ISG Y GTO 00
14 LBL 01
15 RCL 09
16 FRC
17 14.9
18 STO 10
19 INT
20 +
21 STO 08
22 RCL 06
23 ST+ 10
24 ST+ 10
25 +
26 STO 09
27 CLX
28 STO 13
29 LBL 02
30 RCL 09
31 FRC
32 14
33 +
34 STO 11
35 RCL 06
36 +
37 STO 12
38 CLST
39 LBL 03
40 X<>Y
41 RCL IND 10
42 RCL IND 11
43 *
44 +
45 X<>Y
46 RCL IND 10
47 RCL IND 12
48 *
49 +
50 ISG 10
51 ISG 11
52 ISG 12
53 GTO 03
54 RCL 03
55 +
56 X<>Y
57 RCL 02
58 +
59 RCL 04
60 RCL IND 10
61 *
62 RCL 01
63 +
64 XEQ IND 00
65 RCL 04
66 ST* Z
67 *
68 ENTER^
69 X<> IND 09
70 -
71 ABS
72 X<>Y
73 ENTER^
74 X<> IND 08
75 -
76 ABS
77 +
78 ST+ 13
79 ISG 10
80 ISG 08
81 ISG 09
82 GTO 02
83 RCL 13
84 VIEW X
85 E-9
or another "small" number
86 X<Y?
87 GTO 01
88 RCL 06
89 ST- 11
90 ST- 12
91 CLST
92 LBL 04
93 X<>Y
94 RCL IND 10
95 RCL IND 11
96 *
97 +
98 X<>Y
99 RCL IND 10
100 RCL IND 12
101 *
102 +
103 ISG 10
104 ISG 11
105 ISG 12
106 GTO 04
107 ST+ 03
108 X<>Y
109 ST+ 02
110 RCL 04
111 ST+ 01
112 DSE 07
113 GTO 01
a three-byte GTO
114 RCL 03
115 RCL 02
116 RCL 01
117 CLD
118 END
( 177 bytes / SIZE n2+4n+14 )
| STACK | INPUTS | OUTPUTS |
| Z | / | z(x0+N.h) |
| Y | / | y(x0+N.h) |
| X | / | x0+N.h |
-If, for instance, you are using the 6-stage 12-order method described
above, the same 48 coefficients are to be stored into registers R26 thru
R73
Example: y(0) = 1 ; z(0)
= 0 ; dy/dx = z , dz/dx = -2 x.z -2 y ;
calculate y(1) and z(1)
01 LBL "U"
02 RCL Z
03 *
04 +
05 ST+ X
06 CHS
07 RTN
-Initialize registers R26 thru R73.
"U" ASTO 00
0 STO 01
1 STO 02
0 STO 03 and if we choose
h = 0.5
0.5 STO 04
2 STO 05 ( the number of
steps )
6 STO 06 ( the number of
stages )
XEQ "IRKB" yields
x = 1
execution time = 10mn43s
RDN
y(1) = 0.3678794411 ( error
= -10-10 )
RDN
z(0.5) = -0.7357588823 ( error
= 0 )
-The HP-41 displays | k1 - k'1 | + ............
+ | k2n - k'2n | where ki
and k'i are 2 successive k-approximations:
-At each step, this sum must tend to zero. Otherwise, reduce the stepsize
h and start over!
-Press R/S to continue with the same h and N.
-Start over with a smaller stepsize to obtain an error estimate.
3°) A slighly less general Implicit-Runge-Kutta
Program
a) First order differential
equations: dy/dx = f (x;y)
-Implicit Runge-Kutta methods are usually more "stable" than any explicit
method, but n-stage (2n-1) or (2n-2) order methods are even more stable
than n-stage 2n-order methods. ( For a detailed description
of stability properties of Runge-Kutta methods, cf Butcher's work )
-That's precisely what happens with the Radau IIA methods ( n-stage
(2n-1)-order methods ) and Lobatto IIIC methods ( n-stage (2n-2)-order
methods )
which also have another feature in common: ci
= an,i ( i = 1 , 2 , .... , n ) , thus, less data registers
are needed.
-The 2 following programs take advantage of this fact but the 2 programs
presented above would work too.
Note: The programs listed in paragraph 1 use
a 5-stage 8th-order Lobatto IIIC method.
Data Registers: ( Registers R00 to R05 and Rn+12 to Rn2+2n+11 are to be initialized before executing "I2RK" )
• R00: f name • R01
= x0
R06 to R11: temp
• Rn+12 to Rn2+2n+11 = the ( n2 + n )
coefficients ai,j ; bi
• R02 = y0
R12 = k1
which are to be stored as shown below:
• R03 = h = stepsize
R13 = k2
• R04 = N = number of steps
............
• R05 = n = number of stages
Rn+11 = kn
• Rn+12 = a1,1
• Rn+13 = a1,2 ...........................
• R2n+11 = a1,n •
R2n+12 = b1
• R2n+13 = a2,1
• R2n+14 = a2,2 ...........................
• R3n+12 = a2,n •
R3n+13 = b2
..................................................................................................................................................
• Rn2+n+11 = an,1 •
Rn2+n+12 = an,2 ........................
• Rn2+2n+10 = an,n • Rn2+2n+11
= bn
Flags: /
Subroutine: a program calculating f(x;y)
assuming x is in X-register and y is in Y-register upon entry.
01 LBL "I2RK"
02 RCL 04
03 STO 06
04 RCL 05
05 11
06 +
07 .1
08 %
09 12
10 +
11 STO 07
12 CLRGX
if you don't have an HP-41CX , replace line 12 by the 5 lines:
0 LBL 00 STO IND Y
ISG Y GTO 00
13 LBL 01
14 RCL 05
15 STO 08
16 RCL 07
17 FRC
18 12
19 ST+ 08
20 +
21 STO 07
22 CLX
23 ST0 10
24 LBL 02
25 RCL 07
26 FRC
27 12
28 +
29 STO 09
30 CLX
31 LBL 03
32 RCL IND 08
33 RCL IND 09
34 *
35 +
36 ISG 08
37 CLX
38 ISG 09
39 GTO 03
40 RCL 02
41 +
42 STO 11
43 RCL 03
44 RCL IND 08
45 *
46 RCL 01
47 +
48 XEQ IND 00
49 RCL 03
50 *
51 ENTER^
52 X<> IND 07
53 -
54 ABS
55 ST+ 10
56 ISG 08
57 CLX
58 ISG 07
59 GTO 02
60 RCL 10
61 VIEW X
62 E-9
or another "small" number
63 X<Y?
64 GTO 01
65 RCL 03
66 ST+ 01
67 RCL 11
68 STO 02
69 DSE 06
70 GTO 01
71 RCL 01
72 CLD
73 END
( 109 bytes / SIZE n2+2n+12 )
| STACK | INPUTS | OUTPUTS |
| Y | / | y(x0+N.h) |
| X | / | x0+N.h |
-For instance, here are the 56 coefficients of a 7-stage 13-order method
( a Radau IIA method )
-These coefficients are given with 20 decimals in case you want to
use them with a more accurate calculator:
R19 = a1,1
= 0.03754 62649 93921 33133
R27 = a2,1 = 0.08014 75965 15618 96780
R35 = a3,1 = 0.07206 38469 41881 90211
R20 = a1,2
= -0.01403 93345 56460 40154
R28 = a2,2 = 0.08106 20639 85891 53668
R36 = a3,2 = 0.17106 83549 83886 61942
R21 = a1,3
= 0.01035 27896 00742 30094
R29 = a2,3 = -0.02123 79921 20711 03494
R37 = a3,3 = 0.10961 45640 40072 10923
R22 = a1,4
= -0.00815 83225 40275 01191
R30 = a2,4 = 0.01400 02912 38817 11898
R38 = a3,4 = -0.02461 98717 28984 05386
R23 = a1,5
= 0.00638 84138 79534 68494
R31 = a2,5 = -0.01023 41857 30090 16383
R39 = a3,5 = 0.01476 03770 43950 81707
R24 = a1,6
= -0.00460 23267 79148 65550
R32 = a2,6 = 0.00715 34651 51364 59050
R40 = a3,6 = -0.00957 52593 96791 40056
R25 = a1,7
= 0.00182 89425 61470 64370
R33 = a2,7 = -0.00281 26393 72406 72334
R41 = a3,7 = 0.00367 26783 97138 30567
R26 = b1
= 0.02931 64271 59784 89197
R34 = b2 = 0.14807 85996 68484 29185
R42 = b3 = 0.33698 46902 81154
29910
R43 = a4,1
= 0.07570 51258 19824 42042
R51 = a5,1 = 0.07391 23421 63191 84654
R59 = a6,1 = 0.07470 55620 59796 23017
R44 = a4,2
= 0.15409 01551 42171 14465
R52 = a5,2 = 0.16135 56076 15942 43219
R60 = a6,2 = 0.15830 72238 72468 70066
R45 = a4,3
= 0.22710 77366 73202 38641
R53 = a5,3 = 0.20686 72415 52104 19782
R61 = a6,3 = 0.21415 34232 67200 03111
R46 = a4,4
= 0.11747 81870 37024 78199
R54 = a5,4 = 0.23700 71153 42694 23476
R62 = a6,4 = 0.21987 78470 31860 03999
R47 = a4,5
= -0.02381 08271 53044 17358
R55 = a5,5 = 0.10308 67935 33813 44662
R63 = a6,5 = 0.19875 21216 80635 26980
R48 = a4,6
= 0.01270 99855 33661 20563
R56 = a5,6 = -0.01885 41391 52580 44884
R64 = a6,6 = 0.06926 55016 05509 13323
R49 = a4,7
= -0.00460 88442 81289 63344
R57 = a5,7 = 0.00585 89009 74888 79182
R65 = a6,7 = -0.00811 60081 97728 29011
R50 = b4
= 0.55867 15187 71550 13208
R58 = b5 = 0.76923 38620 30054 50092
R66 = b6 = 0.92694 56713 19741 11485
R67 = a7,1
= 0.07449 42355 56010 31793
R68 = a7,2
= 0.15910 21157 33650 74087
R69 = a7,3
= 0.21235 18895 02977 80420
R70 = a7,4
= 0.22355 49145 07283 23475
( and ci = a7,i
)
R71 = a7,5
= 0.19047 49368 22115 57690
R72 = a7,6
= 0.11961 37446 12656 20289
R73 = a7,7
= 0.02040 81632 65306 12245 = 1/49
R74 = b7
= 1
Example: y(0) = 1 and dy/dx
= -2 x.y Evaluate y(1)
01 LBL "T"
02 *
03 ST+ X
04 CHS
05 RTN
-Initialize registers R19 thru R74
"T" ASTO 00
0 STO 01
1 STO 02 and if we choose
h = 0.5
0.5 STO 03
2 STO 04 ( the number of
steps )
7 STO 05 ( we are using a
7-stage method )
XEQ "I2RK" yields
x = 1
( in 7mn13s )
X<>Y
y(1) = 0.3678794410 ( error = -2 10-10
)
( With h = N = 1 the error is only -1.5 10-9 )
-The HP-41 displays | k1 - k'1 | + ............
+ | kn - k'n | where ki
and k'i are 2 successive k-approximations:
-At each step, this sum must tend to zero. Otherwise, reduce the stepsize
h and start over!
-Press R/S to continue with the same h and N.
-Start over with a smaller stepsize to obtain an error estimate.
b) Systems of 2 first
order differential equations: dy/dx = f (x;y;z) ;
dz/dx = g(x;y;z)
Data Registers: ( Registers R00 to R06 and R2n+14 to Rn2+4n+13 are to be initialized before executing "I2RKB" )
• R00: f name • R01
= x0
R07 to R15: temp
• R2n+16 to Rn2+3n+15 = the ( n2 + n
) coefficients ai,j ; bi
• R02 = y0
R16 to R2n+15 = ki
which are to be stored as shown below:
• R03 = z0
• R04 = h = stepsize
• R05 = N = number of steps
• R06 = n = number of stages
• R2n+16 = a1,1
• R2n+17 = a1,2 ...........................
• R3n+15 = a1,n
• R3n+16 = b1
• R3n+17 = a2,1
• R3n+18 = a2,2 ...........................
• R4n+16 = a2,n
• R4n+17 = b2
.....................................................................................................................................................
• Rn2+2n+15 = an,1 • Rn2+2n+16
= an,2 ...................... •
Rn2+3n+14 = an,n • Rn2+3n+15
= bn
Flags: /
Subroutine: a program calculating 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 "I2RKB"
02 RCL 05
03 STO 07
04 RCL 06
05 ST+ X
06 15
07 +
08 .1
09 %
10 16
11 +
12 STO 09
13 CLRGX
if you don't have an HP-41CX , replace line 13 by the 5 lines:
0 LBL 00 STO IND Y
ISG Y GTO 00
14 LBL 01
15 RCL 09
16 FRC
17 16
18 STO 10
19 +
20 STO 08
21 RCL 06
22 ST+ 10
23 ST+ 10
24 +
25 STO 09
26 CLX
27 STO 13
28 LBL 02
29 RCL 09
30 FRC
31 16
32 +
33 STO 11
34 RCL 06
35 +
36 STO 12
37 CLST
38 LBL 03
39 X<>Y
40 RCL IND 10
41 RCL IND 11
42 *
43 +
44 X<>Y
45 RCL IND 10
46 RCL IND 12
47 *
48 +
49 ISG 10
50 CLX
51 ISG 11
52 ISG 12
53 GTO 03
54 RCL 03
55 +
56 STO 15
57 X<>Y
58 RCL 02
59 +
60 STO 14
61 RCL 04
62 RCL IND 10
63 *
64 RCL 01
65 +
66 XEQ IND 00
67 RCL 04
68 ST* Z
69 *
70 ENTER^
71 X<> IND 09
72 -
73 ABS
74 X<>Y
75 ENTER^
76 X<> IND 08
77 -
78 ABS
79 +
80 ST+ 13
81 ISG 10
82 CLX
83 ISG 08
84 ISG 09
85 GTO 02
86 RCL 13
87 VIEW X
88 E-9
or another "small" number
89 X<Y?
90 GTO 01
91 RCL 04
92 ST+ 01
93 RCL 15
94 STO 03
95 RCL 14
96 STO 02
97 DSE 07
98 GTO 01
a three-byte GTO
99 RCL 01
100 CLD
101 END
( 146 bytes / SIZE n2+3n+16 )
| STACK | INPUTS | OUTPUTS |
| Z | / | z(x0+N.h) |
| Y | / | y(x0+N.h) |
| X | / | x0+N.h |
-If, for instance, you are using the 7-stage 13-order method described
above, the same 56 coefficients are to be stored into registers R30 thru
R85.
Example: y(0) = 1 ; z(0)
= 0 ; dy/dx = z , dz/dx = -2 x.z -2 y ;
calculate y(1) and z(1)
01 LBL "U"
02 RCL Z
03 *
04 +
05 ST+ X
06 CHS
07 RTN
-Initialize registers R30 thru R85.
"U" ASTO 00
0 STO 01
1 STO 02
0 STO 03 and if we choose
h = 0.5
0.5 STO 04
2 STO 05 ( the number of
steps )
7 STO 06 ( the number of
stages )
XEQ "I2RKB" yields
x = 1
execution time = 13mn54s
RDN
y(1) = 0.3678794411 ( error
= -10-10 )
RDN
z(0.5) = -0.7357588822 ( error
= 10-10 )
-The HP-41 displays | k1 - k'1 | + ............
+ | k2n - k'2n | where ki
and k'i are 2 successive k-approximations:
-At each step, this sum must tend to zero. Otherwise, reduce the stepsize
h and start over!
-Press R/S to continue with the same h and N.
-Start over with a smaller stepsize to obtain an error estimate.
Reference:
J.C. Butcher - "The Numerical Analysis of Ordinary Differential Equations"
- John Wiley & Sons - ISBN 0-471-91046-5
Go back to the HP-41 software library
Go back to the general software library
Go
back to the main exhibit hall