The Museum of HP Calculators

# Systems of Differential Equations for the HP-41

This program is Copyright © 1999 by Jean-Marc Baillard and may be freely used 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.

## Usage

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

dy1/dx = f1(x,y1, ... ,yn)
..................

dyn/dx = fn(x,y1, ... ,yn)

with the initial values:  yi(x0) = yi;0        i = 1 , ... , n

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

dy/dx = z
dz/dx = -xz-y3  )

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 fi 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 =  x0
R11 =  y1;0 , ..... ,  Rnn+10 = yn;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 fi ( x, y1, ... , yn )  using the values of  x, y1, ... , yn in regiters R10, R11, ... , Rnn+10
and store f1 ( x, y1 , ... , yn) , ... , fn ( x, y1 , ... , yn)  in Rnn+11, ... , R2n+10  respectively.

>>  At the end of the execution of "RK" ,
the new x is in registers X and R10
y1(x) ----------Y and R11
y2(x) ----------Z and R12
y3(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       dy1/dx = - y1y2y3
dy2/dx =  x ( y1+y2-y3 )
dy3/dx =  xy1- y2y3
with the initial values:   y1(0) = 1 ; y2(0) = 1 ; y3(0) = 2

x ; y1 ; y2 ; y3  will be in R10 ; R11 ; R12 ; R13 respectively, and we have to write a program to compute the 3 fi
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
y1   in R11                      1 STO 11
y2   in R12                      1 STO 12
y3   in R13                      2 STO 13

Suppose we want to know   y1(1) ; y2(1) ;  y3(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
y1(1) =  0.2582093855  in Y and R11
y2(1) =  1.157619553    in Z and R12
y3(1) =  0.8421786516  in T and R13
To continue with the same h and N, simply press R/S and you'll have y1(2) y2(2) y3(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:       y1(1) = 0.2582079989    y2(1) = 1.157623732    y3(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!