+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Calculators (and very old HP Computers) (/forum-3.html)
+--- Forum: General Forum (/forum-4.html)
+--- Thread: HP-41 Challenge: Double Integrals by INTEG Recursion (/thread-6311.html)
Pages: 1 2
HP-41 Challenge: Double Integrals by INTEG Recursion - Ángel Martin - 05-27-2016 02:35 PM
Well folks, it's time to shake all those dormant brain cells and put them to a worthy task... if you'd agree.
Your mission is to write a program to calculate double integrals, i.e. those where the integrand is a function of two variables, say f(x,y), and where each of them is integrated along an interval - say [y1, y2] and [x1,x2] respectively.
Furthermore, allow for the possibility that the inner integral limits (x1, x2) could be a function of the outer variable (y).
And here is the key requirement: your program must use the INTEG function from the HP41-Advantage, and do it in a RECURSIVE manner - yes, what according to the manual is not possible.. oh well.
Input: the function name in ALPHA, and the four integration limits in the stack - plus a user program to define the function of course.
You can use as many data registers as you want. You can (and will need to) use functions from other modules, such as the AMC_OS/X and (big hint!) the RamPAGE...
I'll post my solution in a few days; it does all the work in just 31 program steps (not counting the function definition).
Six data registers are used, R00-R05. Can you beat that?
Happy recursion!
ÁM
RE: HP-41 Challenge: Double Integrals by INTEG Recursion - Ángel Martin - 05-30-2016 03:25 PM
Well, based on the overwhelming response so far (or lack thereof :-) I'll post an article with this subject in case somebody is interested at some point in time.
Cheers,
ÁM
RE: HP-41 Challenge: Double Integrals by INTEG Recursion - hth - 05-30-2016 04:33 PM
I have not calculated an integral for 35 years, and never even heard of the existence of a double.
Sounds like an interesting programming challenge, I assume you want to swap in and out the buffer of the Advantage module, maintaining more than instance of the buffer. I look forward to your article, maybe I will (as you say) find time to take a stab at it, at some point in the future.
Håkan
RE: HP-41 Challenge: Double Integrals by INTEG Recursion - Csaba Tizedes - 05-30-2016 06:41 PM
(05-30-2016 03:25 PM)Ángel Martin Wrote: Well, based on the overwhelming response so far (or lack thereof :-) I'll post an article with this subject in case somebody is interested at some point in time.
Cheers,
ÁM
OK, I'm interested. It is possible to discuss the method on 15C?!
Thanks!
Csaba
RE: HP-41 Challenge: Double Integrals by INTEG Recursion - Ángel Martin - 05-31-2016 05:03 AM
(05-30-2016 04:33 PM)hth Wrote: I have not calculated an integral for 35 years, and never even heard of the existence of a double.
I know what you mean, but it gets even stranger: rumor has it there are even triples!!
(05-30-2016 04:33 PM)hth Wrote: Sounds like an interesting programming challenge, I assume you want to swap in and out the buffer of the Advantage module, maintaining more than instance of the buffer. I look forward to your article, maybe I will (as you say) find time to take a stab at it, at some point in the future.
You're of course on the right track. The trick consists of changing the buffer-14 id# that is created by the first call to INTEG, so that when the second call happens it can create another buffer-14# and do its job on the second variable (X). Changing buffer id's is what function REIDBF does for a living, so there's a match made in heaven.
Timing and sync up are the only details to pay attention to. For instance there cannot be any key assignments (buffer-14# is placed BELOW those!) - The solution is simply to save them in XMEM at the beginning , clear them all, and restore them upon completion.
Cheers,
'AM
RE: HP-41 Challenge: Double Integrals by INTEG Recursion - Ángel Martin - 05-31-2016 05:06 AM
(05-30-2016 06:41 PM)Csaba Tizedes Wrote: OK, I'm interested. It is possible to discuss the method on 15C?!
afraid not, at least not using this trick. I don't know the insights of the 15C design but the "buffer" idea is surely not implemented in the same way, if used at all. It's more likely that the OS manages the memory directly using some other approach, not accessible to the user?
RE: HP-41 Challenge: Double Integrals by INTEG Recursion - Csaba Tizedes - 05-31-2016 10:02 AM
(05-31-2016 05:06 AM)Ángel Martin Wrote:(05-30-2016 06:41 PM)Csaba Tizedes Wrote: OK, I'm interested. It is possible to discuss the method on 15C?!
afraid not, at least not using this trick. I don't know the insights of the 15C design but the "buffer" idea is surely not implemented in the same way, if used at all. It's more likely that the OS manages the memory directly using some other approach, not accessible to the user?
OK, please post here one or two example what you want to show on HP-41 and I can see what kind of multiple integrals you want to calculate.
Thank you!
Csaba
RE: HP-41 Challenge: Double Integrals by INTEG Recursion - Ángel Martin - 05-31-2016 03:53 PM
(05-31-2016 10:02 AM)Csaba Tizedes Wrote: OK, please post here one or two example what you want to show on HP-41 and I can see what kind of multiple integrals you want to calculate.
The PPC article in V8N4p31 includes some simple examples - which I'm glad to say are also solved by the recursive approach.
It is embedded in the manual posted below:
Solve/Integrate ROM Manual
RE: HP-41 Challenge: Double Integrals by INTEG Recursion - gjmcclure - 05-31-2016 04:52 PM
Yes, there are triple, quadruple, even quintuple integral examples available on the net. There is no end to the order of integrals, as there is in theory no end to the number of dimensions one can mathematically describe (physically there may be, but then again string / loop theory seems to have come up with quite a few)...
I would like to see the integral program enhanced to do N-order integrals, and would suggest the following quintuple integral be used as a test (since it has been discussed at http://mathfaculty.fullerton.edu/mathews/n2003/montecarlo/MonteCarloMod/Links/MonteCarloMod_lnk_15.html and they suggest several methods for the solution...
With f(x,y,z,u,w) = sqrt(6-x^2-y^2-z^2-u^2-w^2) evaluate
integ(0,0.7) [ integ(0,0.8) [ integ(0,0.9) [ integ(0,1.0) [ integ(0,1.1) f(x,y,z,u,w) dw ]
du ] dz ] dy ] dx. The answer seems to lie around 1.189.
Greg
RE: HP-41 Challenge: Double Integrals by INTEG Recursion - Valentin Albillo - 06-01-2016 04:12 AM
(05-31-2016 04:52 PM)gjmcclure Wrote: I would like to see the integral program enhanced to do N-order integrals, and would suggest the following quintuple integral be used as a test (since it has been discussed at http://mathfaculty.fullerton.edu/mathews/n2003/montecarlo/MonteCarloMod/Links/MonteCarloMod_lnk_15.html and they suggest several methods for the solution...
With f(x,y,z,u,w) = sqrt(6-x^2-y^2-z^2-u^2-w^2) evaluate
integ(0,0.7) [ integ(0,0.8) [ integ(0,0.9) [ integ(0,1.0) [ integ(0,1.1) f(x,y,z,u,w) dw ]
du ] dz ] dy ] dx. The answer seems to lie around 1.189.
This quintuple integral is easy as pie for the HP-71B w/Math ROM using straight out-of-the-box code with no fancy programming or buffer juggling needed.
Assorted results for increasing precision (1E-1, 1E-2, ..., 1E-5) are as follows:
>LIST
10 DEF FNF(X,Y,Z,U,W)=SQR(6-X*X-Y*Y-Z*Z-U*U-W*W)
20 DEF FNG(X,Y,Z,U)=INTEGRAL(0,1.1,K,FNF(X,Y,Z,U,IVAR))
30 DEF FNH(X,Y,Z)=INTEGRAL(0,1,K,FNG(X,Y,Z,IVAR))
40 DEF FNI(X,Y)=INTEGRAL(0,.9,K,FNH(X,Y,IVAR))
50 DEF FNJ(X)=INTEGRAL(0,.8,K,FNI(X,IVAR))
60 FOR I=1 TO 5 @ K=1/10^I @ DISP K,INTEGRAL(0,.7,K,FNJ(IVAR)) @ NEXT I
>DESTROY ALL
>RUN
.1 1.18887862667
.01 1.18887862667
.001 1.18882510429
.0001 1.18878513051
.00001 1.18878333625
so we get from 5 to 8 correct digits give or take a couple ulps, as compared to Mathematica's 1.18878359.
Regards.
V.
.
RE: HP-41 Challenge: Double Integrals by INTEG Recursion - Ángel Martin - 06-01-2016 05:32 AM
(06-01-2016 04:12 AM)Valentin Albillo Wrote:(05-31-2016 04:52 PM)gjmcclure Wrote: With f(x,y,z,u,w) = sqrt(6-x^2-y^2-z^2-u^2-w^2) evaluate
integ(0,0.7) [ integ(0,0.8) [ integ(0,0.9) [ integ(0,1.0) [ integ(0,1.1) f(x,y,z,u,w) dw ]
du ] dz ] dy ] dx. The answer seems to lie around 1.189.
This quintuple integral is easy as pie for the HP-71B w/Math ROM using straight out-of-the-box code with no fancy programming or buffer juggling needed.
Glad this thread pulled you in from your greener pastures Valentín, always a pleasure to read your comments. Sure enough the 71B/MathPac is a vastly superior engine and the recursion functionality there is impressive - yet for a much humbler platform like the 41's the "buffer juggling" is a very elegant work-around - notwithstanding its inherent design limitations of course.
Saludos,
ÁM
RE: HP-41 Challenge: Double Integrals by INTEG Recursion - Ángel Martin - 06-01-2016 05:35 AM
FOCAL code shown below. Includes 6 program steps to preserve and restore the Key assignments so they'll be unmodified at the end of the calculations.
Code:
1 LBL "2DITG"
2 ASTO 05 function name
3 "KEYS"
4 SAVEKA sake keys
5 CLKEYS clear ka space
6 STO 04 x2
7 RDN
8 STO 03 x1
9 RDN
10 STO 02 y2
11 RDN
12 STO O1 y1
13 13
14 B? leftovers?
15 CLB yes, delete
16 RCL 01 y1
17 RCL 02 y2
18 "*2D" F(y)
19 INTEG outer integral
20 "KEYS"
21 GETKA restore keys
22 CLA
23 ARCL 05 function name
24 RTN final result in X.
25 LBL "*2D"
26 STO 00 y
27 14,013
28 REIDBF change buffer id#
29 RCL 03 x1
30 RCL 04 x2
31 CLA
32 ARCL 05 function name
33 INTEG inner integral
34 13,014
35 REIDBF change buffer id#
36 RDN solution to X
37 ENDRE: HP-41 Challenge: Double Integrals by INTEG Recursion - Valentin Albillo - 06-01-2016 06:51 PM
.
Hola, Ángel:
(06-01-2016 05:32 AM)Ángel Martin Wrote: Glad this thread pulled you in from your greener pastures Valentín, always a pleasure to read your comments.
Thanks for your interest & kind comment, Ángel. Matter of fact I do read the forum from time to time and still use my HP calculators (mostly HP-71B and HP-15C) and write code for them. Alas, sharing it is a different matter as I'm still looking for an adequate venue for my many unpublished articles, routines, tips, challenges and assorted stuff.
Quote:yet for a much humbler platform like the 41's the "buffer juggling" is a very elegant work-around
Indeed it is. I hope you didn't take my comment as somewhat derogatory or unappreciative, far from it, I was just stating the raw fact that the HP-71B w/Math ROM makes computing multiple integrals a cinch. I've seen your 41C w/enhancements solution to your own challenge and I think it's as elegant and short as possible, given the limitations you mention. Congratulations.
For the record, this is the HP-71B code to compute the quintuple integral using a Monte-Carlo approach, a straightforward 3-line affair with no Math ROM needed.
>LIST
10 DESTROY ALL @ RANDOMIZE 1 @ B=.7*.8*.9*1.1 @ FOR K=1 TO 4 @ N=10^K
20 S=0 @ FOR I=1 TO N @ X=RND*.7 @ Y=RND*.8 @ Z=RND*.9 @ U=RND @ W=RND*1.1
30 S=S+SQR(6-X*X-Y*Y-Z*Z-U*U-W*W) @ NEXT I @ DISP N,S*B/N @ NEXT K
>RUN
10 1.17631976176
100 1.19138525241
1000 1.18851295632
10000 1.18896583878
The output are the results for 10, 100, 1000 and 10000 random samples, the last having 5 correct digits (save one ulp) and about as fast as the INTEGRAL approach.
Best regards.
V.
RE: HP-41 Challenge: Double Integrals by INTEG Recursion - Ángel Martin - 06-02-2016 08:59 AM
(06-01-2016 06:51 PM)Valentin Albillo Wrote: For the record, this is the HP-71B code to compute the quintuple integral using a Monte-Carlo approach, a straightforward 3-line affair with no Math ROM needed.
>LIST
10 DESTROY ALL @ RANDOMIZE 1 @ B=.7*.8*.9*1.1 @ FOR K=1 TO 4 @ N=10^K
20 S=0 @ FOR I=1 TO N @ X=RND*.7 @ Y=RND*.8 @ Z=RND*.9 @ U=RND @ W=RND*1.1
30 S=S+SQR(6-X*X-Y*Y-Z*Z-U*U-W*W) @ NEXT I @ DISP N,S*B/N @ NEXT K
>RUN
10 1.17631976176
100 1.19138525241
1000 1.18851295632
10000 1.18896583878
The output are the results for 10, 100, 1000 and 10000 random samples, the last having 5 correct digits (save one ulp) and about as fast as the INTEGRAL approach.
Best regards.
V.
Clear and elegant; sometimes I wonder why I bother with RPN and FOCAL - so cumbersome in comparison...
Cheers,
ÁM
RE: HP-41 Challenge: Double Integrals by INTEG Recursion - Dieter - 06-02-2016 12:14 PM
(06-02-2016 08:59 AM)Ángel Martin Wrote: Clear and elegant; sometimes I wonder why I bother with RPN and FOCAL - so cumbersome in comparison...
Remember these JFK words:
"We choose to go to the moon in this decade and do the other things, not because they are easy, but because they are hard, because that goal will serve to organize and measure the best of our energies and skills, because that challenge is one that we are willing to accept, one we are unwilling to postpone, and one which we intend to win, and the others, too."
In a way, this sometimes also applies to FOCAL programming. ;-)
Dieter
RE: HP-41 Challenge: Double Integrals by INTEG Recursion - gjmcclure - 06-03-2016 03:02 AM
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.190You can see from the results that Monte Carlo integration on an HP41 is not ideal, but it does seem to work.
Greg.
RE: HP-41 Challenge: Double Integrals by INTEG Recursion - Namir - 06-03-2016 01:32 PM
Valentin,
Your use of Monte Carlo integration method is ingenious! The BASIC code can be translated to other BASIC dialects, calculator programming code, and even other PC programming languages!!
Hats off to you!!!
Namir
PS: I think I am adding your Monte Carlo solution to my list of tricks for HHC 2016.
RE: HP-41 Challenge: Double Integrals by INTEG Recursion - gjmcclure - 06-03-2016 02:46 PM
Agreed! Both Valentin's solutions are simple and elegant solutions! Actually, his 2nd solution is what lead to my HP41 solution! I just wish it was faster on the HP41.
Oh, I remember reading that some integrator solutions use what is called "Quasi-Monte Carlo Integration" (Mathematica is one). Anyone know a practical way to implement that? It has to do with what random numbers are chosen (random number sets in multiple dimentions don't always yield an even distribution sample)...
Greg
RE: HP-41 Challenge: Double Integrals by INTEG Recursion - Jake Schwartz - 06-03-2016 07:26 PM
[quote='Valentin Albillo' pid='56627' dateline='1464807065']
Alas, sharing it is a different matter as I'm still looking for an adequate venue for my many unpublished articles, routines, tips, challenges and assorted stuff.
[quote]
I don't mean to be controversial here, but is the HPCC Datafile still off limits?
Thanks,
Jake
RE: HP-41 Challenge: Double Integrals by INTEG Recursion - Valentin Albillo - 06-03-2016 10:22 PM
.
Hi, Namir:
(06-03-2016 01:32 PM)Namir Wrote: Valentin,
Your use of Monte Carlo integration method is ingenious! The BASIC code can be translated to other BASIC dialects, calculator programming code, and even other PC programming languages!!
Hats off to you!!!
Thanks for your continued appreciation of my humble efforts, Namir, I'm glad you like them and all the better if you find further uses for them.
Quote:PS: I think I am adding your Monte Carlo solution to my list of tricks for HHC 2016.
You're welcome to add it as you please, I'm sure your HHC 2016 talks will be as deservedly successful as they've been on previous years.
Best regards.
V.
.