nthorderODE

# Nth Order Differential Equations for the HP-41

Overview

1°) Second order Differential Equations

a) Explicit Method of Order 4
b) Explicit Method of Order 6
c) Implicit Method of Order 6

2°) Third order Differential Equations

a) Explicit Method of Order 4
b) Explicit Method of Order 6

3°) Nth order Differential Equations

a) Explicit Method of Order 4     ( X-Functions Module required )
b) Explicit Method of Order 6

-Nth order differential equations can be re-written as a system of first order differential equations, so these programs may seem redundant !
-They are, however, much easier to use, especially for the general case.

Recently added:   §1b §1c §2b & §3b

1°) Second Order Differential Equations

a) Explicit Method of Order 4

-We want to solve 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 "SRK"  )

•  R01 = x0
•  R02 = y0        •  R04 = h/2   where h = stepsize
•  R03 = y'0       •  R05 =  m  = number of steps             R06 thru R09: 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

 01  LBL "SRK" 02  RCL 05 03  STO 08 04  LBL 01 05  RCL 03 06  RCL 02 07  RCL 01 08  XEQ IND 00 09  RCL 04 10  ST+ 01 11  * 12  STO 07 13  RCL 03 14  + 15  STO 09 16  LASTX 17  RCL 04 18  * 19  STO 06 20  RCL 02 21  + 22  RCL 01 23  XEQ IND 00 24  RCL 04 25  ST* 09 26  * 27  ST+ 07 28  ST+ 07 29  RCL 03 30  + 31  ENTER^ 32  X<> 09 33  ST+ 06 34  ST+ 06 35  RCL 02 36  + 37  RCL 01 38  XEQ IND 00 39  RCL 04 40  ST+ 01 41  ST+ X 42  ST* 09 43  * 44  ST+ 07 45  RCL 03 46  + 47  ENTER^ 48  X<> 09 49  ST+ 06 50  RCL 02 51  + 52  RCL 01 53  XEQ IND 00 54  RCL 04 55  ST* 09 56  * 57  RCL 07 58  + 59  3 60  / 61  ST+ 03 62  RCL 09 63  RCL 06 64  + 65  3 66  / 67  ST+ 02 68  DSE 08 69  GTO 01  70  RCL 03         71  RCL 02 72  RCL 01 73  END

( 103 bytes / SIZE 010 )

 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"      ST+ X      3              RTN             CHS                (  LBL "T" or any other global LBL , maximum 6 characters )
X=0?          X<>Y      Y^X         LBL 00         RTN
GTO 00      /               +              3
RCL Z        X<>Y      CHS         1/X

-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.05  STO 04   ( h/2 )
1   STO 02                            10    STO 05

XEQ "SRK"   >>>>       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/2 = 0.025  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

b) Explicit Method of Order 6

-"RKN6" employs a 7-stage 6th-order formula ( cf Runge-Kutta -> Runge-Kutta-Nystrom-General ODEs for the HP-41, example2 )

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

•  R01 = x          •  R04 = h                R06 to R12: temp
•  R02 = y          •  R05 = N
•  R03 = y'
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 327 is a three-byte  GTO 01

 01 LBL "RKN6"  02 RCL 05  03 STO 06  04 LBL 01  05 RCL 03  06 RCL 02  07 RCL 01  08 XEQ IND 00  09 RCL 04  10 *  11 STO 07  12 3  13 /  14 RCL 03  15 +  16 RCL 07  17 18  18 /  19 RCL 03  20 3  21 /  22 +  23 RCL 04  24 *  25 RCL 02  26 +  27 RCL 04  28 3  29 /  30 RCL 01  31 +  32 XEQ IND 00  33 RCL 04  34 *  35 STO 08  36 ST+ X  37 3  38 /  39 RCL 03  40 +  41 RCL 08  42 ST+ X  43 9  44 /  45 RCL 03  46 ST+ X  47 3 48 /  49 +  50 RCL 04  51 *  52 RCL 02  53 +  54 RCL 04  55 ST+ X  56 3  57 /  58 RCL 01  59 +  60 XEQ IND 00  61 RCL 04  62 *  63 STO 09  64 CHS  65 RCL 08  66 4  67 *  68 +  69 RCL 07  70 +  71 12  72 /  73 RCL 03  74 +  75 RCL 07  76 RCL 09  77 +  78 36  79 /  80 RCL 03  81 3  82 /  83 +  84 RCL 04  85 *  86 RCL 02  87 +  88 RCL 04  89 3  90 /  91 RCL 01  92 +  93 XEQ IND 00  94 RCL 04 95 *  96 STO 10  97 90  98 *  99 RCL 09 100 35 101 * 102 + 103 RCL 08 104 110 105 * 106 - 107 RCL 07         108 25 109 * 110 + 111 48 112 / 113 RCL 03 114 + 115 RCL 10 116 270 117 * 118 RCL 09 119 35 120 * 121 + 122 RCL 08 123 330 124 * 125 - 126 RCL 07 127 125 128 * 129 + 130 288 131 / 132 RCL 03 133 5 134 * 135 6 136 / 137 + 138 RCL 04 139 * 140 RCL 02 141 + 142 RCL 04 143 5 144 * 145 6 146 / 147 RCL 01 148 + 149 XEQ IND 00 150 RCL 04 151 * 152 STO 11 153 10 154 / 155 RCL 10 156 2 157 / 158 + 159 RCL 09 160 8 161 / 162 - 163 RCL 08 164 11 165 * 166 24 167 / 168 - 169 RCL 07 170 3 171 * 172 20 173 / 174 + 175 RCL 03 176 + 177 RCL 11 178 15 179 / 180 CHS 181 RCL 10 182 12 183 / 184 - 185 RCL 09 186 16 187 / 188 + 189 RCL 08 190 11 191 * 192 144 193 / 194 + 195 RCL 07 196 40 197 / 198 + 199 RCL 03         200 6 201 / 202 + 203 RCL 04 204 * 205 RCL 02 206 + 207 RCL 04 208 6 209 / 210 RCL 01 211 + 212 XEQ IND 00 213 RCL 04 214 * 215 STO 12 216 1600 217 * 218 RCL 11 219 128 220 * 221 + 222 RCL 10 223 2360 224 * 225 - 226 RCL 09 227 215 228 * 229 + 230 RCL 08 231 1980 232 * 233 + 234 RCL 07 235 783 236 * 237 - 238 780 239 / 240 RCL 03 241 + 242 RCL 12 243 4 E3 244 * 245 RCL 11 246 64 247 * 248 + 249 RCL 10         250 4720 251 * 252 - 253 RCL 09  254 215 255 * 256 + 257 RCL 08 258 3960 259 * 260 + 261 RCL 07 262 2349 263 * 264 - 265 2340 266 / 267 RCL 03 268 + 269 RCL 04 270 * 271 RCL 02 272 + 273 RCL 04 274 RCL 01 275 + 276 XEQ IND 00 277 RCL 04 278 * 279 RCL 12 280 80 281 * 282 RCL 11 283 16 284 * 285 + 286 RCL 10 287 110 288 * 289 + 290 RCL 09 291 55 292 * 293 + 294 RCL 07          295 39 296 * 297 + 298 600 299 / 300 RCL 03  301 + 302 RCL 04 303 ST+ 01 304 * 305 ST+ 02 306 CLX 307 RCL 07 308 + 309 13 310 * 311 RCL 09 312 RCL 10 313 + 314 55 315 * 316 + 317 RCL 11 318 RCL 12 319 + 320 32 321 * 322 + 323 200 324 / 325 ST+ 03 326 DSE 06 327 GTO 01 328 RCL 03 329 RCL 02 330 RCL 01 331 END

( 433 bytes / SIZE 013 )

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

Example:     y" = - 2 y - 2 x y'          y(0) = 1    y'(0) = 0         Calculate  y(1)

-Load for example this routine in main memory:

 01 LBL "T" 02 ST* Z 03 RDN 04 + 05 ST+ X 06 CHS 07 RTN

"T"  ASTO 00

0  STO 01  STO 03
1  STO 02

-With  h = 0.1  &  N = 10

0.1  STO 04  10  STO 05  XEQ "RKN6"   >>>>   x    =  1                     = R01
RDN  y(1) = 0.367879455   = R02
RDN y'(1) = -0.735758821  = R03

Notes:

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

c) Implicit Method of Order 6

"IRKN6" employs the implicit 3-stage 6th-order formula described in reference [2]

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

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

-Lines 215-252 are three-byte  GTO 01

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

( 322 bytes / SIZE 014 )

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

Example:     y" = - 2 y - 2 x y'          y(0) = 1    y'(0) = 0         Calculate  y(1)

-Load for example this routine in main memory:

 01 LBL "T" 02 ST* Z 03 RDN 04 + 05 ST+ X 06 CHS 07 RTN

"T"  ASTO 00

0  STO 01  STO 03
1  STO 02

-With  h = 0.1  &  N = 10

0.1  STO 04  10  STO 05  XEQ "IRKN6"   >>>>   x    =   1                     = R01
RDN  y(1) =  0.367879441   = R02
RDN y'(1) = -0.735758882   = R03

Notes:

-The HP41 displays the successive differences between 2 approximations of the implicit formula ( line 213 )
-Line 212 ( E-9 ) may lead to an infinite loop: change this number if need be.

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

-This formula is the best one listed on this page for second order ODEs, but of course not very fast with an HP41.
-The error for y(1) is only about 2 E-10 in the above example.

2°) Third Order Differential Equations

a) Explicit Method of Order 4

-Likewise, we want to solve 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 "TRK"  )

•  R01 = x0
•  R02 = y0
•  R03 = y'0       •  R05 = h/2   where h = stepsize
•  R04 = y"0      •  R06 =  m  = number of steps             R07 thru R12: 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

 01  LBL "TRK"   02  RCL 06   03  STO 10   04  LBL 16   05  RCL 04   06  RCL 03   07  RCL 02   08  RCL 01   09  XEQ IND 00   10  RCL 05   11  ST+ 01   12  *   13  STO 09   14  RCL 04   15  +   16  STO 11   17  RCL 05   18  LASTX   19  *   20  STO 08   21  RCL 03 22  +   23  STO 12   24  LASTX   25  RCL 05   26  *   27  STO 07   28  RCL 02   29  +   30  RCL 01   31  XEQ IND 00   32  RCL 05   33  ST* 11   34  ST* 12   35  *   36  ST+ 09   37  ST+ 09   38  RCL 04   39  +   40  ENTER^   41  X<> 11   42  ST+ 08 43  ST+ 08   44  RCL 03   45  +   46  ENTER^   47  X<> 12   48  ST+ 07   49  ST+ 07   50  RCL 02   51  +   52  RCL 01   53  XEQ IND 00   54  RCL 05   55  ST+ 01   56  ST+ X   57  ST* 11   58  ST* 12   59  *   60  ST+ 09   61  RCL 04   62  +   63  ENTER^ 64  X<> 11   65  ST+ 08   66  RCL 03   67  +   68  ENTER^   69  X<> 12   70  ST+ 07   71  RCL 02   72  +   73  RCL 01   74  XEQ IND 00   75  RCL 05   76  ST* 11   77  ST* 12   78  *   79  RCL 09   80  +   81  3   82  /   83  ST+ 04   84  RCL 11 85  RCL 08   86  +   87  3   88  /   89  ST+ 03   90  RCL 12   91  RCL 07    92  +   93  3   94  /   95  ST+ 02   96  DSE 10   97  GTO 16   98  RCL 04          99  RCL 03  100  RCL 02 101  RCL 01 102  END

( 143 bytes / SIZE 013 )

 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

LBL "S"     ST+ X     -                         (  LBL "S" or any other global LBL , maximum 6 characters )
X^2           ST* T      -
ST* Z        RDN       RTN
X<> L       X^2

-With  h =  0.1  and  m = 10

"S"  ASTO 00
0   STO 01  STO 03                  0.05  STO 05   ( h/2 )
1   STO 02  CHS  STO 04         10    STO 06

XEQ "TRK"   >>>>       x   =  1                            ( in 49 seconds )
RDN     y(1)  =  0.595434736
RDN     y'(1) = -0.776441445
RDN    y"(1) =  -0.791715205

-With  h/2 = 0.025  and  m = 20,  we would find  y(1) = 0.595431304  ,  y'(1) = -0.776444326  ,  y"(1) = -0.791718300

b) Explicit Method of Order 6

-We can get a better precision with a higher order Runge-Kutta formula:

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

•  R01 = x0
•  R02 = y0
•  R03 = y'0       •  R05 =  h  = stepsize
•  R04 = y"0      •  R06 =  N = number of steps             R07 thru R25: 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

-Line 433 is a three-byte GTO 01

 01 LBL "TRK6"  02 RCL 06  03 STO 25  04 LBL 01  05 RCL 04  06 STO 08  07 RCL 03  08 STO 07  09 RCL 02  10 RCL 01  11 XEQ IND 00  12 RCL 05  13 ST* 07  14 ST* 08  15 *  16 STO 09  17 3  18 /  19 RCL 04  20 +  21 STO 11  22 RCL 08  23 3  24 /  25 RCL 03  26 +  27 STO 10  28 RCL 07  29 3  30 /  31 RCL 02  32 +  33 RCL 05  34 3  35 /  36 RCL 01  37 +  38 RCL 11  39 RDN  40 XEQ IND 00  41 RCL 05  42 ST* 10  43 ST* 11  44 *  45 STO 12  46 1.5  47 /  48 RCL 04  49 +  50 STO 14  51 RCL 11  52 1.5  53 /  54 RCL 03  55 +  56 STO 13  57 RCL 10  58 1.5  59 /  60 RCL 02  61 +  62 RCL 05  63 1.5 64 /  65 RCL 01  66 +  67 RCL 14  68 RDN  69 XEQ IND 00  70 RCL 05  71 ST* 13  72 ST* 14  73 *  74 STO 15  75 CHS  76 RCL 12  77 4  78 *  79 +  80 RCL 09  81 +  82 12  83 /  84 RCL 04  85 +  86 STO 17  87 RCL 11  88 4  89 *  90 RCL 14  91 -  92 RCL 08  93 +  94 12  95 /  96 RCL 03  97 +  98 STO 16  99 RCL 10 100 4 101 * 102 RCL 13 103 - 104 RCL 07 105 + 106 12 107 / 108 RCL 02 109 + 110 RCL 05 111 3 112 / 113 RCL 01 114 + 115 RCL 17 116 RDN 117 XEQ IND 00 118 RCL 05 119 ST* 16 120 ST* 17 121 * 122 STO 18 123 90 124 * 125 RCL 15 126 35 127 * 128 + 129 RCL 12 130 110 131 * 132 - 133 RCL 09 134 25 135 * 136 + 137 48 138 / 139 RCL 04  140 + 141 STO 20         142 RCL 17 143 90 144 * 145 RCL 14 146 35 147 * 148 + 149 RCL 11 150 110 151 * 152 - 153 RCL 08 154 25 155 * 156 + 157 48 158 / 159 RCL 03 160 + 161 STO 19 162 RCL 16 163 90 164 * 165 RCL 13 166 35 167 * 168 + 169 RCL 10 170 110 171 * 172 - 173 RCL 07 174 25 175 * 176 + 177 48 178 / 179 RCL 02 180 + 181 RCL 05 182 1.2 183 / 184 RCL 01 185 + 186 RCL 20 187 RDN 188 XEQ IND 00 189 RCL 05 190 ST* 19 191 ST* 20 192 * 193 STO 21 194 10 195 / 196 RCL 18 197 2 198 / 199 + 200 RCL 15 201 8 202 / 203 - 204 RCL 12         205 11 206 * 207 24 208 / 209 - 210 RCL 09 211 .15 212 * 213 + 214 RCL 04 215 + 216 STO 23 217 RCL 20 218 10 219 / 220 RCL 17 221 2 222 / 223 + 224 RCL 14 225 8 226 / 227 - 228 RCL 11 229 11 230 * 231 24 232 / 233 - 234 RCL 08 235 .15 236 * 237 + 238 RCL 03  239 + 240 STO 22 241 RCL 19 242 10 243 / 244 RCL 16 245 2 246 / 247 + 248 RCL 13 249 8 250 / 251 - 252 RCL 10 253 11 254 * 255 24 256 / 257 - 258 RCL 07 259 .15 260 * 261 + 262 RCL 02 263 + 264 RCL 05 265 6 266 / 267 RCL 01 268 + 269 RCL 23 270 RDN 271 XEQ IND 00 272 RCL 05 273 ST* 22 274 ST* 23 275 * 276 STO 24 277 40 278 X^2 279 * 280 RCL 21 281 128 282 * 283 + 284 RCL 18 285 2360 286 * 287 - 288 RCL 15 289 215 290 * 291 + 292 RCL 12 293 1980 294 * 295 + 296 RCL 09 297 783 298 * 299 - 300 780 301 / 302 RCL 04 303 + 304 X<> 11 305 1980 306 * 307 RCL 23 308 40 309 X^2 310 * 311 + 312 RCL 20 313 128 314 * 315 + 316 RCL 17 317 2360 318 * 319 - 320 RCL 14 321 215 322 * 323 + 324 RCL 08 325 783 326 * 327 - 328 780 329 / 330 RCL 03         331 + 332 ENTER 333 X<> 10 334 1980 335 * 336 RCL 22  337 40 338 X^2 339 * 340 + 341 RCL 19 342 128 343 * 344 + 345 RCL 16 346 2360 347 * 348 - 349 RCL 13 350 215 351 * 352 + 353 RCL 07 354 783 355 * 356 - 357 780 358 / 359 RCL 02 360 + 361 RCL 01 362 RCL 05 363 + 364 RCL 11 365 RDN 366 XEQ IND 00 367 RCL 05 368 ST+ 01 369 ST* 10 370 ST* 11 371 * 372 STO 12 373 RCL 09 374 + 375 13 376 * 377 RCL 15 378 RCL 18 379 + 380 55 381 * 382 + 383 RCL 21 384 RCL 24 385 + 386 32 387 * 388 + 389 .5 390 % 391 ST+ 04 392 RCL 08          393 RCL 11 394 + 395 13 396 * 397 RCL 14  398 RCL 17 399 + 400 55 401 * 402 + 403 RCL 20 404 RCL 23 405 + 406 32 407 * 408 + 409 .5 410 % 411 ST+ 03 412 RCL 07 413 RCL 10 414 + 415 13 416 * 417 RCL 13 418 RCL 16 419 + 420 55 421 * 422 + 423 RCL 19 424 RCL 22  425 + 426 32 427 * 428 + 429 .5 430 % 431 ST+ 02 432 DSE 25 433 GTO 01 434 RCL 04 435 RCL 03 436 RCL 02 437 RCL 01 438 END

( 617 bytes / SIZE 026 )

 STACK INPUTS OUTPUTS T / y"(x0+N.h) Z / y'(x0+N.h) Y / y(x0+N.h) X / x0+N.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

LBL "S"     ST+ X     -
X^2           ST* T      -
ST* Z        RDN       RTN
X<> L       X^2

-With  h =  0.1  and  N = 10

"S"  ASTO 00
0   STO 01  STO 03                   0.1   STO 05
1   STO 02  CHS  STO 04         10    STO 06

XEQ "TRK6"   >>>>       x   =  1
RDN     y(1)  =  0.595431073
RDN     y'(1) = -0.776444516
RDN    y"(1)  = -0.791718501

-With  h = 0.05  &  N = 20,  it yields  y(1) = 0.5954310715  ,  y'(1) = -0.7764445234  ,  y"(1) = -0.7917185205

Note:

-The exact results, rounded to 10 decimals are:

y(1) = 0.5954310718  ,  y'(1) = -0.7764445229  ,  y"(1) = -0.7917185202

3°) N-th Order Differential Equations

a) Explicit Method of Order 4

-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
-"NRK" employs the "classical" fourth-order Runge-Kutta formula.

Data Registers:           •  R00 = function name          ( Registers R00 thru R03 and R09 thru Rn+9 are to be initialized before executing "NRK"  )

•  R01 = n      ( order )
•  R02 = h/2    where 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 .... Rn+9

-If you don't have an HP-41CX, replace line 24 with  ENTER^  CLX  LBL 00  STO IND Y  ISG Y  GTO 00

 01  LBL "NRK"  02  RCL 03 03  STO 08 04  RCL 01 05  9.009 06  + 07  STO 04 08  INT 09  RCL 01 10  + 11  STO 05 12  LASTX 13  + 14  STO 06 15  RCL 04 16  INT 17  RCL 05 18  1 19  ST+ Z 20  + 21   E3 22  / 23  + 24  CLRGX  25  10 26  LASTX 27  + 28  RCL 01 29   E6 30  / 31  + 32  STO 07 33  GTO 03 34  LBL 01 35  XEQ IND 00 36  LBL 02 37  RCL 02 38  * 39  ST+ IND 05 40  FS? 05 41  ST+ IND 05 42  FC? 06 43  ST+ X 44  RCL IND 06 45  + 46  X<> IND 04 47  DSE 06 48  DSE 05 49  DSE 04 50  GTO 02 51  RCL 01 52  ST+ 04 53  ST+ 05 54  ST+ 06 55  RTN 56  LBL 03 57  RCL 07 58  REGMOVE  59  CF 05 60  SF 06 61  XEQ 01 62  SF 05 63  RCL 02 64  ST+ 09 65  XEQ 01 66  CF 06 67  XEQ 01 68  CF 05 69  RCL 02 70  ST+ 09 71  XEQ 01 72  RCL 07 73  REGSWAP  74  RCL 04 75  RCL 05 76  3 77  SIGN 78  LBL 04 79  CLX 80  X<> IND Y  81  LASTX 82  / 83  ST+ IND Z  84  DSE Y 85  DSE Z 86  GTO 04 87  DSE 08 88  GTO 03 89  RCL 12         90  RCL 11 91  RCL 10 92  RCL 09 93  END

( 150 bytes / SIZE 3n+10 )

 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

LBL "W"       RCL 13        +                  -                        (  LBL "W" or any other global LBL , maximum 6 characters )
RCL 14         *                  RCL 10        RTN
RCL 09         -                  RCL 11
ST+ X           RCL 12       *

-With  h =  0.1  and  m = 10

"W"   ASTO 00                                                       and the initial values:     0  STO 09   STO 11  STO 13  STO 14
5    STO 01        ( fifth order equation )                                                    1  STO 10  CHS  STO 12
0.05  STO 02        ( h/2 = 0.05 )
10    STO 03        ( m = 10 )

XEQ "NRK"   >>>>       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""(1) = -0.897595630

-With  h/2 = 0.025  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

b) Explicit Method of Order 6

Data Registers:           •  R00 = function name          ( Registers R00 to R03 and R19 to Rn+19 are to be initialized before executing "NRK6"  )

•  R01 =  n      ( order )
•  R02 =  h  =  stepsize                       R04 thru R14 & Rn+20 thru R8n+19: temp               R15-R16-R17-R18:  unused
•  R03 =  m  =  number of steps

•  R19 = x0     •  R20 = y0     •  R21 = y'0    •  R22 = y"0   .....................  •  Rn+19 = 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(n-1)  are in registers  R19 R20 R21 .... Rn+19

-Line 367 is a three-byte GTO 10

 01 LBL "NRK6"  02 RCL 03  03 STO 12  04 RCL 01  05 RCL 01  06 RCL 01  07 STO 13  08 19.019  09 +  10 STO 04  11 +  12 STO 05  13 +  14 STO 06  15 +  16 STO 07  17 +  18 STO 08  19 +  20 STO 09  21 +  22 STO 10  23 +  24 STO 11  25 LBL 10  26 RCL 19  27 STO 14  28 LBL 01  29 RCL IND 04  30 STO IND 05  31 RCL 02  32 *  33 DSE 06  34 DSE 13  35 STO IND 06  36 DSE 05  37 DSE 04  38 GTO 01  39 RCL 01  40 ST+ 04  41 ST+ 05  42 ST+ 06  43 STO 13  44 XEQ IND 00  45 RCL 02  46 *  47 STO IND 06  48 LASTX  49 3  50 /  51 RCL 14  52 +  53 STO 19 54 LBL 02  55 RCL IND 06  56 3  57 /  58 RCL IND 05  59 +  60 STO IND 04  61 RCL 02  62 *  63 DSE 07  64 DSE 13  65 STO IND 07  66 DSE 06  67 DSE 05  68 DSE 04  69 GTO 02  70 RCL 01  71 ST+ 04  72 ST+ 05  73 ST+ 06  74 ST+ 07  75 STO 13  76 XEQ IND 00  77 RCL 02  78 *  79 STO IND 07  80 RCL 14  81 RCL 02  82 1.5  83 /  84 +  85 STO 19  86 LBL 03  87 RCL IND 07  88 1.5  89 /  90 RCL IND 05  91 +  92 STO IND 04  93 RCL 02  94 *  95 DSE 08  96 DSE 13  97 STO IND 08  98 DSE 07  99 DSE 05 100 DSE 04 101 GTO 03 102 RCL 01 103 ST+ 04 104 ST+ 05 105 ST+ 07 106 ST+ 08 107 STO 13 108 XEQ IND 00 109 RCL 02 110 * 111 STO IND 08 112 RCL 14 113 RCL 02 114 3 115 / 116 + 117 STO 19 118 LBL 04 119 RCL IND 07 120 4 121 * 122 RCL IND 08 123 - 124 RCL IND 06 125 + 126 12 127 / 128 RCL IND 05 129 + 130 STO IND 04 131 RCL 02 132 * 133 DSE 09 134 DSE 13 135 STO IND 09 136 DSE 08 137 DSE 07 138 DSE 06 139 DSE 05 140 DSE 04 141 GTO 04 142 RCL 01 143 ST+ 04 144 ST+ 05 145 ST+ 06 146 ST+ 07 147 ST+ 08 148 ST+ 09 149 STO 13 150 XEQ IND 00 151 RCL 02 152 * 153 STO IND 09 154 RCL 14 155 RCL 02 156 1.2 157 / 158 + 159 STO 19 160 LBL 05 161 RCL IND 09 162 90 163 * 164 RCL IND 08 165 35 166 * 167 + 168 RCL IND 07 169 110 170 * 171 - 172 RCL IND 06 173 25 174 * 175 + 176 48 177 / 178 RCL IND 05 179 + 180 STO IND 04 181 RCL 02 182 * 183 DSE 10 184 DSE 13 185 STO IND 10 186 DSE 09 187 DSE 08 188 DSE 07 189 DSE 06 190 DSE 05 191 DSE 04 192 GTO 05 193 RCL 01 194 ST+ 04 195 ST+ 05 196 ST+ 06 197 ST+ 07 198 ST+ 08 199 ST+ 09 200 ST+ 10 201 STO 13 202 XEQ IND 00 203 RCL 02 204 * 205 STO IND 10 206 RCL 14 207 RCL 02 208 6 209 / 210 + 211 STO 19 212 LBL 06 213 RCL IND 10 214 10 215 / 216 RCL IND 09 217 2 218 / 219 + 220 RCL IND 08 221 8 222 / 223 - 224 RCL IND 07 225 11 226 * 227 24 228 / 229 - 230 RCL IND 06 231 .15 232 * 233 + 234 RCL IND 05 235 + 236 STO IND 04 237 RCL 02 238 * 239 DSE 11 240 DSE 13 241 STO IND 11 242 DSE 10 243 DSE 09 244 DSE 08 245 DSE 07 246 DSE 06 247 DSE 05 248 DSE 04 249 GTO 06 250 RCL 01 251 ST+ 04 252 ST+ 05 253 ST+ 06 254 ST+ 07 255 ST+ 08 256 ST+ 09 257 ST+ 10 258 ST+ 11 259 STO 13 260 XEQ IND 00 261 RCL 02 262 * 263 STO IND 11 264 RCL 14 265 RCL 02 266 + 267 STO 19 268 RCL IND 07 269 LBL 07 270 1980 271 * 272 RCL IND 11 273 40 274 X^2 275 * 276 + 277 RCL IND 10 278 128 279 * 280 + 281 RCL IND 09 282 2360 283 * 284 - 285 RCL IND 08 286 215 287 * 288 + 289 RCL IND 06 290 783 291 * 292 - 293 780 294 / 295 RCL IND 05 296 + 297 STO IND 04 298 RCL 02 299 * 300 DSE 07 301 DSE 13 302 X<> IND 07 303 DSE 11 304 DSE 10 305 DSE 09 306 DSE 08 307 DSE 06 308 DSE 05 309 DSE 04 310 GTO 07 311 RCL 01 312 ST+ 04 313 ST+ 05 314 ST+ 06 315 ST+ 07 316 ST+ 08 317 ST+ 09 318 ST+ 10 319 ST+ 11 320 STO 13 321 XEQ IND 00 322 RCL 02 323 * 324 STO IND 07 325 LBL 08 326 RCL IND 06 327 RCL IND 07 328 + 329 13 330 * 331 RCL IND 08 332 RCL IND 09 333 + 334 55 335 * 336 + 337 RCL IND 10 338 RCL IND 11 339 + 340 32 341 * 342 + 343 .5 344 % 345 RCL IND 05 346 + 347 STO IND 04 348 DSE 11 349 DSE 10 350 DSE 09 351 DSE 08 352 DSE 07 353 DSE 06 354 DSE 05 355 DSE 04 356 GTO 08 357 RCL 01 358 ST+ 04 359 ST+ 05 360 ST+ 06 361 ST+ 07 362 ST+ 08 363 ST+ 09 364 ST+ 10 365 ST+ 11 366 DSE 12 367 GTO 10 368 RCL 22 369 RCL 21 370 RCL 20 371 RCL 19 372 END

( 611 bytes / SIZE 8n+20 )

 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

LBL "W"       RCL 23        +                  -                        (  LBL "W" or any other global LBL , maximum 6 characters )
RCL 24         *                  RCL 20        RTN
RCL 19         -                  RCL 21
ST+ X           RCL 22       *

-With  h =  0.1  and  m = 10

"W"   ASTO 00

5    STO 01        ( fifth order equation )                 and the initial values:     0  STO 19   STO 21  STO 23  STO 24
0.1   STO 02        ( h = 0.1 )                                                                      1  STO 20  CHS  STO 22
10    STO 03        ( m = 10 )

SIZE 060 ( or larger )

XEQ "NRK6"   >>>>       x   =  1                      = R19
RDN     y(1)  =  0.491724179   = R20
RDN     y'(1) = -1.041200379   = R21       and     RCL 23 =  y"'(1)  = -0.479804016
RDN    y"(1) =  -1.163353546  = R22                  RCL 24 = y""(1)  = -0.897594395

-With  h = 0.05  and  m = 20,  it yields:

y(1) = 0.491724180 ,  y'(1) = -1.041200380 ,  y"(1) = -1.163353546 ,  y'''(1)  = -0.479804017 ,  y""(1) = -0.897594397

References:

[1]  J.C. Butcher  - "The Numerical Analysis of Ordinary Differential Equations" - John Wiley & Sons  -   ISBN  0-471-91046-5
[2]  Sa Agam, Ya Yahaya, Sc Osuala - "An Implicit Runge-Kutta Method for General Second Order Ordinary Differential Equations"