Threaded Mode | Linear Mode

Ángel Martin · (This post was last modified: 05-27-2016 05:33 PM by Ángel Martin.)

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

Ángel Martin · (This post was last modified: 05-30-2016 03:25 PM by Ángel Martin.)

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

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. Smile

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

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

Á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

Ángel Martin · (This post was last modified: 05-31-2016 05:06 AM by Ángel Martin.)

(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?

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

Ángel Martin · (This post was last modified: 05-31-2016 03:57 PM by Ángel Martin.)

(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

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...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

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...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.
.

Ángel Martin · (This post was last modified: 06-01-2016 05:32 AM by Ángel Martin.)

(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

Ángel Martin · (This post was last modified: 06-01-2016 05:38 AM by Ángel Martin.)

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

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.

Á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

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

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.190

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

Greg.

Namir · (This post was last modified: 06-03-2016 01:34 PM by Namir.)

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.

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

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

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.
.