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 x² - 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(-x²)  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 x⁴ - 48 x² + 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(-x²)  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"