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)
with the initial values: 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: • R00: f name ( Registers R00 to R04 and R16 to R45 are to be initialized before executing "IRK8" )
• 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 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: • R00: f name ( Registers R00 to R05 and R25 to R54 are to be initialized before executing "IRK8B" )
• 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.
-If you don't have an HP-41CX,
replace lines 04-05 with CLX STO 15 STO 16
STO 17 STO 18 STO 19 STO 20 STO 21 STO 22
STO 23 STO 24
-Line 80 is a three-byte GTO 01
01 LBL "IRK8B"
02 RCL 05 03 STO 12 04 15.024 05 CLRGX 06 LBL 01 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 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 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: • R00: f name ( Registers R00 to R05 and Rn+11 to Rn2+3n+10 are to be initialized before executing "IRK" )
• 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.
-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
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 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 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: • R00: f name ( Registers R00 to R06 and R2n+14 to Rn2+4n+13 are to be initialized before executing "IRKB" )
• 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.
-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
-Line 113 is a three-byte GTO 01.
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 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 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 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: • R00: f name ( Registers R00 to R05 and Rn+12 to Rn2+2n+11 are to be initialized before executing "I2RK" )
• 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.
-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
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 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 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: • R00: f name ( Registers R00 to R06 and R2n+14 to Rn2+4n+13 are to be initialized before executing "I2RKB" )
• 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.
-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
-Line 98 is a three-byte GTO 01
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 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 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 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:
[1] J.C. Butcher - "The Numerical Analysis of Ordinary Differential
Equations" - John Wiley & Sons - ISBN 0-471-91046-5