DM 41X: Doble integral using different modules fails

[rawi]

There are two modules for the HP41 that offer the possiblility for numeric integration of a function: The Advantage and the Math module. For each of the modules there is a warning that no function can be integrated that uses the integration routine of that module; so it is not possible to evaluate double integrals with one module. But there is no warning that the function for the Advantage-integration routine INTEG should not use the Math-Module integration routine INTG or vice veersa. So I came to the idea to test whether ist is possible to combine these two routines from the Advantage and the Math module as an easy way to evaluate double integrals.

I tested the evaluation of the double integral of the following function:
e^(-x^2/(y^2+y*x+1)) with limits for x and y from 0 to 1.
HP Prime delivers the result very fast: 0.822832

For the DM 41X I wrote two little programs:

LBL 'F1 @ for integration of the y variable
STO 00 @ y-Variable
0
ENTER
1 @ limits for y-variable
'F2
INTEG @ calls integragtion routine from Advantage module to integrate x-Variable
RTN
END

LBL 'F2 @ for the integration of the x variable
ENTER
ENTER @ x-value in stack
x^2 @ x^2
CHS
x<>y
RCL 00 @ y-value
* @ x*y
RCL 00
x^2
+
1
+ @ x*y + y^2 +1
/
E^X
RTN
END

Then I tried to evaluate the the double integral with the following commands:
XEQ 'INTG @ routine for Math module
0
ENTER
1 @ limits for y
f A
24 @ 24 numbers of interval used which is a lot
f B >> Calulator asks for function name >> 'F1 >> R/S
Result is 0.062505 which is far away from the real value.

Did I something wrong? Or simply shoudln't I combine those two routines from different modules? I could not find anything about that.

Thank you!

[Thomas Klemm]

From Math Pac:

Quote:26 Numerical Integration

The program uses Registers 00-07; the remaining registers are available for defining \(f_{i}(x)\).

Thus you have to use e.g. register 08 instead of 00:

Code:

01▸LBL "F1"     @ for integration of the y-variable
02 STO 08       @ y-variable
03 0
04 1            @ limits for y-variable
05 "F2"
06 INTEG        @ calls integration routine from Advantage module to integrate x-variable
07 END

Also this program can be made a bit shorter:

Code:

01▸LBL "F2"     @ for the integration of the x-variable
02 X↑2          @ x^2
03 CHS
04 LASTX        @ x-value
05 RCL 08       @ y-value
06 +            @ x+y
07 LASTX        @ y-value
08 *            @ (x+y)*y
09 1
10 +            @ x*y + y^2 + 1
11 /
12 E↑X
13 END

---

With FIX 4 I got after a while:

0.822832788

[rawi]

Thomas, thank you very much for your very helpful advice. Yes, this was the problem. Somehow I thought it is with the Advantage module which needs a certain number of registers but does not prohibit the use of certain registers. In the future I will read the instructions more carefully. Thank you as well for the tips to shorten the programs. I feel I could learn a lot from you.
Best
Raimund

[Thomas Klemm]

(09-30-2024 07:53 AM)rawi Wrote:  I feel I could learn a lot from you.

Aren't we all learning a lot in this great forum?
At least I am.

[rawi]

Thomas Klemm wrote:

Quote:Aren't we all learning a lot in this great forum?
At least I am.

Yes, so do I. That is why I nearly daily browsw this forum.

[Albert Chan]

(09-29-2024 08:08 PM)rawi Wrote:  I tested the evaluation of the double integral of the following function:
e^(-x^2/(y^2+y*x+1)) with limits for x and y from 0 to 1.
HP Prime delivers the result very fast: 0.822832

I was curious if there is a way to avoid double integral, and still get a good estimate.

For x=0..1, y=0..1 --> Denominator k = y^2+y*x+1 = 1..3

Assume k is just a constant. This reduced ∫∫ to ∫

∫(e^(-x^2/k), x=0..1) = √k * ∫e^(-t^2), t=0 .. 1/√k) = √(k*pi)/2 * erf(1/√k) = √(k*pi) * cdf(0 .. √(2/k))

I = ∫(e^(-x²/k), x=0..1) = √(2*pi) * cdf(0 .. z)/z      , where z = √(2/k)

Mean denominator
= ∫∫(y²+y*x+1), x=0..1, y=0..1)
= (∫(y²+1), y=0..1) + (∫x, x=0..1)*(∫y, y=0..1)
= (1/3 + 1) + 1/2*1/2
= 19/12

z = √(2/k) = √(24/19) = 1.1239

Since this is only rough estimate, I don't bother with interpolation.
We just look up from statistic book, cdf(0 .. z) where z = 1.12 and 1.13

I(z=1.12) = 2.5066 * 0.3686 / 1.12 = 0.8249
I(z=1.13) = 2.5066 * 0.3708 / 1.13 = 0.8225