hp41programs

Special Runge-Kutta

Runge-Kutta Methods for Special ODEs for the HP-41


Overview
 

 1°)  Third-order ODEs  y''' = f(x,y,y')
 2°)  Fourth-order ODEs  y'''' = f(x,y)
 

-All these differential equations may be solved by usual Runge-Kutta formulae.
-But the following programs use special methods with less stages to get the same accuracy
 or the same number of stages to get a better accuracy.
 
 

1°)  Third-order ODEs  y''' = f(x,y,y')
 
 

-This routine employs a 4-stage 5th-order Runge-Kutta method given in reference [1]
 
 

Data Registers:           •  R00 = function name              ( Registers R00 thru R06 are to be initialized before executing "RKTG" )

                                      •  R01 = x          •  R04 = y"        R07 to R11: temp
                                      •  R02 = y          •  R05 = h
                                      •  R03 = y'         •  R06 = N
Flags: /
Subroutine:  A program that takes  x , y , y'  in registers  X , Y , Z  respectively and returns  f(x,y,y')  in X-register.
 

-Line 209 is a three-byte GTO 01
 
 

 01 LBL "RKTG"
 02 RCL 06
 03 STO 11
 04 LBL 01
 05 RCL 03
 06 RCL 02
 07 RCL 01        
 08 XEQ IND 00
 09 STO 07
 10 RCL 05
 11 *
 12 50
 13 /
 14 RCL 04
 15 5
 16 /
 17 +
 18 RCL 05
 19 *
 20 RCL 03
 21 +
 22 RCL 04
 23 RCL 05
 24 *
 25 10
 26 /
 27 RCL 03
 28 +
 29 RCL 05
 30 *
 31 5		  32 /
 33 RCL 02
 34 +
 35 RCL 05
 36 5
 37 /
 38 RCL 01
 39 +
 40 XEQ IND 00
 41 STO 08
 42 7
 43 *
 44 RCL 07
 45 -
 46 27
 47 /
 48 RCL 05
 49 *
 50 RCL 04
 51 ST+ X
 52 3
 53 /
 54 +
 55 RCL 05
 56 *
 57 RCL 03
 58 +
 59 RCL 08
 60 301
 61 *
 62 RCL 07		  63 49
 64 *
 65 -
 66 RCL 05
 67 *
 68 4860
 69 /
 70 RCL 04        
 71 4
 72 *
 73 18
 74 /
 75 +
 76 RCL 05
 77 *
 78 RCL 03
 79 ST+ X
 80 3
 81 /
 82 +
 83 RCL 05
 84 *
 85 RCL 02
 86 +
 87 RCL 05
 88 ST+ X
 89 3
 90 /
 91 RCL 01 
 92 +
 93 XEQ IND 00		  94 STO 09
 95 9
 96 *
 97 RCL 08
 98 ST+ X
 99 -
100 35
101 /
102 RCL 07        
103 .3
104 *
105 +
106 RCL 05
107 *
108 RCL 04
109 +
110 RCL 05
111 *
112 RCL 03
113 +
114 RCL 09
115 ST+ X
116 RCL 08
117 -
118 RCL 07
119 7
120 *
121 +
122 RCL 05 
123 *
124 50		 125 /
126 RCL 04
127 2
128 /
129 +
130 RCL 05
131 *
132 RCL 03        
133 +
134 RCL 05
135 *
136 RCL 02
137 +
138 RCL 01
139 RCL 05
140 +
141 XEQ IND 00
142 STO 10
143 RCL 07
144 7
145 *
146 STO 07
147 RCL 08
148 40
149 *
150 +
151 RCL 09
152 9
153 *
154 +		 155 RCL 05
156 *
157 336
158 /
159 RCL 04
160 2
161 /
162 +
163 RCL 05        
164 *
165 RCL 03
166 +
167 RCL 05
168 *
169 ST+ 02
170 RCL 07
171 RCL 08 
172 50
173 *
174 +
175 RCL 09
176 27
177 *
178 +
179 RCL 05
180 *
181 168
182 /
183 RCL 04
184 +		 185 RCL 05
186 *
187 ST+ 03
188 RCL 07
189 ST+ X
190 RCL 08        
191 125
192 *
193 +
194 RCL 09
195 162
196 *
197 +
198 RCL 10 
199 35
200 *
201 +
202 RCL 05
203 ST+ 01
204 *
205 336
206 /
207 ST+ 04
208 DSE 11
209 GTO 01
210 RCL 04
211 RCL 03
212 RCL 02
213 RCL 01
214 END

 
      ( 267 bytes / SIZE 012 )
 
 

      STACK        INPUTS      OUTPUTS
           T             /          y"(x)
           Z             /          y'(x)
           Y             /          y(x)
           X             /            x

 
Example:     y''' = 8 x y + ( 4 x2 - 2 ) y'        y(0) = 1    y'(0) = 0   y"(0) = -2          Calculate  y(1)

-Load for example this routine in main memory:
 
 

01 LBL "T"
02 ST* Y
03 X^2
04 ST+ X
05 1
06 -
07 ST+ X
08 R^
09 *
10 X<>Y
11 8
12 *
13 +
14 RTN

 
 "T"  ASTO 00

  0  STO 01  STO 03
  1  STO 02  -2  STO 04

-If you choose  h = 0.05  &  N = 20

  0.05  STO 05  20  STO 06   XEQ "RKTG"  >>>>   x    =  1                     = R01
                                                                       RDN  y(1) = 0.367879436   = R02
                                                                       RDN y'(1) = -0.735758901  = R03
                                                                       RDN y"(1) = 0.735758853   = R04

Notes:

-The exact solution is  y = exp(-x2)  so  y(1) = 0.367879441....
-Simply press R/S to continue with the same h & N
 

2°)  Fourth-order ODEs  y'''' = f(x,y)
 
 

-This routine employs a 3-stage 4th-order Runge-Kutta method given in reference [2]
 

Data Registers:           •  R00 = function name              ( Registers R00 thru R07 are to be initialized before executing "RKFD" )

                                      •  R01 = x          •  R04 = y"          •  R07 = N        R08 to R11: temp
                                      •  R02 = y          •  R05 = y'''
                                      •  R03 = y'         •  R06 = h
Flags: /
Subroutine:  A program that takes  x , y  in registers  X , Y  respectively and returns  f(x,y)  in X-register.
 

-Line 182 is a three-byte GTO 01
 
 

 01 LBL "RKFD"
 02 RCL 07
 03 STO 11
 04 LBL 01
 05 RCL 02
 06 RCL 01
 07 XEQ IND 00
 08 STO 08
 09 CHS
 10 RCL 06
 11 *
 12 5
 13 /
 14 RCL 05
 15 64
 16 *
 17 7986
 18 /
 19 +
 20 RCL 06
 21 *
 22 RCL 04
 23 8
 24 *
 25 121
 26 /
 27 +		  28 RCL 06
 29 *
 30 RCL 03
 31 4
 32 *
 33 11
 34 /
 35 +
 36 RCL 06
 37 *
 38 RCL 02
 39 +
 40 RCL 06
 41 4
 42 *
 43 11
 44 /
 45 RCL 01
 46 +
 47 XEQ IND 00
 48 STO 09
 49 RCL 08
 50 +
 51 19
 52 *
 53 RCL 06
 54 *		  55 125
 56 /
 57 RCL 05 
 58 4913
 59 *
 60 48 E3
 61 /
 62 +
 63 RCL 06        
 64 *
 65 RCL 04
 66 289
 67 *
 68 800
 69 /
 70 +
 71 RCL 06
 72 *
 73 RCL 03
 74 .85
 75 *
 76 +
 77 RCL 06
 78 *
 79 RCL 02
 80 +
 81 RCL 06		  82 .85
 83 *
 84 RCL 01
 85 +
 86 XEQ IND 00
 87 STO 10
 88 20
 89 /
 90 RCL 09
 91 7
 92 *
 93 75
 94 /
 95 -
 96 RCL 08
 97 .085
 98 *
 99 +
100 RCL 06
101 *
102 RCL 05
103 6
104 /
105 +
106 RCL 06
107 *
108 RCL 04		 109 2
110 /
111 +
112 RCL 06
113 *
114 RCL 03        
115 +
116 RCL 06
117 *
118 ST+ 02
119 RCL 10
120 5
121 *
122 RCL 09
123 209
124 *
125 +
126 1926
127 /
128 RCL 08
129 18
130 /
131 +
132 RCL 06
133 *
134 RCL 05
135 2		 136 /
137 +
138 RCL 06        
139 *
140 RCL 04 
141 +
142 RCL 06
143 *
144 ST+ 03
145 RCL 10
146 2400
147 *
148 RCL 09
149 14399
150 *
151 +
152 RCL 08
153 5029
154 *
155 STO 08
156 +
157 43656
158 /
159 RCL 06
160 *
161 RCL 05		 162 +
163 RCL 06        
164 *
165 ST+ 04
166 RCL 10 
167 16 E3
168 *
169 RCL 09
170 22627
171 *
172 +
173 RCL 08
174 +
175 43656
176 /
177 RCL 06
178 ST+ 01
179 *
180 ST+ 05
181 DSE 11
182 GTO 01
183 RCL 04
184 RCL 03
185 RCL 02
186 RCL 01
187 END

 
    ( 268 bytes / SIZE 012 )
 
 

      STACK        INPUTS      OUTPUTS
           T             /          y"(x)
           Z             /          y'(x)
           Y             /          y(x)
           X             /            x

 
Example:       y''''(x) = ( 16 x4 - 48 x2 + 12 ) y        y(0) = 1    y'(0) = 0   y"(0) = -2   y'''(0) = 0          Calculate  y(1)
 

-Load for example this routine in main memory:
 
 

01 LBL "T"
02 X^2
03 RCL X
04 16
05 *
06 48
07 -
08 *
09 12
10 +
11 *
12 RTN

 
 "T"  ASTO 00

  0  STO 01  STO 03  STO 05
  1  STO 02  -2  STO 04

-If you choose  h = 0.05  &  N = 20

  0.05  STO 06  20  STO 07   XEQ "RKFD"  >>>>   x    =  1                     = R01
                                                                       RDN  y(1) = 0.367879077   = R02
                                                                       RDN y'(1) = -0.735759984  = R03
                                                                       RDN y"(1) = 0.735756781   = R04

  and   R05 = y'''(1) = 1.471518787
 

Notes:

-The exact solution is  y = exp(-x2)  so  y(1) = 0.367879441....
-Simply press R/S to continue with the same h & N
 
 

References:

[1]  Firas Adel Fawzi, Norazak Senu, Fudziah Ismail, Zanariah Abd. Majid - "An Efficient of Direct Integrator of Runge-Kutta Type Method
       for solving y''' = f(x,y,y') with Application to Thin Film Flow Problem"
[2]  Kasim Hussain, Fudziah Ismail, Norazak Senu - "Runge-Kutta Type Methods for Directly Solving Special Fourth-Order Ordinary Differential Equations"