The Museum of HP Calculators

Systems of 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.

Usage

The following program solves any system of n differential equations ( n<75 ) of the type:

dy₁/dx = f₁(x,y₁, ... ,y_n)
..................

dy_n/dx = f_n(x,y₁, ... ,y_n)

with the initial values: y_i(x₀) = y_i;0 i = 1 , ... , n

( For example, a second order equation like y''+xy'+y³=0 must be replaced by the equivalent system:

dy/dx = z
dz/dx = -xz-y³ )

It uses the "classical" fourth order Runge-Kutta method.
There are 2 instructions from the X-function module ( REGMOVE and REGSWAP )
and 1 from the HP-41 CX ( CLRGX at line 24 ).
If you don't have an HP-41 CX, line 24 can be replaced by the 5 lines:

    0
    LBL 00
    STO IND Y
    ISG Y
    GTO 00

Data Registers:

Registers R00 to R4n+10 are used and some of them must be initialized before executing "RK", namely:

R00 = the name of a subroutine which must calculate the n functions f_i and store the results in registers Rnn+11 to R2n+10
R01 = h/2 = half of the step size
R02 = N = the number of steps
R09 = n = the number of equations ( don't confuse it with the number of steps ... )

R10 = x₀
R11 = y_1;0 , ..... , Rnn+10 = y_n;0 ( the initial values )

Flags:

Flags F06 and F07 are used ( they are cleared at the end of the program )

Subroutine:

The subroutine called by "RK" must calculate the f_i ( x, y₁, ... , y_n ) using the values of x, y₁, ... , y_n in regiters R10, R11, ... , Rnn+10
and store f₁ ( x, y₁ , ... , y_n) , ... , f_n ( x, y₁ , ... , y_n) in Rnn+11, ... , R2n+10 respectively.

>> At the end of the execution of "RK" ,
         the new x is in registers X and R10
                     y₁(x) ----------Y and R11
                     y₂(x) ----------Z and R12
                     y₃(x) ----------T and R13 and so on ...

-The example below should make all these things clear.

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

(155 bytes / SIZE 4n+11)

Example:

Let's solve the system       dy₁/dx = - y₁y₂y₃
                                        dy₂/dx = x ( y₁+y₂-y₃ )
                                        dy₃/dx = xy₁- y₂y₃
with the initial values:   y₁(0) = 1 ; y₂(0) = 1 ; y₃(0) = 2

x ; y₁ ; y₂ ; y₃ will be in R10 ; R11 ; R12 ; R13 respectively, and we have to write a program to compute the 3 f_i
and store the results in R14 ; R15 ; R16 ( just after the "initial value registers" ).
For instance, the following subroutine will do the job ( let's call it "T" ):

01 LBL "T"
02 RCL 11
03 RCL 12
04 RCL 13
05 *
06 *
07 CHS
08 STO 14
09 RCL 11
10 RCL 12
11 +
12 RCL 13
13 -
14 RCL 10
15 *
16 STO 15
17 RCL 10
18 RCL 11
19 *
20 RCL 12
21 RCL 13
22 *
23 -
24 STO16
25 END ( 32 bytes )

We store the name of the subroutine in R00 "T" ASTO 00
and the numbers of equations in R09: 3 STO 09

then the initial values:      x    in R10                      0 STO 10
                                     y₁   in R11                      1 STO 11
                                     y₂   in R12                      1 STO 12
                                     y₃   in R13                      2 STO 13

Suppose we want to know y₁(1) ; y₂(1) ; y₃(1) with a step size h = 0.1 and therefore N = 10

h/2 is stored in R01 0.05 STO 01
N is stored in R02 10 STO 02 and XEQ "RK"

about 2mn24s later, the HP-41 gives:            x = 1                        in X and R10
                                                               y₁(1) = 0.2582093855 in Y and R11
                                                               y₂(1) = 1.157619553    in Z and R12
                                                               y₃(1) = 0.8421786516 in T and R13
To continue with the same h and N, simply press R/S and you'll have y₁(2) y₂(2) y₃(2)

-An estimation of the accuracy may be obtained by doing the calculation again with a smaller step size like h=0.05

It yields: y₁(1) = 0.2582079989 y₂(1) = 1.157623732 y₃(1) = 0.8421783424

Theoretically, when h tends to zero, the computed results tend to the exact solution,
but naturally, roundoff errors ( and execution time ) will limit the accuracy attainable!

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