HP-41 Challenge: Double Integrals by INTEG Recursion

[gjmcclure]

Well, this is not the solution I want, but it will work. It is based on comments above on using a Monte Carlo approach. For the HP41, this approach takes a lot of time, unless you have a 50x speed simulator! I am not satisfied with the immense time vs. lack of precision it yields. I could convert it to MCODE, but I feer that would not be worth the effort.

Nevertheless, here is my program for N-dimension Monte Carlo integration, it requires the AMC OS/X module (E3/E+, SEEDT and RAND) and my GJMV2 module (X/E3, which simply divides X by 1000):

Code:


;
; MONTE CARLO INTEGRATION (N-TH ORDER)
;
; ORDER SHOULD BE IN R00
; NUMBER OF ITERATIONS (N) SHOULD BE IN R01
; POINTER TO RANDOM REGISTERS WILL BE SAVED IN R02 FOR USERS
; NAME OF USER FUNCTION SHOULD BE IN ALPHA, IT WILL BE SAVED IN R03
; INTEGRAL SUM WILL BE SAVED IN R04 (AS WILL BE FINAL RESULT)
; LIMITS IN R05-R2N+4
; RANDOM VALUES FOR USER FUNCTION IN R2N+5 TO R3N+4
;
LBL "MCINT"
CF 21        ; ALLOWS AVIEW TO SHOW STEP DOWN OF ITERATIONS
CLX
STO 04        ; CLEAR SUM
SEEDT        ; RANDOMIZE
ASTO 03        ; SAVE USER FUNCTION NAME
RCL 01        ; SAVE COUNT IN O
STO O
RCL 00        ; GET DIMENSION
ST+ X        ; DOUBLE FOR NUMBER OF REGISTERS FOR LIMITS
E3/E+        ; CONVERT TO ISG VALUE
4.004        ; BUMP TO POINT TO FIRST LL
+
STO M        ; SAVE IN M
RCL 00        ; CREATE POINTER TO RANDOM REGISTERS
X/E3
RCL 00
ST+ X
+
+
STO N        ; SAVE IN N
STO 02        ; FOR USER
;
; PRODUCE RANDOM VALUES
;
LBL 00
RCL IND M        ; LL(N) IN X
ENTER^        ; LL(N) IN X AND Y
ISG M
RCL IND M        ; UL(N) IN X, LL(N) IN Y AND Z
X<>Y
-            ; UL(N) - LL(N) IN X, LL(N) IN Y
RAND
*
+            ; RAND VALUE BETWEEN LL(N) AND UL(N) IN X
STO IND N
ISG M        ; CONTINUE
STO X
ISG N        ; BUMP RANDOM REGISTER POINTER
GTO 00        ; CONTINUE UNTIL DONE
RCL 00
ST- N            ; RESET LIMITS AND RANDOM VALUES COUNTERS
ST- M
ST- M
;
; CALL USER FUNCTION, SUM IN R04
;
XEQ IND 03
ST+ 04
VIEW O
DSE O
GTO 00
;
; FINAL RESULT
; CALC MULTIPLIER
;
1            ; INIT MULTIPLIER
LBL 01
RCL IND M        ; LL(N) IN X
ISG M
RCL IND M        ; UL(N) IN X, LL(N) IN Y
X<>Y
-            ; UL(N) - LL(N) IN X
*            ; NEW MULTIPLIER
ISG M
GTO 01
;
; MULTIPLY RESULT, DIVIDE BY NUMBER OF POINTS
;
ST* 04
RCL 02
ST/ 04
;
; RESTORE ALPHA NAME, DISPLAY RESULT
;
CLA
ARCL 03
RCL 04
CLD
END

; EXAMPLE OF QUINTUPLE INTEGRAL OF SQRT(6-X*X-Y*Y-Z*Z-U*U-V*V)
; X FROM 0 TO 0.7
; Y FROM 0 TO 0.8
; Z FROM 0 TO 0.9
; U FROM 0 TO 1.0
; V FROM 0 TO 1.1

; ALPHA = "5DINT"
; REGISTER 00 = 5
; REGISTER 01 = N (10, 100, AND 1000 USED FOR RUNS BELOW)
; REGISTERS 05, 07, 09, 11, 13 = 0
; REGISTER 06 = 0.7
; REGISTER 08 = 0.8
; REGISTER 10 = 0.9
; REGISTER 12 = 1.0
; REGISTER 14 = 1.1
LBL "5DINT"
6
RCL 15
X^2
RCL 16
X^2
RCL 17
X^2
+
+
-
RCL 18
X^2
RCL 19
X^2
+
-
SQRT
END

; 9 RUNS FOR N=10:   1.150 1.229 1.193 1.179 1.193 1.194 1.204 1.189 1.174
; 6 RUNS FOR N=100:  1.192 1.190 1.179 1.193 1.187 1.192
; 4 RUNS FOR N=1000: 1.186 1.192 1.192 1.190

You can see from the results that Monte Carlo integration on an HP41 is not ideal, but it does seem to work.

Greg.