Overview
| XROM | Function | Desciption |
| 15,00
15,01 15,02 15,03 15,04 15,05 15,06 15,07 15,08 15,09 15,10 15,11 15,12 15,13 15,14 15,15 15,16 15,17 15,18 15,19 15,20 15,21 15,22 15,23 15,24 15,25 15,26 15,27 15,28 15,29 |
-DIFFEQTN
"C41" "C51" "RK4" "RK4B" "RK4C" "RK4N" "RK6" "RK6B" "RK6N" "RK8" "RK8B" "RK8N" "RKF" "RKFB" "RKFN" "RKS2" "RKS3" "RKSN" "RKV" "RKVB" -BULIRSHSTR "BST" "BST2" "BSTN" "STRM" -NTH-ORDER "NRK4" "SRK4" "TRK4" |
Section Header
Constants for "RK8" Constants for "RKV" Runge-Kutta 4th order formulae RK4 for Systems of 2 Differential Equations RK4 for Systems of 3 Differential Equations RK4 for Systems of n Differential Equations Runge-Kutta 6th order formulae RK6 for Systems of 2 Differential Equations RK6 for Systems of n Differential Equations Runge-Kutta 8th order formulae RK8 for Systems of 2 Differential Equations RK8 for Systems of n Differential Equations Runge-Kutta-Felhberg 4-5 method RKF for Systems of 2 Differential Equations RKF for Systems of n Differential Equations RK4 for Systems of 2 Second order Conservative Eq RK4 for Systems of 3 Second order Conservative Eq RK4 for Systems of n Second order Conservative Eq Runge-Kutta-Verner 5-6 method RKV for Systems of 2 Differential Equations Section Header Bulirsh-Stoer Method BST for Systems of 2 Differential Equations BST for Systems of n Differential Equations Stoermer's rule for Conservative 2nd-order Equation Section Header N-th order differential equations ( 4th order method ) 2-nd order differential equations ( 4th order method ) 3-rd order differential equations ( 4th order method ) |
"C41" stores these 41 constants into registers R11 thru R51
They are used by an 11-stage method of orde 8 ( RK8 )
| R11 = 0.8273268354
R12 = 0.5 R13 = 0.1726731646 R14 = 9 R15= 14 R16 = 0.25 R17 = 0.8961811934 R18 = -0.2117115009 R19 = 7 R20 = 0.06514850915 R21 = 0.5766714727 |
R22 = 0.1855068535
R23 = 0.3865246267 R24 = -0.4634553896 R25 = 0.3772937693 R26 = 0.1996369936 R27 = 0.09790045616 R28 = 0.3285172131 R29 = -0.3497052863 R30 = -0.03302551131 R31 = 0.1289862930 R32 = -0.01186868389 |
R33 = 0.002002165993
R34 = 0.9576053440 R35 = -0.6325461607 R36 = 0.1527777778 R37 = -0.009086961101 R38 = 0.03125 R39 = 1.062498447 R40 = -1.810863083 R41 = 2.031083139 R42 = -0.6379313502 R43 = 0.9401094520 |
R44 = -2.229158210
R45 = 7.553840442 R46 = -7.164951553 R47 = 2.451380432 R48 = -0.5512205631 R49 = 0.05 R50 = 0.2722222222 R51 = 2.8125 |
"C51" stores 51 constants into R19 thru R69
These constants are used by "RKV" & "RKVB" which employs
an 8-stage method with order 5 with a 6th-order error-estimating formula.
| R19 = 1/18
R21 = -1/12 R24 = -2/81 R28 = 40/33 R33 = -369/73 R39 = -8716/891 R46 = 3015/256 |
R22 = 1/4 R25 = 4/27 R29 = -4/11 R34 = 72/73 R40 = 656/297 R47 = -9/4 |
R26 = 8/81 R30 = -56/11 R35 = 5380/219 R41 = 39520/891 R48 = -4219/78 |
R31 = 54/11 R36 = -12285/584 R42 = -416/11 R49 = 5985/128 |
R37 = 2695/1752 R43 = 52/27 R50 = -539/384 |
R44 = 0 R51 = 0 |
R52 = 693/3328 |
R20 = 1/18
R23 = 1/6 R27 = 2/9 R32 = 2/3 R38 = 1 R45 = 8/9 R53 = 1 |
| R54 = 3/80
R62 = 33/640 |
R55 = 0
R63 = 0 |
R56 = 4/25
R64 = -132/325 |
R57 = 243/1120
R65 = 891/2240 |
R58 = 77/160
R66 = -33/320 |
R59 = 73/700
R67 = -73/700 |
R60 = 0
R68 = 891/8320 |
R61 = 0
R69 = 2/35 |
Runge-Kutta 4th order formulae
-We solve dy/dx = f(x,y) with the initial value: y(x0) = y0
RK4 uses the following formulae:
k1 = h.f(x;y)
k2 = h.f(x+h/2;y+k1/2)
k3 = h.f(x+h/2;y+k2/2)
k4 = h.f(x+h;y+k3)
and then, y(x+h) = y(x) + ( k1 + 2.k2
+ 2.k3 + k4 ) / 6 This
formula duplicates the Taylor series up to the term in h4
Data Registers: • R00: f name ( Registers R00 to R04 are to be initialized before executing "RK4" )
• R01 = x0 •
R03 = h = stepsize
R05 & R07: temp
• R02 = y0 •
R04 = N = number of steps
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)
| STACK | INPUTS | OUTPUTS |
| Y | / | y(x0+N.h) |
| X | / | x0+N.h |
Example: y(0) = 1 and dy/dx
= -2 x.y Evaluate y(1)
( The solution is y = exp ( -x2 ) )
| 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
10 STO 04 XEQ "RK4"
yields x = 1
( in 20s )
X<>Y
y(1) = 0.367881
( the exact result is 1/e = 0.367879441 )
-Press R/S to continue with the same h and N.
-Start over with a smaller stepsize to obtain an error estimate ( when
h is divided by 2 , errors are approximately divided by 24 =
16 )
RK4B: RK4 for Systems of 2 Differential Equations
-"RK4B" solves dy/dx = f(x,y,z) , dz/dx = g(x,y,z) with the initial values: y(x0) = y0 , z(x0) = z0
-The preceding formulae are now generalized to a system of 2 differential equations:
k1 = h.f(x;y;z)
l1 = h.g(x;y;z)
k2 = h.f(x+h/2;y+k1/2;z+l1/2)
l2 = h.g(x+h/2;y+k1/2;z+l1/2)
k3 = h.f(x+h/2;y+k2/2;z+l2/2)
l3 = h.g(x+h/2;y+k2/2;z+l2/2)
k4 = h.f(x+h;y+k3;z+l3)
l4 = h.g(x+h;y+k3;z+l3)
then, y(x+h) = y(x) + ( k1 + 2.k2
+ 2.k3 + k4 ) / 6
and z(x+h) = z(x) + ( l1 + 2.l2
+ 2.l3 + l4 ) / 6
duplicate the Taylor series through the terms in h4
Data Registers: • R00: f name ( Registers R00 to R05 are to be initialized before executing "RK4B" )
• R01 = x0 •
R04 = h = stepsize
R06 to R09: temp
• R02 = y0 •
R05 = N = number of steps
• 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.
>>>>> In other words, this subroutine has to change the stack from
Z = z
Y = y to
Y' = f(x;y;z)
X = x
X' = g(x;y;z)
| STACK | INPUTS | OUTPUTS |
| Z | / | z(x0+N.h) |
| Y | / | y(x0+N.h) |
| X | / | x0+N.h |
Example: y(0) = 1 ; y'(0) = 0
; y'' + 2.x.y' + 2.y = 0 Evaluate y(1) and y'(1)
( The solution is y = exp ( -x2 ) ; y' = -2.x exp
( -x2 ) )
-This second order equation is equivalent to the system:
y(0) = 1 ; z(0) = 0 ; dy/dx = z , dz/dx = -2 x.z -2 y
where z = y'
| 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
10 STO 05 XEQ "RK4B"
yields x = 1
( in 32s )
RDN
y(1) = 0.367881
the exact result is y(1) = 1/e = 0.367879441
RDN
z(1) = -0.735762
the exact result is z(1) = -2/e = -0.735758883
-Press R/S to continue with the same h and N.
N.B:
-To obtain an error estimate, start over again with a smaller stepsize:
when h is divided by 2 , errors are approximately divided by
16
RK4C: RK4 for Systems of 3 Differential Equations
-"RK4C" solves dy/dx = f(x,y,z,u) , dz/dx
= g(x,y,z,u) , du/dx = h(x,y,z,u) with the initial
values: y(x0) = y0 ,
z(x0) = z0 , u(x0) = u0
Data Registers: • R00: f name ( Registers R00 to R06 are to be initialized before executing "RK4C" )
• R01 = x0 •
R05 = h = stepsize
R07 to R11: temp
• R02 = y0 •
R06 = N = number of steps
• R03 = z0
• R04 = u0
Flags: /
Subroutine: a program calculating f(x;y;z;u)
in Z-register , g(x;y;z;u) in Y-register and h(x;y;z;u) in X-register
assuming x is in X-register , y is in Y-register , z is in Z-register and
u is in T-register upon entry.
>>>>> In other words, this subroutine has to change the stack from
T = u
Z = z
Z' = f(x;y;z;u)
Y = y to
Y' = g(x;y;z;u)
X = x
X' = h(x;y;z;u)
| STACK | INPUTS | OUTPUTS |
| T | / | u(x0+N.h) |
| Z | / | z(x0+N.h) |
| Y | / | y(x0+N.h) |
| X | / | x0+N.h |
Example:
dy/dx = -y.z.u
y(0) = 1 Evaluate
y(1) , z(1) , u(1) with h = 0.1 ( whence
N = 10 )
dz/dx = x.( y + z - u )
z(0) = 1
du/dx = x.y - z.u
u(0) = 2
| LBL "T"
R^ CHS STO L R^ ST* Y ST+ L CLX RCL Z R^ ST+ L ST* Y X<> Z ST+ Y ST* Z X<> T LASTX * X<>Y RTN |
"T" ASTO 00
0 STO 01
1 STO 02 STO 03
2 STO 04
0.1 STO 05
10 STO 06 XEQ "RK4C"
>>>> x = 1
( in 53 seconds )
RDN y(1) = 0.258209
RDN z(1) = 1.157620
RDN u(1) = 0.842179
-Press R/S to continue with the same h and N.
RK4N: RK4 for Systems of n Differential Equations
-"RK4N" solves the following system with the "classical" 4th order Runge-Kutta method.
dy1/dx = f1(x,y1,
... ,yn)
dy2/dx = f2(x,y1,
... ,yn)
..................
with the initial values: yi(x0)
i = 1 , ... , n
dyn/dx = fn(x,y1, ... ,yn)
Data Registers: • R00 = subroutine name ( Registers R00 thru R03 and R10 thru R10+n are to be initialized before executing "RK4N" )
• R01 = n = number of equations = number of functions
R04 thru R09: temp
• R02 = h = stepsize
• R03 = N = number of steps
• R10 = x0
• R11 = y1(x0)
R11+n thru R10+4n: temp
• R12 = y2(x0)
............... ( initial
values )
• R10+n = yn(x0)
Flags: F06-F07
Subroutine: A program which calculates
and stores f1 ( x, y1 , ... , yn)
, ... , fn ( x, y1 , ... , yn) in
registers R11+n, ... , R10+2n respectively
with x, y1, ... , yn in regiters R10, R11,
... , R10+n
| STACK | INPUTS | OUTPUTS |
| T | / | y3(x0+N.h) |
| Z | / | y2(x0+N.h) |
| Y | / | y1(x0+N.h) |
| X | / | x0+N.h |
Example: Let's solve the system:
dy1/dx = - y1y2y3
dy2/dx = x ( y1+y2-y3
)
dy3/dx = xy1-
y2y3
with the initial values: y1(0) = 1 ; y2(0) = 1 ; y3(0) = 2
x ; y1 ; y2 ; y3 will be
in R10 ; R11 ; R12 ; R13 respectively, and we have to write a program to
compute the 3 functions
and store the results in R14 ; R15 ; R16 ( just after
the "initial-value registers" ).
-For instance, the following subroutine will do the job - let's call
it "T"
| 01 LBL "T"
02 RCL 11 03 RCL 12 04 RCL 13 05 * 06 * 07 CHS 08 STO 14 09 RCL 11 10 RCL 12 11 + 12 RCL 13 13 - 14 RCL 10 15 * 16 STO 15 17 RCL 10 18 RCL 11 19 * 20 RCL 12 21 RCL 13 22 * 23 - 24 STO 16 25 RTN |
We store the name of the subroutine in R00 "T"
ASTO 00
the numbers of equations in R01 3 STO 01
then the initial values:
x in R10
0 STO 10
y1 in R11
1 STO 11
y2 in R12
1 STO 12
y3 in R13
2 STO 13
-Suppose we want to know y1(1) , y2(1) , y3(1) with a step size h = 0.1 and therefore N = 10
h is stored in R02
0.1 STO 02
N is stored in R03
10 STO 03 and XEQ "RK4N"
about 2mn24s later, the HP-41 returns:
x = 1
in X and R10
y1(1) =
0.2582093855 in Y and R11
y2(1) =
1.157619553 in Z and R12
y3(1) =
0.8421786516 in T and R13
-To continue with the same h and N, simply press R/S and you'll have y1(2) y2(2) y3(2)
-An estimation of the accuracy may be obtained by doing the calculation again with a smaller step size like h=0.05
It yields: y1(1) = 0.2582079989 y2(1) = 1.157623732 y3(1) = 0.8421783424
-Theoretically, the results tend to the exact solution as h tends to
zero.
Runge-Kutta 6th order formulae
dy/dx = f(x,y) with the initial value: y(x0) = y0
-The formula duplicates the Taylor series up to the term in h6
-Methods of this order require 7 stages ( 7 evaluations of the function
f within each step)
Data Registers: • R00: f name ( Registers R00 to R04 are to be initialized before executing "RK6" )
• R01 = x0 •
R03 = h = the stepsize
R05: counter
• R02 = y0 •
R04 = N = the number of steps R06 , ......
, R11 = k1 , ..... , k6
Flags: /
Subroutine: a program calculating f(x;y)
assuming x is in X-register and y is in Y-register upon entry.
| STACK | INPUTS | OUTPUTS |
| Y | / | y(x0+N.h) |
| X | / | x0+N.h |
Example: y(0) = 1 and dy/dx
= -2 x.y Evaluate y(1)
| 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
10 STO 04 XEQ "RK6"
yields x = 1
( in 84s )
X<>Y
y(1) = 0.367879436 ( error = -5 10-9
)
-Press R/S to continue with the same h and N.
-Start over with a smaller stepsize to obtain an error estimate ( when
h is divided by 2 , errors are approximately divided by 26 =
64 )
RK6B: RK6 for Systems 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
Data Registers: • R00: f name ( Registers R00 to R05 are to be initialized before executing "RK6B" )
• R01 = x0 •
R04 = h = stepsize
R06 to R16: temp
• R02 = y0 •
R05 = N = number of steps
• 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.
| 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(1) and z(1)
| 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
10 STO 05 XEQ "RK6B"
yields x = 1
( in 2mn30s )
RDN
y(1) = 0.367879433
( y-error = -9 10-9 )
RDN
z(1) = -0.735758865
( z-error = 17 10-9 )
-Press R/S to continue with the same h and N.
-Start over with a smaller stepsize to obtain an error estimate ( when
h is divided by 2 , errors are approximately divided by 26 =
64 )
RK6N: RK6 for Systems of n Differential Equations
dy1/dx = f1(x,y1,
... ,yn)
dy2/dx = f2(x,y1,
... ,yn)
..................
with the initial values: yi(x0)
i = 1 , ... , n
dyn/dx = fn(x,y1,
... ,yn)
Data Registers: • R00 = subroutine name ( Registers R00 thru R03 and R20 thru R20+n are to be initialized before executing "RK6N" )
• R01 = n = number of equations = number of functions
R04 thru R13: temp R14 thru R19:
unused
• R02 = h = stepsize
• R03 = N = number of steps
• R20 = x0
• R21 = y1(x0)
R21+n thru R20+8n: temp
• R22 = y2(x0)
............... ( initial
values )
• R20+n = yn(x0)
Flags: /
Subroutine: A program which calculates
and stores f1 ( x, y1 , ... , yn)
, ... , fn ( x, y1 , ... , yn) in
registers R21+n, ... , R20+2n respectively
with x, y1, ... , yn in regiters R20, R21,
... , R20+n
| STACK | INPUTS | OUTPUTS |
| T | / | y3(x0+N.h) |
| Z | / | y2(x0+N.h) |
| Y | / | y1(x0+N.h) |
| X | / | x0+N.h |
Example:
dy/dx = -y.z.u
y(0) = 1 Evaluate
y(1) , z(1) , u(1) with h = 0.1 ( whence
N = 10 )
dz/dx = x.( y + z - u )
z(0) = 1
du/dx = x.y - z.u
u(0) = 2
| 01 LBL "T"
02 RCL 21 03 RCL 22 04 RCL 23 05 * 06 * 07 CHS 08 STO 24 09 RCL 21 10 RCL 22 11 + 12 RCL 23 13 - 14 RCL 20 15 * 16 STO 25 17 LASTX 18 RCL 21 19 * 20 RCL 22 21 RCL 23 22 * 23 - 24 STO 26 25 RTN |
"T" ASTO 00
Initial values:
0
STO 20
1
STO 21 STO 22
2
STO 23
n = 3 STO 01 (
3 functions )
h = 0.1 STO 02
N = 10 STO 03 XEQ "RK6N"
>>>> x = 1
= R20
( in 6mn28s )
RDN y(1) = 0.258207889 = R21
y(1) = 0.258207906
RDN z(1) = 1.157623948 = R22
the exact results, rounded to 9D are
z(1) = 1.157623981
RDN u(1) = 0.842178329 = R23
u(1) = 0.842178312
-To continue the calculations, simply press R/S ( after changing
h & N in registers R02 & R03 if needed )
Runge-Kutta 8th order formulae
-Methods of this order require 11 stages. The following one was discovered
in 1972 by Cooper and Verner.
-The formula duplicates the Taylor series up to the term in h8
-"RK8" uses 41 coefficients which are of the form ( a + b 211/2
)/c where a , b , c are integers.
-They are stored in registers R11 thru R51 by "C41"
Data Registers: • R00: f name ( Registers R00 to R04 are to be initialized before executing "RK8" )
• R01 = x0 •
R03 = h = the stepsize
R05 to R10 and R52: temp
• R02 = y0 •
R04 = N = the number of steps
R11 to R51: constants
Flags: /
Subroutine: a program calculating f(x;y)
assuming x is in X-register and y is in Y-register upon entry.
| STACK | INPUTS | OUTPUTS |
| Y | / | y(x0+N.h) |
| X | / | x0+N.h |
Example: y(0) = 1 and dy/dx
= -2 x.y Evaluate y(1)
| 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
10 STO 04 XEQ "RK8" yields
x = 1
( in 2mn )
X<>Y
y(1) = 0.3678794412 correct to 10 D!
-Press R/S to continue with the same h and N.
-Start over with a smaller stepsize to obtain an error estimate ( when
h is divided by 2 , errors are approximately divided by 28 =
256 )
RK8B: RK8 for Systems of 2 Differential Equations
Data Registers: • R00: f name ( Registers R00 to R05 are to be initialized before executing "RK8B" )
• R01 = x0 •
R04 = h = stepsize
R06 to R17: temp R59: counter
• R02 = y0 •
R05 = N = number of steps
R18 thru R58 = the same 41 numbers used in "RK8".
• 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.
| 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(1) and z(1)
| 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
10 STO 05 XEQ "RK8B" yields
x = 1
in 3mn10s
RDN
y(1) = 0.3678794412 correct to 10D.
RDN
z(1) = -0.7357588824 z-error = -10-10
-Press R/S to continue with the same h and N.
-Start over with a smaller stepsize to obtain an error estimate ( when
h is divided by 2 , errors are approximately divided by 28 =
256 )
RK8N: RK8 for Systems of n Differential Equations
dy1/dx = f1(x,y1,
... ,yn)
dy2/dx = f2(x,y1,
... ,yn)
..................
with the initial values: yi(x0)
i = 1 , ... , n
dyn/dx = fn(x,y1,
... ,yn)
Data Registers: • R00 = subroutine name ( Registers R00 to R03 & R50 to R50+n are to be initialized before executing "RK8N" )
• R01 = n = number of equations = number of functions
R04 thru R13: temp R49: temp
• R02 = h = stepsize
• R03 = N = number of steps
R14 thru R45 = Constants
R46-R47-R48 are unused
• R50 = x0
• R51 = y1(x0)
R51+n thru R50+9n: temp
• R52 = y2(x0)
............... ( initial
values )
• R50+n = yn(x0)
Flags: /
Subroutine: A program which calculates
and stores f1 ( x, y1 , ... , yn)
, ... , fn ( x, y1 , ... , yn) in
registers R51+n, ... , R50+2n respectively
with x, y1, ... , yn in regiters R50, R51,
... , R50+n
| STACK | INPUTS | OUTPUTS |
| T | / | y3(x0+N.h) |
| Z | / | y2(x0+N.h) |
| Y | / | y1(x0+N.h) |
| X | / | x0+N.h |
Example:
dy/dx = -y.z.u
y(0) = 1 Evaluate
y(1) , z(1) , u(1) with h = 0.1 and N =
10
dz/dx = x.( y + z - u )
z(0) = 1
du/dx = x.y - z.u
u(0) = 2
| 01 LBL "T"
02 RCL 51 03 RCL 52 04 RCL 53 05 * 06 * 07 CHS 08 STO 54 09 RCL 51 10 RCL 52 11 + 12 RCL 53 13 - 14 RCL 50 15 * 16 STO 55 17 LASTX 18 RCL 51 19 * 20 RCL 52 21 RCL 53 22 * 23 - 24 STO 56 25 RTN |
"T" ASTO 00
Initial values:
0 STO 50
1 STO 51
STO 52
2 STO 53
n = 3 STO 01 (
3 functions )
h = 0.1 STO 02
N = 10 STO 03 XEQ "RK8N"
>>>> x = 1
= R50
( in 9mn31s )
RDN y(1) = 0.258207906 = R51
y(1) = 0.258207906
RDN z(1) = 1.157623979 = R52
the exact results, rounded to 9D are
z(1) = 1.157623981
RDN u(1) = 0.842178313 = R53
u(1) = 0.842178312
-To continue the calculations, simply press R/S ( after changing
h & N in registers R02 & R03 if needed )
Runge-Kutta-Felhberg 4-5 method
-In this Runge-Kutta-Fehlberg method, a 4th-order formula is used to
compute the solution
and the difference between this 4th-order formula and a 5th-order
formula provides an error estimate.
Data Registers: • R00: f name ( Registers R00 to R04 are to be initialized before executing "RKF" )
• R01 = x0 •
R03 = h = stepsize
R05 = error estimate
• R02 = y0 •
R04 = N = number of steps
R06 thru R10: temp
Flags: /
Subroutine: a program calculating f(x;y)
assuming x is in X-register and y is in Y-register upon entry.
| STACK | INPUTS | OUTPUTS |
| Y | / | y(x0+N.h) |
| X | / | x0+N.h |
Example: y(0) = 1 and dy/dx
= -2 x.y Evaluate y(1)
| 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
10 STO 04 XEQ "RKF" yields
x = 1
X<>Y
y(1) = 0.367879263
and the error estimate RCL 05 >>> -9.7 10-8
-The error is actually -1.8 10-7
-The successive errors are simply added to R05, so they are sometimes
underestimated...
-Press R/S to continue with the same h and N.
RKFB: RKF for Systems of 2 Differential Equations
Data Registers: • R00: f name ( Registers R00 to R05 are to be initialized before executing "RKFB" )
• R01 = x0 •
R04 = h = stepsize
R06 = y-error estimate
R08 to R16: temp
• R02 = y0 •
R05 = N = number of steps
R07 = z-error estimate
• 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.
| 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(1) and z(1)
| 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
10 STO 05 XEQ "RKFB" yields
x = 1
( in 2mn16s )
RDN
y(1) = 0.367879517
RDN
z(1) = -0.735759034
and the error RCL 06 >>>
-8.7 10-8
actually y-error = 7.6 10-8
estimates: RCL 07
>>> -2.1 10-7
actually z-error = -1.5 10-7
-Press R/S to continue with the same h and N.
RKFN: RKF for Systems of n Differential Equations
dy1/dx = f1(x,y1,
... ,yn)
dy2/dx = f2(x,y1,
... ,yn)
..................
with the initial values: yi(x0)
i = 1 , ... , n
dyn/dx = fn(x,y1,
... ,yn)
Data Registers: • R00 = subroutine name ( Registers R00 thru R03 and R20 thru R20+n are to be initialized before executing "RKFN" )
• R01 = n = number of equations = number of functions
R04 thru R14: temp R15 thru R19:
unused
• R02 = h = stepsize
• R03 = N = number of steps
• R20 = x0
at the end:
• R21 = y1(x0)
R21+n = err(y1)
R21+n thru R20+8n: temp
• R22 = y2(x0)
R22+n = err(y2)
...............
......................
• R20+n = yn(x0) R20+2n = err(yn)
Where err(yi) is an error-estimate of yi ( during the calculations, these error-estimates are moved in registers R21+2n .... R20+3n )
Flags: /
Subroutine: A program which calculates
and stores f1 ( x, y1 , ... , yn)
, ... , fn ( x, y1 , ... , yn) in
registers R21+n, ... , R20+2n respectively
with x, y1, ... , yn in regiters R20, R21,
... , R20+n
| STACK | INPUTS | OUTPUTS |
| T | / | y3(x0+N.h) |
| Z | / | y2(x0+N.h) |
| Y | / | y1(x0+N.h) |
| X | / | x0+N.h |
Example:
dy/dx = -y.z.u
y(0) = 1 Evaluate
y(1) , z(1) , u(1) with h = 0.1 & N
= 10
dz/dx = x.( y + z - u )
z(0) = 1
du/dx = x.y - z.u
u(0) = 2
| 01 LBL "T"
02 RCL 21 03 RCL 22 04 RCL 23 05 * 06 * 07 CHS 08 STO 24 09 RCL 21 10 RCL 22 11 + 12 RCL 23 13 - 14 RCL 20 15 * 16 STO 25 17 LASTX 18 RCL 21 19 * 20 RCL 22 21 RCL 23 22 * 23 - 24 STO 26 25 RTN |
"T" ASTO 00
Initial values:
0
STO 20
1
STO 21 STO 22
2
STO 23
n = 3 STO 01 (
3 functions )
h = 0.1 STO 02
N = 10 STO 03 XEQ "RKF" >>>>
x = 1
= R20
( in 5mn40s )
RDN y(1) = 0.258207319 = R21
err(y) = -603 E-9 = R24
RDN z(1) = 1.157624973 = R22
and the error estimates
err(z) = 1062 E-9 = R25
RDN u(1) = 0.842178529 = R23
err(u) = 1768 E-9 = R26
-Actually, the true errors are -587 E-9 for y(1) 992 E-9 for z(1) 217 E-9 for u(1)
-To continue the calculations, simply press R/S ( after changing
h & N in registers R02 & R03 if needed )
RKS2: RK4 for Systems of 2 Second order Conservative
Equations
-"RKS2" solves the system:
y" = f(x,y,z)
y(x0) = y0 y'(x0)
= y'0
z" = g(x,y,z)
z(x0) = z0 z'(x0)
= z'0
by the 4th order formulae:
y(x+h) = y(x) + h ( y'(x) + k1/6 + 2k2/3
)
y'(x+h) = y'(x) + k1/6 + 2k2/3 +
k3/6
where k1 = h.f (x,y) ; k2 = h.f(x+h/2,y+h.y'/2+h.k1/8)
; k3 = h.f(x+h,y+h.y'+h.k2/2)
Data Registers: • R00: subroutine name ( Registers R00 thru R07 are to be initialized before executing "RKS2" )
• R01 = x0 •
R04 = h = stepsize
• R06 = y'0
R08 to R13: temp
• R02 = y0 •
R05 = N = number of steps
• R07 = z'0
• R03 = z0
Flags: /
Subroutine: a program which calculates 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.
| STACK | INPUTS | OUTPUTS |
| Z | / | z(x0+N.h) |
| Y | / | y(x0+N.h) |
| X | / | x0+N.h |
Example:
d2y/dx2 = -y.z
y(0) = 2 y'(0) = 1
Evaluate y(1) , z(1) , y'(1) , z'(1) with
h = 0.1 & N = 10
d2z/dx2 = x.(
y + z ) z(0) = 1
z'(0) = 1
| LBL "T"
RDN STO Z X<>Y CHS ST* Z - R^ * RTN |
"T" ASTO 00
0 STO 01
1 STO 06
2 STO 02
1 STO 07
1 STO 03
0.1 STO 04
10 STO 05
XEQ "RKS2" >>>> x =
1
( in 41 seconds )
RDN y(1) = 1.531358015 and
R06 = y'(1) = -2.312838895
RDN z(1) = 2.620254480
R07 = z'(1) = 2.941751649
-With h = 0.05 it yields y(1) =
1.531356736 y'(1) = -2.312840085
z(1) = 2.620254295 z'(1) =
2.941748608
-Actually, the exact results - rounded to 9D - are
y(1) = 1.531356646
y'(1) = -2.312840137
z(1) = 2.620254281
z'(1) = 2.941748399
-To continue the calculations, simply press R/S ( after changing
h/2 & N in registers R04 & R05 if needed )
RKS3: RK4 for Systems of 3 Second order Conservative
Equations
-"RKS3" solves the system:
y" = f(x,y,z,u)
y(x0) = y0 y'(x0)
= y'0
z" = g(x,y,z,u)
z(x0) = z0 z'(x0)
= z'0
u" = h(x,y,z,u)
u(x0) = u0 u'(x0)
= u'0
Data Registers: • R00: subroutine name ( Registers R00 thru R09 are to be initialized before executing "RKS3" )
• R01 = x0 •
R05 = h = stepsize
• R07 = y'0
• R02 = y0 •
R06 = N = number of steps
• R08 = z'0
• R03 = z0
R10 to R17: temp
• R09 = u'0
• R04 = u0
Flags: /
Subroutine: a program which calculates f(x;y;z;u)
in Z-register , g(x;y;z;u) in Y-register and h(x;y;z;u) in X-register
assuming x is in X-register , y is in Y-register , z is in Z-register and
u is in T-register upon entry.
| STACK | INPUTS | OUTPUTS |
| T | / | u(x0+N.h) |
| Z | / | z(x0+N.h) |
| Y | / | y(x0+N.h) |
| X | / | x0+N.h |
Example:
d2y/dx2 = -y.z.u
y(0) = 1 y'(0) = 1
Evaluate y(1) , z(1) , u(1) , y'(1) , z'(1) , u'(1)
with h = 0.1 & N = 10
d2z/dx2 = x.( y + z
- u ) z(0) = 1
z'(0) = 1
d2u/dx2 = x.y - z.u
u(0) = 2 u'(0) = 1
| LBL "T"
R^ CHS STO L R^ ST* Y ST+ L CLX RCL Z R^ ST+ L ST* Y X<> Z ST+ Y ST* Z X<> T LASTX * X<>Y RTN |
"T" ASTO 00
0 STO 01
1 STO 02 STO 03 STO 07
STO 08 STO 09
2 STO 04
0.1 STO 05
10 STO 06 XEQ "RKS3"
>>>> x = 1
( in 65 seconds )
RDN y(1) = 0.439528419
y'(1) = -2.101120400 = R07
RDN z(1) = 2.070938499
z'(1) = 1.269599239 = R08
RDN u(1) = 1.744522976
u'(1) = -1.704232092 = R09
-With h = 0.05 it yields
y(1) = 0.439524393
y'(1) = -2.101122784
z(1) = 2.070940521
z'(1) = 1.269597110
u(1) = 1.744524843
u'(1) = -1.704234567
-Actually, the exact results rounded to 9D are:
y(1) = 0.439524100
y'(1) = -2.101122880
z(1) = 2.070940654
z'(1) = 1.269596950
u(1) = 1.744524964
u'(1) = -1.704234756
-Press R/S to continue the calculations ( after changing h &
N in registers R05 & R06 if needed )
-The subroutine may be difficult to write, but "RKS3" is faster than
"RKSN" below, for 3 equations.
RKSN: RK4 for Systems of n Second order Conservative
Equations
-"RKSN" solves the system
y"1 = f1 ( x, y1 , ... , yn)
y1(x0) = y1,0
y'1(x0) = y'1,0
..................................
........................................
y"n = fn ( x, y1 , ... , yn)
yn(x0) = yn,0
y'n(x0) = y'n,0
Data Registers: • R00 = subroutine name ( Registers R00 thru R03 and R10 thru R10+2n are to be initialized before executing "RKSN" )
• R01 = n = number of equations = number of functions
R04 thru R09: temp
• R02 = h = stepsize
• R03 = N = number of steps
• R10 = x0
• R11 = y1(x0)
• R11+n = y'1(x0)
R11+2n thru R10+6n: temp
• R12 = y2(x0)
• R12+n = y'2(x0)
...............
......................
• R10+n = yn(x0) • R10+2n = y'n(x0)
( during the calculations, y'i are moved in registers R11+2n .... R10+3n )
Flags: /
Subroutine: A program which calculates
and stores f1 ( x, y1 , ... , yn)
, ... , fn ( x, y1 , ... , yn) in
registers R11+n, ... , R10+2n respectively
with x, y1, ... , yn in regiters R10, R11,
... , R10+n
| STACK | INPUTS | OUTPUTS |
| T | / | y3(x0+N.h) |
| Z | / | y2(x0+N.h) |
| Y | / | y1(x0+N.h) |
| X | / | x0+N.h |
Example:
d2y/dx2 = -y.z.u
y(0) = 1 y'(0) = 1
Evaluate y(1) , z(1) , u(1) , y'(1) , z'(1) , u'(1)
with h = 0.1 & N = 10
d2z/dx2 = x.( y + z - u )
z(0) = 1 z'(0) = 1
d2u/dx2 = x.y - z.u
u(0) = 2 u'(0) = 1
| 01 LBL "T"
02 RCL 11 03 RCL 12 04 RCL 13 05 * 06 * 07 CHS 08 STO 14 09 RCL 11 10 RCL 12 11 + 12 RCL 13 13 - 14 RCL 10 15 * 16 STO 15 17 RCL 10 18 RCL 11 19 * 20 RCL 12 21 RCL 13 22 * 23 - 24 STO 16 25 RTN |
"T" ASTO 00
Initial values: 0 STO 10
1 STO 11 STO 12 STO 14
STO 15 STO 16
2 STO 13
n = 3 STO 01 (
3 functions )
h = 0.1 STO 02
N = 10 STO 03 XEQ "RKSN"
>>>> x = 1
= R10
( in 2mn19s )
RDN y(1) = 0.439528419 = R11
y'(1) = -2.101120401 = R14
RDN z(1) = 2.070938499 = R12
z'(1) = 1.269599239 = R15
RDN u(1) = 1.744522977 = R13
u'(1) = -1.704232091 = R16
-Press R/S to continue the calculations ( after changing h &
N in registers R02 & R03 if need be )
-All n stage explicit Runge-Kutta formulae with error estimate are determined by ( n2 + 5.n - 2 ) / 2 coefficients ai,j ; bi ; ci ; di
k1 = h.f(x;y)
k2 = h.f(x+b2.h;y+a2,1.k1)
k3 = h.f(x+b3.h;y+a3,1.k1+a3,2.k2)
.............................................................
( actually, ai,1 + ai,2
+ .......... + ai,i-1 = bi )
kn = h.f(x+bn.h;y+an,1.k1+an,2.k2+ .......... +an,n-1.kn-1)
then, y(x+h) = y(x) + c1.k1+ ................ + cn.kn and the difference between the higher-order and the lower-order formulae gives an error estimate:
err-estim = d1.k1+ ................ + dn.kn
-"RKV" uses an 8-stage method with order 5 with a 6th-order error-estimating
formula.
Data Registers: • R00: f name ( Registers R00 to R05 are to be initialized before executing "RKV" )
• R01 = x0
R05 = y-error
R19 to R69 = the 51 coefficients given by "C51"
• R02 = y0
R06 to R10: temp
• R03 = h = stepsize
R11 to R18 = ki
• R04 = N = number of steps
Flags: /
Subroutine: a program calculating f(x;y)
assuming x is in X-register and y is in Y-register upon entry.
| STACK | INPUTS | OUTPUTS |
| Y | / | y(x0+N.h) |
| X | / | x0+N.h |
Example: Our usual equation:
y(0) = 1 and dy/dx = -2 x.y Evaluate
y(1) with h = 0.1 and N = 10
| 01 LBL "T"
02 * 03 ST+ X 04 CHS 05 RTN |
"T" ASTO 00
0 STO 01
1 STO 02
0.1 STO 03
10 STO 04 XEQ "RKV" yields
x = 1
( in 3mn56s )
X<>Y
y(1) = 0.367879457
and the error estimate >>> R05 = -13 10-9
-The error is actually 16 10-9
-Press R/S to continue ( after changing h and N in registers
R03 and R04 if need be ).
RKVB: RKV for Systems of 2 Differential Equations
Data Registers: • R00: f name ( Registers R00 to R06 are to be initialized before executing "RKVB" )
• R01 = x0
R06 = y-error
• R30 to R80 = the 51 coefficients given by "C51"
• R02 = y0
R07 = z-error
• R03 = y0
• R04 = h = stepsize
R08 to R13: temp
• R05 = N = number of steps
R14 to R21 = ki(y)
R22 to R29 = ki(z)
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.
| 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(1) and z(1) with h = 0.1 and
N = 10
| 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
0.1 STO 04
10 STO 05 XEQ "RKVB" yields
x = 1
( in 6mn34s )
RDN
y(1) = 0.367879378
RDN
z(1) = -0.735758757
and the error estimates RCL 06 =
-85 10-9 actually
y-error = -63 10-9
RCL 07 = 153 10-9
actually z-error = 126 10-9
-Press R/S to continue ( after changing h and N in registers R04
and R05 if needed ).
-"BST" solves the differential equation dy/dx = f(x,y)
with the initial value y(x0) = y0
-It uses an extrapolation to the limit and the modified midpoint formula
to compute y(x0+H)
-Let
z0 = y0
z1 = z0 + h.f(x0,z0)
z2 = z0 + 2h.f(x0+h,z1)
z3 = z1 + 2h.f(x0+2h,z2)
where h = H/n
..................................
zn = zn-2 + 2h.f(x0+(n-1).h,zn-1)
y(x0+H) ~ yn = (1/2).[ zn + zn-1 + h.f(x0+H,zn)
-The numbers yn are computed for n = 2 , 4 , 6 , ..... ,
16 ( at most ) and the Lagrange extrapolation polynomial is
used
until the estimated error is smaller than a specified tolerance
( tol )
-If the estimated error is still greater than tol with n = 16
, H is halved
-On the other hand, if the desired accuracy is satisfied with n <
8 , H is doubled
-Other ( and much better ) methods for stepsize control may be used
( cf "Numerical Recipes" )
-The Bulirsch-Stoer technique is quite similar to the Romberg integration:
the modified midpoint formula is second-order only,
but, for instance, with n = 16 and the deferred approach to
the limit, we actually use a 16th-order method !
Data Registers: • R00 = f name ( Registers R00 thru R03 are to be initialized before executing "BST" )
• R01 = x0 R04 = x1
R07 = h = H/n
R13 = next to last H
• R02 = y0 R05 = H
R08 to R12: temp R14, R15, ........
, R21 = y-approximations
• R03 = tol R06 = n = 2 , 4 , .... , 16
>>>> where tol is a small positive number that defines the desired accuracy
Flags: F09 - F10
Subroutine: A program that computes
f(x,y) assuming x is in X-register and y is in Y-register upon entry.
| STACK | INPUTS | OUTPUTS |
| Y | / | y(x1) |
| X | x1 | x1 |
Example: dy/dx =
x.(y/2)2
y(0) = 1 Evaluate
y(2) and y(2.5)
| LBL "T"
X<>Y 2 / X^2 * RTN |
"T" ASTO 00
0 STO 01
1 STO 02 and if we choose
tol = 10 -7 STO 03
2 XEQ "BST" >>>>
x = 2
( in 1mn49s )
y(1) has been computed with H = 1 , n = 10
RDN y(2) = 2.000000018
the exact result is 2
y(2) ------------------------ H = 1 , n = 14
2.5 R/S >>>>
x = 2.5
( in 3mn17s )
RDN y(2.5) = 4.571428682
the last 3 decimals should be 571
H = 1 has been replaced by 0.5 but the tolerance 10
-7 cannot be
the error is 1.11 10 -7
satisfied with this stepsize. So, H finally became 0.25 ( and n
= 14 )
-In fact, the solution is y = 1/(1-x2/8) so we'd need to increase tol in R03 and smaller and smaller stepsizes as x get closer to sqrt(8)
Note:
-Do not choose a too small tolerance , it could produce an
infinite loop.
BST2: BST for Systems of 2 Differential Equations
-The same method is employed to solve the system
dy/dx = f(x,y,z)
dz/dx = g(x,y,z)
with initial values y(x0) = y0
, z(x0) = z0
Data Registers: • R00 = f name ( Registers R00 thru R04 are to be initialized before executing "BST2" )
• R01 = x0 R05 = H
R09 to R17: temp
• R02 = y0 R06 = x1
R18 = next to last H-value
• R03 = z0 R07 = n = 2 ,
4 , .... , 16 R19, R20, ........ , R26 =
y-approximations
• R04 = tol R08 = h = H/n
R27, R28, ........ , R34 = z-approximations
>>>> where tol is a small positive number that defines the desired accuracy
Flags: F09 - F10
Subroutine: A program that computes
f(x,y,z) in Y-register
and g(x,y,z) in Z-register assuming x is in X-register
, y is in Y-register and z is in Z-register upon entry.
| STACK | INPUTS | OUTPUTS |
| Z | / | z(x1) |
| Y | / | y(x1) |
| X | x1 | x1 |
Example: d2y/dx2
= -2y - 2x.dy/dx y(0) = 1 y'(0) = 0
Calculate y(1) and y'(1)
[ the exact solution is y = exp ( -x2 ) ]
This second-order equation is equivalent to the system
dy/dx = z
y(0) = 1
dz/dx = -2y - 2x.z
z(0) = 0
| LBL "T"
RCL Z * + ST+ X CHS RTN |
"T" ASTO 00
0 STO 01
1 STO 02
0 STO 03 and with tol = 10
-7 STO 04
1 XEQ "BST2" >>>>
x = 1
( in 2mn08s )
RDN
y(1) = 0.367879446
the last decimal should be a 1
RDN z(1) = y'(1) = -0.735758909
the last 3 decimals should be 882 z-error
~ 27 E-9 < E-7
BSTN: BST for Systems of n Differential Equations
dy1/dx = f1(x,y1,
... ,yn)
dy2/dx = f2(x,y1,
... ,yn)
..................
with the initial values: yi(x0)
i = 1 , ... , n
dyn/dx = fn(x,y1,
... ,yn)
Data Registers: • R00 = subroutine name ( Registers R00-R01-R02 and R20 thru R20+n are to be initialized before executing "BSTN" )
• R01 = n = number of equations = number of functions R03 thru R18: temp R19: unused
• R02 = tolerance = a small positive number that specifies the desired accuracy
• R20 = x0
• R21 = y1(x0)
R21+n thru R20+14n: temp
• R22 = y2(x0)
............... ( initial
values )
• R20+n = yn(x0)
Flags: F09-F10
Subroutine: A program which calculates
and stores f1 ( x, y1 , ... , yn)
, ... , fn ( x, y1 , ... , yn) in
registers R21+n, ... , R20+2n respectively
with x, y1, ... , yn in regiters R20, R21,
... , R20+n
| STACK | INPUTS | OUTPUTS |
| T | / | y3(x1) |
| Z | / | y2(x1) |
| Y | / | y1(x1) |
| X | x1 | x1 |
Example:
dy/dx = -y.z.u
y(0) = 1 Evaluate
y(1) , z(1) , u(1) with a tolerance = 10 -7
dz/dx = x.( y + z - u )
z(0) = 1
du/dx = x.y - z.u
u(0) = 2
| 01 LBL "T"
02 RCL 21 03 RCL 22 04 RCL 23 05 * 06 * 07 CHS 08 STO 24 09 RCL 21 10 RCL 22 11 + 12 RCL 23 13 - 14 RCL 20 15 * 16 STO 25 17 LASTX 18 RCL 21 19 * 20 RCL 22 21 RCL 23 22 * 23 - 24 STO 26 25 RTN |
"T" ASTO 00
Initial values: 0 STO 20
1 STO 21 STO 22
2 STO 23
n = 3 STO 01
( 3 functions )
tol = E-7 STO 02
1 XEQ "BSTN" >>>>
x = 1
= R20
( in 12mn06s )
RDN y(1) = 0.258207909 = R21
y(1) = 0.258207906
RDN z(1) = 1.157623986 = R22
the exact results, rounded to 9D are
z(1) = 1.157623981
RDN u(1) = 0.842178304 = R23
u(1) = 0.842178312
-The long execution time is due to the fact that the method is first
used with H = 1
but the desired accuracy is not satisfied. So the stepsize is
divided by 2.
-To compute y(2) z(2) u(2) key in 2 R/S it yields ( in 6m46s )
y(2) = 0.106363294
y(2) = 0.106363329
z(2) = 3.886706181
the exact results are
z(2) = 3.886706159
( rounded to 9D )
u(2) = 0.196515847
u(2) = 0.196515847
Notes:
-The initial stepsize is always +/-1 ( line 03 )
-The modified midpoint rule ( and Richarson extrapolation ) are used
with h/2 , h/4 , h/6 , ..... , h/16 ( at most ) until the desired
accuracy is satisfied.
Stoermer's rule for Conservative 2nd-order Equations
-"STRM" solves the differential equation d2y/dx2
= f(x,y) with initial values
y(x0) = y0 , y'(x0)
= y'0
where the first derivative does not appear explicitly.
-This kind of equations may be re-written as a system which could be solved by "BST2" but the following method is faster.
-Like "BST", "STRM" also uses the deferred approach to the limit but
- instead of the modified midpoint method - the Stoermer's rule is employed:
-Let
d0 = h.[ y'0 + (h/2).f(x0,y0)
]
y1 = y0 + d0
dk = dk-1 + h2f(x0+k.h,yk)
dk = yk+1 - yk
k = 1 , 2 , ..... , n-1
h = H/n
y'n = dn-1/h + (h/2).f(x0+H,yn)
Data Registers: • R00 = f name ( Registers R00 thru R04 are to be initialized before executing "STRM" )
• R01 = x0 R05 = H
R09 to R14: temp
• R02 = y0 R06 = x1
R15 = next to last H-value
• R03 = y'0 R07 = n = 2 , 4 ,
.... , 16 R16, R17, ........ , R23 = y-approximations
• R04 = tol R08 = h = H/n
R24, R25, ........ , R31 = y'-approximations
>>>> where tol is a small positive number that defines the desired accuracy
Flags: F09 - F10
Subroutine: A program that computes
f(x,y) assuming x is in X-register and y is in Y-register upon
entry.
| STACK | INPUTS | OUTPUTS |
| Z | / | y'(x1) |
| Y | / | y(x1) |
| X | x1 | x1 |
Example: d2y/dx2
= -y.( x2 + y2 )1/2
y(0) = 1 , y'(0) = 0 Evaluate
y(1) and y(PI)
| LBL "T"
X^2 RCL Y X^2 + SQRT * CHS RTN |
"T" ASTO 00
0 STO 01
1 STO 02
0 STO 03 and with tol = 10
-7 STO 04
1 XEQ "STRM" >>>>
x = 1
( in 53 seconds )
RDN y(1) = 0.536630616
the 9 decimals are correct
RDN y'(1) = -0.860171925
the last decimal should be a 7
PI R/S
>>>> x = 3.141592654
( in 2mn24s )
RDN y(PI) = -0.411893053
the 9 decimals are exact
RDN y'(PI) = 1.018399901
the last decimal should be a 3
2nd order differential equations ( 4th order method )
-"SRK4" solves the equation y" = f(x,y,y')
with the initial values: y(x0) = y0
, y'(x0) = y'0
Data Registers: • R00 = function name ( Registers R00 thru R05 are to be initialized before executing "SRK4" )
• R01 = x0
• R02 = y0 •
R04 = h = stepsize
• R03 = y'0 •
R05 = m = number of steps
R06 thru R10: temp
Flags: /
Subroutine:
A program which computes y" = f(x,y,y') assuming x ,
y , y' are in registers X , Y , Z ( respectively ) upon entry
| STACK | INPUTS | OUTPUTS |
| Z | / | y'(x0+m.h) |
| Y | / | y(x0+m.h) |
| X | / | x0+m.h |
-Simply press R/S to continue with the same h and m
Example: Let's consider the Lane-Emden equation of index 3 y" = -(2/x) y' - y3 with the initial values y(0) = 1 , y'(0) = 0
-There is already a difficulty because x = 0 is a singular point!
-But if we use a series expansion, we find after substitution that
y = 1 + a.x2 + .... will satisfy the LEE if
a = -1/6 whence y"(0) = -1/3
-So, we can load the following subroutine into program memory:
| LBL "T"
X=0? GTO 00 RCL Z ST+ X X<>Y / X<>Y 3 Y^X + CHS RTN LBL 00 3 1/X CHS RTN |
-If we want to estimate y(1) with a stepsize h = 0.1
( whence m = 10 )
"T" ASTO 00
0 STO 01 STO 03
0.1 STO 04 = h
1 STO 02
10 STO 05
XEQ "SRK4" >>>>
x = 1
( in 56 seconds )
RDN y(1) = 0.855057170
RDN y'(1) = -0.252129561
-These new "initial" values are also stored in registers R01 R02 R03
-With h = 0.05 and m = 20 we would have found
y(1) = 0.855057539 & y'(1) = -0.252129276
Notes: The solution of the Lane-Emden Equation of index 3 looks like this:
y
1|*
I
|
*
|
*
|
*
|
*
|-----------------------------------------------------------*x1--------
x
0
y(x1) = 0 for x1 = 6.896848619
and y'(x1) = -0.0424297576
and there is an inflexion point I with xI = 1.495999168
, yI = 0.720621687 and y'(xI) = -0.279913175
-The solutions of the Lane-Emden Equations of index n (
y" + (2/x).y' + yn = 0 , y(0)= 1 , y'(0) = 0 )
can be expressed by elementary functions for only 3 n-values:
a) n = 0 y(x) = 1 - x2/6
b) n = 1 y(x) = (sin x)/x
c) n = 5 y(x) = ( 1+x2/3)
-1/2
3rd order differential equations ( 4th order method )
-Likewise, "TRK4" solves the equation y"' = f(x,y,y',y") with the initial values: y(x0) = y0 , y'(x0) = y'0 , y"(x0) = y"0
Data Registers: • R00 = function name ( Registers R00 thru R06 are to be initialized before executing "TRK4" )
• R01 = x0
• R02 = y0
• R03 = y'0 •
R05 = h = stepsize
• R04 = y"0 • R06
= m = number of steps
R07 thru R13: temp
Flags: /
Subroutine:
A program which computes y"' = f(x,y,y',y") assuming
x , y , y' , y" are in registers X , Y , Z , T ( respectively )
upon entry
| STACK | INPUTS | OUTPUTS |
| T | / | y"(x0+m.h) |
| Z | / | y'(x0+m.h) |
| Y | / | y(x0+m.h) |
| X | / | x0+m.h |
-Simply press R/S to continue with the same h and m
Example: y"' = 2xy" - x2y' + y2 with y(0) = 1 , y'(0) = 0 , y"(0) = -1
-We load, for instance:
| LBL "S"
X^2 ST* Z X<> L ST+ X ST* T RDN X^2 - - RTN |
-With h = 0.1 and m = 10
"S" ASTO 00
0 STO 01 STO 03
0.1 STO 05 = h
1 STO 02 CHS STO 04
10 STO 06
XEQ "TRK4" >>>>
x = 1
( in 49 seconds )
RDN y(1) = 0.595434736
RDN y'(1) = -0.776441445
RDN y"(1) = -0.791715205
-With h = 0.05 and m = 20, we would find
y(1) = 0.595431304 , y'(1) = -0.776444326
, y"(1) = -0.791718300
N-th order differential equations ( 4th order method )
-The differential equation is now y(n) = f(x,y,y',y",.....,y(n-1))
with the initial values: y(x0) = y0
, y'(x0) = y'0 , ........
, y(n-1)(x0) = y(n-1)0
Data Registers: • R00 = function name ( Registers R00 thru R03 and R09 thru Rn+9 are to be initialized before executing "NRK4" )
• R01 = n ( order )
• R02 = h = stepsize
R04 thru R08 & Rn+10 thru R3n+9: temp
• R03 = m = number of steps
• R09 = x0 • R10 = y0
• R11 = y'0 • R12 = y"0
..................... • Rn+9 = y(n-1)0
Flags: /
Subroutine: A program which computes
y(n) = f(x,y,y',y",.....,y(n-1)) assuming
x , y , y' , y" , ........ , y(n-1) are in registers
R09 R10 R11 R12 .... Rnn+9
| STACK | INPUTS | OUTPUTS |
| T | / | y"(x0+m.h) |
| Z | / | y'(x0+m.h) |
| Y | / | y(x0+m.h) |
| X | / | x0+m.h |
-Simply press R/S to continue with the same h and m
Example: y(5) = y(4) - 2x.y"' + y" - y.y' with y(0) = 1 , y'(0) = y"'(0) = y(4)(0) = 0 , y"(0) = -1
-We load, for instance:
| LBL "W"
RCL 14 RCL 09 ST+ X RCL 13 * - RCL 12 + RCL 10 RCL 11 * - RTN |
-With h = 0.1 and m = 10
"W" ASTO 00
and the initial values: 0 STO 09
STO 11 STO 13 STO 14
1 STO 10 CHS STO 12
5 STO 01
( fifth order equation )
0.1 STO 02
( h )
10 STO 03
( m = 10 )
XEQ "NRK4" >>>>
x = 1
= R09
( in 2mn48s )
RDN y(1) = 0.491724880
= R10
RDN y'(1) = -1.041200697 = R11
and RCL 13 = y"'(1) = -0.479803795
RDN y"(1) = -1.163353624 = R12
RCL 14 = y(4)(1) = -0.897595630
-With h = 0.05 and m = 20, it yields:
y(1) = 0.491724223 , y'(1) = -1.041200398
, y"(1) = -1.163353549 , y"'(1) = -0.479804004
, y(4)(1) = -0.897594479
References:
[1] Abramowitz and Stegun , "Handbook of Mathematical Functions"
- Dover Publications - ISBN 0-486-61272-4
[2] Press , Vetterling , Teukolsky , Flannery - "Numerical Recipes
in C++" - Cambridge University Press - ISBN 0-521-75033-4
[3] F. Scheid - "Numerical Analysis" - Mc Graw Hill
ISBN 07-055197-9
[4] J.C. Butcher - "The Numerical Analysis of Ordinary
Differential Equations" - John Wiley & Sons - ISBN
0-471-91046-5