HP-41 Challenge: Double Integrals by INTEG Recursion
05-27-2016, 02:35 PM (This post was last modified: 05-27-2016 05:33 PM by Ángel Martin.)
Post: #1
 Ángel Martin Senior Member Posts: 1,239 Joined: Dec 2013
HP-41 Challenge: Double Integrals by INTEG Recursion
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
05-30-2016, 03:25 PM (This post was last modified: 05-30-2016 03:25 PM by Ángel Martin.)
Post: #2
 Ángel Martin Senior Member Posts: 1,239 Joined: Dec 2013
RE: HP-41 Challenge: Double Integrals by INTEG Recursion
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
05-30-2016, 04:33 PM
Post: #3
 hth Senior Member Posts: 403 Joined: Mar 2014
RE: HP-41 Challenge: Double Integrals by INTEG Recursion
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
05-30-2016, 06:41 PM
Post: #4
 Csaba Tizedes Senior Member Posts: 489 Joined: May 2014
RE: HP-41 Challenge: Double Integrals by INTEG Recursion
(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
05-31-2016, 05:03 AM
Post: #5
 Ángel Martin Senior Member Posts: 1,239 Joined: Dec 2013
RE: HP-41 Challenge: Double Integrals by INTEG Recursion
(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
05-31-2016, 05:06 AM (This post was last modified: 05-31-2016 05:06 AM by Ángel Martin.)
Post: #6
 Ángel Martin Senior Member Posts: 1,239 Joined: Dec 2013
RE: HP-41 Challenge: Double Integrals by INTEG Recursion
(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?
05-31-2016, 10:02 AM
Post: #7
 Csaba Tizedes Senior Member Posts: 489 Joined: May 2014
RE: HP-41 Challenge: Double Integrals by INTEG Recursion
(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
05-31-2016, 03:53 PM (This post was last modified: 05-31-2016 03:57 PM by Ángel Martin.)
Post: #8
 Ángel Martin Senior Member Posts: 1,239 Joined: Dec 2013
RE: HP-41 Challenge: Double Integrals by INTEG Recursion
(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
05-31-2016, 04:52 PM
Post: #9
 gjmcclure Junior Member Posts: 28 Joined: Jul 2014
RE: HP-41 Challenge: Double Integrals by INTEG Recursion
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...nk_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
06-01-2016, 04:12 AM
Post: #10
 Valentin Albillo Senior Member Posts: 727 Joined: Feb 2015
RE: HP-41 Challenge: Double Integrals by INTEG Recursion
(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...nk_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.
.

All My Articles & other Materials here:  Valentin Albillo's HP Collection

06-01-2016, 05:32 AM (This post was last modified: 06-01-2016 05:32 AM by Ángel Martin.)
Post: #11
 Ángel Martin Senior Member Posts: 1,239 Joined: Dec 2013
RE: HP-41 Challenge: Double Integrals by INTEG Recursion
(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
06-01-2016, 05:35 AM (This post was last modified: 06-01-2016 05:38 AM by Ángel Martin.)
Post: #12
 Ángel Martin Senior Member Posts: 1,239 Joined: Dec 2013
RE: HP-41 Challenge: Double Integrals by INTEG Recursion
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    END
06-01-2016, 06:51 PM
Post: #13
 Valentin Albillo Senior Member Posts: 727 Joined: Feb 2015
RE: HP-41 Challenge: Double Integrals by INTEG Recursion
.
Hola, Ángel:

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.

All My Articles & other Materials here:  Valentin Albillo's HP Collection

06-02-2016, 08:59 AM
Post: #14
 Ángel Martin Senior Member Posts: 1,239 Joined: Dec 2013
RE: HP-41 Challenge: Double Integrals by INTEG Recursion
(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
06-02-2016, 12:14 PM
Post: #15
 Dieter Senior Member Posts: 2,397 Joined: Dec 2013
RE: HP-41 Challenge: Double Integrals by INTEG Recursion
(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
06-03-2016, 03:02 AM
Post: #16
 gjmcclure Junior Member Posts: 28 Joined: Jul 2014
RE: HP-41 Challenge: Double Integrals by INTEG Recursion
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.
06-03-2016, 01:32 PM (This post was last modified: 06-03-2016 01:34 PM by Namir.)
Post: #17
 Namir Senior Member Posts: 792 Joined: Dec 2013
RE: HP-41 Challenge: Double Integrals by INTEG Recursion
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.
06-03-2016, 02:46 PM
Post: #18
 gjmcclure Junior Member Posts: 28 Joined: Jul 2014
RE: HP-41 Challenge: Double Integrals by INTEG Recursion
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
06-03-2016, 07:26 PM
Post: #19
 Jake Schwartz Member Posts: 299 Joined: Dec 2013
RE: HP-41 Challenge: Double Integrals by INTEG Recursion
[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
06-03-2016, 10:22 PM
Post: #20
 Valentin Albillo Senior Member Posts: 727 Joined: Feb 2015
RE: HP-41 Challenge: Double Integrals by INTEG Recursion
.
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.
.

All My Articles & other Materials here:  Valentin Albillo's HP Collection

 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)