# hp41programs

Dblinteg Double Integrals for the HP-41

We want to evaluate:   òab ò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:   /
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 )

 STACK INPUTS OUTPUTS X / I

Example:   Evaluate  I = ò12 òxx^2Ö(1+x4 y4) 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

n = 1     1  STO 03  XEQ "DIN"  >>>   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)