Optimal Runge-Kutta Method of Order 4 for the HP-41

Overview
 

-This program solves differential equations:   dy/dx = f(x,y)

  systems of 2 differential equations:   dy/dx = f(x,y,z)
                                                        dz/dx = g(x,y,z)

  and systems of 3 differential equations:  dy/dx = f(x,y,z,t)
                                                              dz/dx = g(x,y,z,t)
                                                              dt/dx = h(x,y,z,t)

  by an optimal Runge-Kutta method of order 4 given in reference [1] page 280.

-It is a 4-step method:

    k₁ = h.f(x;y)
    k₂ = h.f(x+b₂.h;y+a_2,1.k₁)
    k₃ = h.f(x+b₃.h;y+a_3,1.k₁+a_3,2.k₂)
    k₄ = h.f(x+b₄.h;y+a_4,1.k₁+a_4,2.k₂+a_4,3.k₃)

then,     y(x+h) = y(x) + c₁.k₁+ ........... + c₄.k₄

-And the coefficients are:

          0.3716151060     0.3716151060
          0.6                     -0.1180444797        0.7180444797
          1                          0.5173871366       -0.5608902997         1.043503163

          0.1474734369     0.3125088197        0.3903768538         0.1496408895
 

Program Listing
 

-The coefficients - except 1 - are stored in registers R16 to R26 ( line 02 to 25 )
-Then the differential equation is solved.

-Clear flags F02 & F03 for a simple differential equation    CF 02  CF 03
-Set F02 & clear F03 for a system of 2 ODEs                   SF 02  CF 03
-Set F03 for a system of 3 differential equations                 SF 03
 

Data Registers:           •  R00 = function name               ( Registers R00 thru R04 or R05 or R06 are to be initialized before executing "RK4OP" )

    1 equation                 •  R01 = x          •  R03 = h                            R05 to R08: temp    R09 to R15: unused   R16 to R26: coeff   R27-R28: unused
                                      •  R02 = y          •  R04 = N
                                      •  R03 = z

    2 equations               •  R01 = x          •  R04 = h                            R06 to R12: temp    R13 to R15: unsed   R16 to R26: coeff    R27-R28: unused
                                      •  R02 = y          •  R05 = N
                                      •  R03 = z

    3 equations               •  R01 = x          •  R05 = h                             R07 to R28: temp   ( R16 to R26: coeff )
                                      •  R02 = y          •  R06 = N
                                      •  R03 = z
                                      •  R04 = t

Flags:   CF 02  CF 03  <->  1 differential equation
              SF 02  CF 03  <->  2 differential equations
                   SF 03         <->  3 differential equations

Subroutine:   CF 02  CF 03  A program that takes x , y in registers  X , Y  and returns  f(x,y)  in X-register
                        SF 02  CF 03  ---------------------  x , y , z  in  X , Y , Z    -----------  f(x,y,z) in Y-register and g(x,y,z) in X-registe
                             SF 03         ---------------------  x , y , z , t  in X , Y , Z , T  -------  f(x,y,z,t) , g(x,y,z,t) , h(x,y,z,t)  in registers  Z , Y , X  respectively.
 
 

-Lines 27-234-239-428-434 are three-byte GTOs
 
 

01 LBL "RK4OP"  02 .371615106  03 STO 16  04 .6  05 STO 17  06 .1180444797  07 CHS  08 STO 18  09 -  10 STO 19  11 .5173871366  12 STO 20  13 .5608902997  14 CHS  15 STO 21  16 1.043503163  17 STO 22  18 .1474734369  19 STO 23  20 .3125088197  21 STO 24  22 .3903768538  23 STO 25  24 .1496408895  25 STO 26  26 FS? 03  27 GTO 03  28 FS? 02  29 GTO 02  30 LBL 01  31 RCL 04  32 STO 08  33 LBL 11  34 RCL 02  35 RCL 01  36 XEQ IND 00  37 RCL 03  38 *  39 STO 05  40 RCL 16  41 *  42 RCL 02  43 +  44 RCL 03  45 RCL 16  46 *  47 RCL 01  48 +  49 XEQ IND 00  50 RCL 03  51 *  52 STO 06  53 RCL 19  54 *  55 RCL 05  56 RCL 18  57 *  58 +  59 RCL 02  60 +  61 RCL 03  62 RCL 17

63 *  64 RCL 01  65 +  66 XEQ IND 00   67 RCL 03  68 *  69 STO 07  70 RCL 22  71 *  72 RCL 06  73 RCL 21  74 *  75 +  76 RCL 05  77 RCL 20  78 *  79 +  80 RCL 02  81 +  82 RCL 03  83 RCL 01  84 +  85 XEQ IND 00  86 RCL 03  87 ST+ 01  88 *  89 RCL 26  90 *  91 RCL 07  92 RCL 25  93 *  94 +  95 RCL 06  96 RCL 24  97 *  98 +  99 RCL 05 100 RCL 23 101 * 102 + 103 ST+ 02 104 DSE 08 105 GTO 11 106 RCL 02 107 RCL 01 108 RTN 109 GTO 01 110 LBL 02 111 RCL 05 112 STO 12 113 LBL 12 114 RCL 03 115 RCL 02 116 RCL 01 117 XEQ IND 00 118 RCL 04 119 ST* Z 120 * 121 STO 09 122 RCL 16 123 * 124 RCL 03

125 + 126 X<>Y 127 STO 06 128 RCL 16 129 * 130 RCL 02 131 + 132 RCL 04 133 RCL 16 134 * 135 RCL 01 136 + 137 XEQ IND 00 138 RCL 04 139 ST* Z 140 * 141 STO 10 142 RCL 19 143 * 144 RCL 09 145 RCL 18 146 * 147 + 148 RCL 03 149 + 150 X<>Y 151 STO 07 152 RCL 19 153 * 154 RCL 06 155 RCL 18 156 * 157 + 158 RCL 02 159 + 160 RCL 04 161 RCL 17 162 * 163 RCL 01 164 + 165 XEQ IND 00 166 RCL 04 167 ST* Z 168 * 169 STO 11 170 RCL 22 171 * 172 RCL 10 173 RCL 21 174 * 175 + 176 RCL 09 177 RCL 20 178 * 179 + 180 RCL 03 181 + 182 X<>Y 183 STO 08 184 RCL 22 185 * 186 RCL 07

187 RCL 21 188 * 189 + 190 RCL 06 191 RCL 20 192 * 193 + 194 RCL 02 195 + 196 RCL 01 197 RCL 04 198 + 199 XEQ IND 00 200 RCL 04 201 ST+ 01 202 RCL 26 203 * 204 ST* Z 205 * 206 RCL 11 207 RCL 25 208 * 209 + 210 RCL 10 211 RCL 24 212 * 213 + 214 RCL 09 215 RCL 23 216 * 217 + 218 ST+ 03 219 X<> 08 220 RCL 25 221 * 222 + 223 RCL 07 224 RCL 24 225 * 226 + 227 RCL 06 228 RCL 23 229 * 230 + 231 ST+ 02 232 DSE 12 233 GTO 12 234 RCL 03 235 RCL 02 236 RCL 01 237 RTN 238 GTO 02 239 LBL 03 240 RCL 06 241 STO 27 242 LBL 13 243 RCL 04 244 RCL 03 245 RCL 02 246 RCL 01 247 XEQ IND 00 248 RCL 05

249 ST* T 250 ST* Z 251 * 252 STO 13 253 RCL 16 254 * 255 RCL 04 256 + 257 X<>Y 258 STO 10 259 RCL 16 260 * 261 RCL 03 262 + 263 R^ 264 STO 07 265 RCL 16 266 * 267 RCL 02 268 + 269 STO 15 270 CLX 271 RCL 05 272 RCL 16 273 * 274 RCL 01 275 + 276 RCL 15 277 X<>Y 278 XEQ IND 00 279 RCL 05 280 ST* T 281 ST* Z 282 * 283 STO 15 284 RDN 285 STO 11 286 RCL 19 287 * 288 RCL 10 289 RCL 18 290 * 291 + 292 RCL 03 293 + 294 X<>Y 295 STO 08 296 RCL 19 297 * 298 RCL 07 299 RCL 18 300 * 301 + 302 RCL 02 303 + 304 X<> 15 305 STO 14 306 RCL 19 307 * 308 RCL 13 309 RCL 18 310 *

311 + 312 RCL 04 313 + 314 X<>Y 315 RCL 05 316 RCL 17 317 * 318 RCL 01 319 + 320 RCL 15 321 X<>Y 322 XEQ IND 00 323 RCL 05 324 ST* T 325 ST* Z 326 * 327 STO 28 328 RDN 329 STO 12 330 RCL 22 331 * 332 RCL 11 333 RCL 21 334 * 335 + 336 RCL 10 337 RCL 20 338 * 339 + 340 RCL 03 341 + 342 X<>Y 343 STO 09 344 RCL 22 345 * 346 RCL 08 347 RCL 21 348 * 349 + 350 RCL 07 351 RCL 20 352 * 353 + 354 RCL 02 355 + 356 X<> 28 357 STO 15 358 RCL 22 359 * 360 RCL 14 361 RCL 21 362 * 363 + 364 RCL 13 365 RCL 20 366 * 367 + 368 RCL 04 369 + 370 X<>Y 371 RCL 05 372 RCL 01

373 + 374 RCL 28 375 X<>Y 376 XEQ IND 00 377 RCL 05 378 ST+ 01 379 ST* T 380 ST* Z 381 * 382 RCL 26 383 ST* T 384 ST* Z 385 * 386 X<> 12 387 RCL 25 388 * 389 + 390 RCL 11 391 RCL 24 392 * 393 + 394 RCL 10 395 RCL 23 396 * 397 + 398 ST+ 03 399 X<> 09 400 RCL 25 401 * 402 + 403 RCL 08 404 RCL 24 405 * 406 + 407 RCL 07 408 RCL 23 409 * 410 + 411 ST+ 02 412 RCL 12 413 RCL 15 414 RCL 25 415 * 416 + 417 RCL 14 418 RCL 24 419 * 420 + 421 RCL 13 422 RCL 23 423 * 424 + 425 ST+ 04 426 DSE 27 427 GTO 13 428 RCL 04 429 RCL 03 430 RCL 02 431 RCL 01 432 RTN 433 GTO 03 434 END

 
     ( 669 bytes / SIZE 027 or 029 )
 
 

STACK	INPUT	OUTPUTS1	OUTPUTS2	OUTPUTS3
T	/	/	/	t
Z	/	/	z	z
Y	/	y	y	y
X	/	x	x	x

 
Example1:    y' = -2 x.y   y(0) = 1  y'(0) = 0      Evaluate  y(1)
 
 

01  LBL "T"  02  *  03  ST+ X  04  CHS  05  RTN

 
  CF 02  CF 03

 "T"  ASTO 00

  0     STO 01
  1     STO 02    and if we choose  h = 0.1
 0.1   STO 03
 10    STO 04    XEQ "RK4OP"   yields          x = 1                                       ---Execution time = 35s---
                                X<>Y                        y(1) = 0.367879270

-The exact result is 1/e =  0.367879441...

-Press R/S to continue with the same h and N or modify R03-R04
 

Example2:     y(0) = 1 ; y'(0) = 0 ;  y'' + 2.x.y' + 2.y = 0   Evaluate y(1) and y'(1)

-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

 
  SF 02  CF 03

 "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 "RK4OP"   yields       x =  1                                               ---Execution time = 55s---
                                  RDN                      y(1) =  0.367879816
                                  RDN                      z(1) = -0.735759631

-The exact results are  y(1) = 1/e =  0.367879441... & z(1) = -2/e = -0.735758883...

-Press R/S to continue with the same h and N or modify  R04-R05

Example3:

     dy/dx = -y.z.t                    y(0) = 1         Evaluate  y(1) , z(1) , t(1)  with  h = 0.1   ( whence N = 10 )
     dz/dx = x.( y + z - t )         z(0) = 1
     dt/dx = x.y - z.t                  t(0) = 2
 
 

01 LBL "V"  02 R^  03 CHS  04 STO L  05 R^  06 ST* Y  07 ST+ L  08 CLX  09 RCL Z  10 R^  11 ST+ L  12 ST* Y  13 X<> Z  14 ST+ Y  15 ST* Z  16 X<> T  17 LASTX  18 *  19 X<>Y  20 RTN  21 END

 
  SF 03

  0     STO 01
  1     STO 02  STO 03
  2     STO 04
 0.1   STO 05
 10    STO 06    XEQ "RK4OP"   >>>>      x   =  1                                                       ---Execution time = 88s---
                                                     RDN     y(1) =  0.258209755
                                                     RDN     z(1) =  1.157620123
                                                     RDN     t(1) =  0.842178165

-The exact values are:

      y(1) = 0.258207906                      ( rounded to 9 decimals )
      z(1) = 1.157623981
      t(1) = 0.842178312

-Press R/S to continue with the same h and N or modify   R05-R06
 

Notes:

-When you press R/S after the 1st execution, registers R16 to R26 are unchanged,
  so the constants are still stored in these registers and lines 02 to 29 are no more executed.

-The results are often - but not always - more accurate than with the "classical" 4th-order Runge-Kutta method !
 
 

References:

[1]  J.C. Butcher  - "The Numerical Analysis of Ordinary Differential Equations" - John Wiley & Sons  -   ISBN  0-471-91046-5
[2]  Anthony Ralston - "Runge-Kutta Methods with Minimum Error Bounds"

-In reference [2], you can find coefficients for another optimal method