 The Museum of HP Calculators

# Nth Order Differential Equations for the HP-41

This program is Copyright © 2006 by Jean-Marc Baillard and is used here by permission.

This program is supplied without representation or warranty of any kind. Jean-Marc Baillard and The Museum of HP Calculators therefore assume no responsibility and shall have no liability, consequential or otherwise, of any kind arising from the use of this program material or any part thereof.

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

Z = y'
-In other words, the subroutine must change the stack from   Y = y   into   X = y" = f(x,y,y')
X = x

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

T = y"
Z = y'
-In other words, the subroutine must change the stack from   Y = y   into   X = y"' = f(x,y,y',y")
X = x

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

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

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                        If you don't have an HP-41CX, replace line 24 with  ENTER^  CLX  LBL 00  STO IND Y  ISG Y  GTO 00
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(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