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
-Two kinds of programs are presented below:
1°) IRK is a general Runge-Kutta program suitable
to all implicit formulae
2°) I2RK is slightly less general.
-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)
------------------- y(x0) = y0
; z(x0) = z0
1°) General implicit Runge-Kutta methods
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)
001 LBL "IRKB"
002 RCL 05
003 STO 07
004 RCL 06
005 ST+ X
006 13
007 +
008 .1
009 %
010 14
011 +
012 STO 09
013 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
014 LBL 01
015 RCL 09
016 FRC
017 14.9
018 STO 10
019 INT
020 +
021 STO 08
022 RCL 06
023 ST+ 10
024 ST+ 10
025 +
026 STO 09
027 CLX
028 STO 13
029 LBL 02
030 RCL 09
031 FRC
032 14
033 +
034 STO 11
035 RCL 06
036 +
037 STO 12
038 CLST
039 LBL 03
040 X<>Y
041 RCL IND 10
042 RCL IND 11
043 *
044 +
045 X<>Y
046 RCL IND 10
047 RCL IND 12
048 *
049 +
050 ISG 10
051 ISG 11
052 ISG 12
053 GTO 03
054 RCL 03
055 +
056 X<>Y
057 RCL 02
058 +
059 RCL 04
060 RCL IND 10
061 *
062 RCL 01
063 +
064 XEQ IND 00
065 RCL 04
066 ST* Z
067 *
068 ENTER^
069 X<> IND 09
070 -
071 ABS
072 X<>Y
073 ENTER^
074 X<> IND 08
075 -
076 ABS
077 +
078 ST+ 13
079 ISG 10
080 ISG 08
081 ISG 09
082 GTO 02
083 RCL 13
084 VIEW X
085 E-9
or another "small" number
086 X<Y?
087 GTO 01
088 RCL 06
089 ST- 11
090 ST- 12
091 CLST
092 LBL 04
093 X<>Y
094 RCL IND 10
095 RCL IND 11
096 *
097 +
098 X<>Y
099 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.
2°) Implicit Runge-Kutta methods, a special
case:
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 last programs listed in "Runge-Kutta methods
for the HP-41" use a 5-stage 8-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.
001 LBL "I2RKB"
002 RCL 05
003 STO 07
004 RCL 06
005 ST+ X
006 15
007 +
008 .1
009 %
010 16
011 +
012 STO 09
013 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
014 LBL 01
015 RCL 09
016 FRC
017 16
018 STO 10
019 +
020 STO 08
021 RCL 06
022 ST+ 10
023 ST+ 10
024 +
025 STO 09
026 CLX
027 STO 13
028 LBL 02
029 RCL 09
030 FRC
031 16
032 +
033 STO 11
034 RCL 06
035 +
036 STO 12
037 CLST
038 LBL 03
039 X<>Y
040 RCL IND 10
041 RCL IND 11
042 *
043 +
044 X<>Y
045 RCL IND 10
046 RCL IND 12
047 *
048 +
049 ISG 10
050 CLX
051 ISG 11
052 ISG 12
053 GTO 03
054 RCL 03
055 +
056 STO 15
057 X<>Y
058 RCL 02
059 +
060 STO 14
061 RCL 04
062 RCL IND 10
063 *
064 RCL 01
065 +
066 XEQ IND 00
067 RCL 04
068 ST* Z
069 *
070 ENTER^
071 X<> IND 09
072 -
073 ABS
074 X<>Y
075 ENTER^
076 X<> IND 08
077 -
078 ABS
079 +
080 ST+ 13
081 ISG 10
082 CLX
083 ISG 08
084 ISG 09
085 GTO 02
086 RCL 13
087 VIEW X
088 E-9
or another "small" number
089 X<Y?
090 GTO 01
091 RCL 04
092 ST+ 01
093 RCL 15
094 STO 03
095 RCL 14
096 STO 02
097 DSE 07
098 GTO 01
a three-byte GTO
099 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