nthorderODE

# Nth Order Differential Equations for the HP-41

Overview

1°) Second order Differential Equations
2°) Third order Differential Equations
3°) Nth order Differential Equations     ( X-Functions Module required )

-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 more easy to use, especially for the general case: less data registers are needed.
-All these routines use the "classical" fourth order Runge-Kutta method.

1°) Second Order Differential Equations

-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

2°) Third Order Differential Equations

-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

-We load, for instance:

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

3°) N-th Order Differential Equations

-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 "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

-We load, for instance:

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(4)(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