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

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 =
ò_{1}^{2
}ò_{x}^{x^2}Ö(1+x^{4}
y^{4}) 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