The Museum of HP Calculators

Double Integrals 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.

Introduction

We want to evaluate: ò_a^bò_u(x)^v(x) f(x;y) dx dy

The idea behind the following program is to use the 3 point Gauss-Legendre formula
in both x-axis and y-axis. There are at least 3 reasons for this choice:
1-The formula is relatively simple.
2-It has a quite good accuracy ( better than Simpson's rule )
3-It can be used even if the integrand has a singularity at the endpoints of the interval,
although in this case, convergence is of course much slower.

However, when the endpoints are infinite, it's necessary to make a change of argument ( like x = tan u )
to reduce the interval of integration to a bounded one.

Data Registers: Registers R00 thru R16 are used and some of them must be initialized before executing "DIN", namely:

         R00 = the name of the subroutine which calculates f(x;y)
         R01 = a
         R02 = b
         R03 = n = the number of subintervals into which the intervals (a;b) and (u(x);v(x)) are to be divided.
         R04 = the name of the subroutine which calculates u(x).
         R05 = the name of the subroutine which calculates v(x).

Flags: none.

Subroutines: 3 subroutines must be keyed into program memory:

      -The first one must take x and y from the X-register and Y-register respectively and calculate f(x;y).
      -The second takes x from the X-register and calculates u(x).
      -The third takes x from the X-register and calculates v(x).

N.B: When the program stops, the result is in X-register and in R06.

01  LBL "DIN"
02  1.6
03  STO 07
04  .15
05  SQRT
06  STO 08
07  RCL 02
08  RCL 01
09  STO 10
10  -
11  RCL 03
12  STO 15
13  /
14  STO 09
15  2
16  /
17  ST+ 10
18  CLX
19  STO 06
20  LBL 01
21  RCL 10
22  RCL 09
23  RCL 08
24  *
25  -
26  XEQ 02
27  ST+ 06
28  RCL 10
29  XEQ 02
30  RCL 07
31  *
32  ST+ 06
33  RCL 10
34  RCL 09
35  ST+ 10
36  RCL 08
37  *
38  +
39  XEQ 02
40  ST+ 06
41  DSE 15
42  GTO 01
43  GTO 04
44  LBL 02
45  STO 13
46  XEQ IND 04
47  STO 12
48  RCL 13
49  XEQ IND 05
50  RCL 12
51  -
52  RCL 03
53  STO 16
54  /
55  STO 11
56  2
57  /
58  ST+ 12
59  CLX
60  STO 14
61  LBL 03
62  RCL 12
63  RCL 11
64  RCL 08
65  *
66  -
67  RCL 13
68  XEQ IND 00
69  ST+ 14
70  RCL 12
71  RCL 13
72  XEQ IND 00
73  RCL 07
74  *
75  ST+ 14
76  RCL 12
77  RCL 11
78  ST+ 12
79  RCL 08
80  *
81  +
82  RCL 13
83  XEQ IND 00
84  ST+ 14
85  DSE 16
86  GTO 03
87  RCL 14
88  RCL 11
89  *
90  RTN
91  LBL 04
92  RCL 06
93  RCL 09
94  *
95  3.6
96  X^2
97  /
98  STO 06
99  END

( 140 bytes / SIZE 017 )

Example: Evaluate I = ò₁²ò_x^{x^2}Ö(1+x⁴ y⁴) dx dy.

We must load the 3 needed subroutines into program memory ( any global labels, maximum of 6 characters ),
for instance:

01 LBL "FF"
02 *
03 X^2
04 X^2
05 SIGN
06 LAST X
07 +
08 SQRT
09 RTN
10 LBL "U"
11 RTN
12 LBL "V"
13 X^2
14 RTN

and then, "FF" ASTO 00
                   1     STO 01
                   2     STO 02
                 "U" ASTO 04
                 "V" ASTO 05

let's try n = 1     1 STO 03 XEQ "DIN" the result is 15.45937082    in the X-register and in R06
             n = 2     2 STO 03      R/S            "      "     "   15.46673275
             n = 4     4 STO 03      R/S            "      "     "   15.46686031
             n = 8     8 STO 03      R/S            "      "     "   15.46686245 ... all the digits are correct!

With n = 8, execution time is about 7mn16s.
(When n is multiplied by 2, execution time is roughly multiplied by 4)

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