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

-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
 

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

-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
 

Go back to the HP-41 software library
Go back to the general software library
Go back to the main exhibit hall