The Museum of HP Calculators

Nth Order Differential Equations for the HP-41

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(x₀) = y₀ , y'(x₀) = y'₀

Data Registers:           • R00 = function name                        ( Registers R00 thru R05 are to be initialized before executing "SRK" )
                                      • R01 = x₀
                                      • R02 = y₀        • R04 = h/2   where h = stepsize
                                      • R03 = y'₀       • 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'(x₀+m.h)

           Y             /      y(x₀+m.h)

           X             /        x₀+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' - y³ 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.x² + .... 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
      |                  *
      |                              *
      |                                           *
      |                                                              *
      |-----------------------------------------------------------*x₁-------- x
     0
              y(x₁) = 0 for x₁ = 6.896848619    and    y'(x₁) = -0.0424297576
and there is an inflexion point I with x_I = 1.495999168 , y_I = 0.720621687 and y'(x_I) = -0.279913175

-The solutions of the Lane-Emden Equations of index n ( y" + (2/x).y' + yⁿ = 0 , y(0)= 1 , y'(0) = 0 )
can be expressed by elementary functions for only 3 n-values:

   a) n = 0   y(x) = 1 - x²/6
   b) n = 1   y(x) = (sinx)/x
   c) n = 5   y(x) = ( 1+x²/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(x₀) = y₀ , y'(x₀) = y'₀ , y"(x₀) = y"₀

Data Registers:           • R00 = function name                           ( Registers R00 thru R06 are to be initialized before executing "TRK" )
                                      • R01 = x₀
                                      • R02 = y₀
                                      • R03 = y'₀       • R05 = h/2   where h = stepsize
                                      • R04 = y"₀      • 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"(x₀+m.h)

           Z             /      y'(x₀+m.h)

           Y             /      y(x₀+m.h)

           X             /        x₀+m.h

-Simply press R/S to continue with the same h and m

Example: y"' = 2xy" - x²y' + y² 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⁽ⁿ⁾ = f(x,y,y',y",.....,y^(n-1)) with the initial values: y(x₀) = y₀ , y'(x₀) = y'₀ , ........ , y^(n-1)(x₀) = y^(n-1)₀

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 = x₀ • R10 = y₀ • R11 = y'₀ • R12 = y"₀ ..................... • Rn+9 = y^(n-1)₀
Flags: /
Subroutine: A program which computes y⁽ⁿ⁾ = 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"(x₀+m.h)

           Z             /      y'(x₀+m.h)

           Y             /      y(x₀+m.h)

           X             /        x₀+m.h

-Simply press R/S to continue with the same h and m

Example: y⁽⁵⁾ = y⁽⁴⁾ - 2x.y"' + y" - y.y' with y(0) = 1 , y'(0) = y"'(0) = y⁽⁴⁾(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⁽⁴⁾(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⁽⁴⁾(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

STACK	INPUTS	OUTPUTS
T	/	y"(x₀+m.h)
Z	/	y'(x₀+m.h)
Y	/	y(x₀+m.h)
X	/	x₀+m.h