HP42s first major program (Double Integral) Best way to approach?

[DM48]

I am trying to convert a program from my HP48G to my HP42s and DM42.

Background
On the HP, I take three numbers from the stack and store them as local variables and have the equation /->NUM evaluate the expression to a number. It's really easy on the 48G, but slow (seconds to minutes depending on the numbers used. I use the iOS app m48+ to evaluate the expression fixed to five decimals. It computes instantly on my iPhone. I hoping to run this on the DM42 which is pocketable.

48G Stack

3:10
2:12
1:5

Stack Line 3 stores to a
Stack Line 2 stores to b
Stack Line 1 stores to c

Question
What is the best method to approach this double integral on the 42s? Keep in mind that a, b, and c are variables and I run this calc a lot when I run it. May run it over 100 times in one setting.

Example numbers are a (24), b (29), and c (12). The final answer would be multiplied by four giving an answer of 2325.62 for the example numbers.

\(\large I = \int_a^b\int_0^c\frac{r\sqrt{c^2-x^2}}{\sqrt{r^2-x^2}}\;\partial x\;\partial r \)

Any help getting me pointed in the correct direction would be appreciated. I am assuming I need to to solve the inner integral and then the outer integral on the 42s separately. This assumes the 42s doesn't have the same ability to evaluate the expression in the same way as the 48G. Correct?

[Thomas Okken]

The HP-42S doesn't have expressions, so you'd have to write a keystroke program instead, like this:

Code:

00 { 40-Byte Prgm }
01▸LBL "FN"
02 MVAR "R"
03 MVAR "C"
04 MVAR "X"
05 RCL "R"
06 RCL "C"
07 X↑2
08 RCL "X"
09 X↑2
10 -
11 SQRT
12 ×
13 RCL "R"
14 X↑2
15 RCL "X"
16 X↑2
17 -
18 SQRT
19 ÷
20 END

HOWEVER... The HP-42S and DM42 don't support double integrals. If you try to use INTEG inside a function that is itself being integrated, you'll get the error message Integ(Integ).

[Paul Dale]

Could the 15C do double integrals?
I don't have a 15C nearby to check...


Pauli

[Ángel Martin]

(05-27-2020 03:46 AM)Paul Dale Wrote:  Could the 15C do double integrals?
I don't have a 15C nearby to check...

I very much doubt it, neither does INTEG in the Advantage Pac - which is supposed to be a direct port from the 15C, or at least based on it. When you try it presents the RECURSION error message.

However recursive integrals are supported by FITG2 in the Advantage_Math ROM, which uses FINTG from the SandMath directly with a custom trick played on the INTEG buffer to allow dual existence while the inner integral is being calculated. Manual is available here, look for ADVG in the listing.

For this example the instructions are:

type "FN" in ALPHA
24, ENTER^, 29 ENTER^, 0, ENTER^, 12
XEQ "FITG2" with the desired accuracy in display mode.

With FIX 4 I get 581.4058 - which agrees with the given answer when multiplied by four.

the integrand function is as follows: inner integral is dx and outer integral is dr
On entry, the Inner variable is in the stack (four levels)
and the Outer variable is in R01

Code:

01 LBL "FN"
02 RCL 03
03 X^2
04 X<>Y  
05 X^2
06 -
07 SQRT
08 RCL 01
09 *
10 RCL 01
11 X^2
12 RCL Z
13 X^2
14  -
15 SQRT
16 /
17 END

[Werner]

In this case, we can convert it to a single integral by reversing the order of integration (how do I include these nice formulae? copy/paste of the LaTex code doesn't work, do I need to include it between codes?).
Then the inner integral becomes

int(a,b,r/(sqrt(r^2-x^2),dr)

which is simply sqrt(b^2-x^2) - sqrt(a^2-x^2), and we can use the following 42S program to evaluate the resulting integral:

Code:

00 { 47-Byte Prgm }
01▸LBL "FX"
02 MVAR "X"
03 MVAR "A"
04 MVAR "B"
05 RCL "B"
06 XEQ 00
07 RCL "A"
08 XEQ 00
09 -
10 RCL "ULIM"
11 XEQ 00
12 ×
13 RTN
14▸LBL 00
15 X↑2
16 RCL "X"
17 X↑2
18 -
19 SQRT
20 END

press shift-integrate, select FX, specify the values for A and B, press X, then put 0 for LLIM, C for ULIM, some small value for ACC (I used 0.05) and press the integrate symbol.

Cheers, Werner

[Paul Dale]

Ángel, the 15C manual talks about available registers for solve and integrate. I'm surprised that these aren't recursive.

The 34C was quite explicit that recursion wasn't possible, but that a SOLVE of an integral was. Perhaps vice-versa.

Pauli

[Thomas Okken]

The data structures used by SOLVE and INTEG are static, which is why those functions don't nest. They are separate from each other, so you can solve a function that uses INTEG or integrate a function that uses SOLVE.

I'm guessing HP didn't bother to make SOLVE and INTEG nestable because it would have been too slow to be useful anyway. In Free42 or the DM42, it would make more sense, but it would require additions to the instruction set since you would have to be able to explicitly distinguish between PGMSLV and PGMINT resetting the existing solver or integrator instance and creating new nested instances.

[J-F Garnier]

(05-27-2020 06:38 AM)Paul Dale Wrote:  Ángel, the 15C manual talks about available registers for solve and integrate. I'm surprised that these aren't recursive.

An attempt to integrate or solve recursively produces the Error 7 on the 15C, which is SOLVE(SOLVE) or INTEG(INTEG) as indicated on the back of the 15C :-)

The HP-71B Math ROM was the first to allow recursive solve or integrate functions, up to 5 levels. It was a major design challenge, as related in the HP Journal article.

J-F

[Ángel Martin]

(05-27-2020 07:15 AM)Thomas Okken Wrote:  I'm guessing HP didn't bother to make SOLVE and INTEG nestable because it would have been too slow to be useful anyway.

Indeed the convergence isn't optimized (zig-zagging as opposed to bee-lining) so you're right they are very slow at normal speed - but nowadays with V41, the 41X and TURBO 50 on the CL it's very decent performance. A bit of brute force can be ok some times...

[Ángel Martin]

(05-27-2020 07:16 AM)J-F Garnier Wrote:  The HP-71B Math ROM was the first to allow recursive solve or integrate functions, up to 5 levels. It was a major design challenge, as related in the HP Journal article.

Thanks for sharing, what an amazing design that was - and a great article that makes you long those yore times!

[J-F Garnier]

(05-27-2020 01:25 AM)DM48 Wrote:  Any help getting me pointed in the correct direction would be appreciated. I am assuming I need to to solve the inner integral and then the outer integral on the 42s separately.

One possible approach is to use the built-in integrate function for the inner integral, and code the Romberg algorithm in RPN for the outer integral. Speed penality should not be too large, since most of the time is spent in the inner integral.
The Romberg algorithm is well known, quite easy to implement and using low resources (the HP-42S has ample memory for it). The algorithm has been discussed recently here.

J-F

[rawi]

Just posted a program to compute double integrals for the HP 42s / DM 42:
https://www.hpmuseum.org/forum/thread-15...t=integral
Perhaps this helps?

Best

Raimund

[Albert Chan]

(05-27-2020 06:15 AM)Werner Wrote:  Then the inner integral becomes

int(a,b,r/(sqrt(r^2-x^2),dr)

which is simply sqrt(b^2-x^2) - sqrt(a^2-x^2) ...

Continued on, we have

\(I = \int_a^b\int_0^c\frac{r\sqrt{c^2-x^2}}{\sqrt{r^2-x^2}}\;dx\;dr
= \int _0 ^c \sqrt{c^2-x^2}\left( \sqrt{b^2-x^2} - \sqrt{a^2-x^2} \right)\;dx\)

The shape look like an ellipse, with semi-axis, c, c*(b-a)

a, b, c = 24, 29, 12
I ≈ pi/4 * (b-a) * c^2 ≈ 565.4866776461628

If integrand can be transformed to bell-shape, integral is easy to integrate.
see https://www.hpmuseum.org/forum/thread-13...#pid127590

Let \(x = c \cos(t)\; → dx = -c \sin(t)\;dt\)

\(I = \int _0 ^{\pi /2} (c^2-x^2) \left( \sqrt{b^2-x^2} - \sqrt{a^2-x^2}\right)\;dt \)

Shape of transformed curve: https://www.wolframalpha.com/input/?i=pl...0+to+pi%29

Code:

function I(a,b,c, n)    -- n = number of trapezoids
    local A, B = (a+c)*(a-c), (b+c)*(b-c)
    local function f(y) return y/(sqrt(B+y)+sqrt(A+y)) end
    local function g(t) t=c*sin(t); return f(t*t) end
    n = n or 8          -- default = 8 trapezoids
    local h, t = pi/(n+n), 0
    for i=1,n-1 do t = t + g(i*h) end
    return h * (b-a) * ((b+a)*t + c*c/2)
end

Convergence with trapezoid rule is amazingly fast:

lua> a, b, c = 24, 29, 12
lua> for n=1,8 do print(n, I(a,b,c, n)) end
1       565.4866776461628
2       581.3714770410786
3       581.40570463529
4       581.405804245352
5       581.4058045759407
6       581.4058045771336
7       581.4058045771384
8       581.4058045771384

[DM48]

Today was a great day. I woke up to so many replies from people who are waaaay stronger in Math and programming than myself. Thank you to all. I sincerely mean it.

My DM42 arrived today as well. This is great news too.

I want to give a little more background on the problem.

This is a 1/4 turn rotation and this is why the final answer is multiplied by four.

The function calculates the volume of the material removed from a hollow right cylinder when a hole is drilled into the side.

Inside diameter = a
Outside diameter = b
Cored hole = c

Wall thickness is the difference between b and a.

I have to study through these responses and will come back. Many thanks again to everyone. I am learning to program the 42s and am under the gun to move this program from my 48 to the 42s.

I am looking to make this the fastest calculation possible with at least five, preferably six decimals.

[DM48]

(05-27-2020 02:13 AM)Thomas Okken Wrote:  The HP-42S doesn't have expressions, so you'd have to write a keystroke program instead, like this:

Code:

00 { 40-Byte Prgm }
01▸LBL "FN"
02 MVAR "R"
03 MVAR "C"
04 MVAR "X"
05 RCL "R"
06 RCL "C"
07 X↑2
08 RCL "X"
09 X↑2
10 -
11 SQRT
12 ×
13 RCL "R"
14 X↑2
15 RCL "X"
16 X↑2
17 -
18 SQRT
19 ÷
20 END

HOWEVER... The HP-42S and DM42 don't support double integrals. If you try to use INTEG inside a function that is itself being integrated, you'll get the error message Integ(Integ).

I am trying to enter in this expression that Albert Chan supplied (tested this on my 48G and it works 33% faster than the double integral, thanks Albert!)

First Question

When entering the integral, how do you know what the calculator interprets to be under the square root house without using parentheses? I am use to the 48G Equation Writer to make sure I have crazy equations like this formatted correctly when entering in a program. The second part of Albert's expression has me puzzled how to key this in with certainty. Any documentation on this I can read or is it just a RPN Order of Operations thing?

\(I = \int _0 ^c \sqrt{c^2-x^2}\left( \sqrt{b^2-x^2} - \sqrt{a^2-x^2} \right)\;dx\)

[grsbanks]

(05-27-2020 09:46 PM)DM48 Wrote:  When entering the integral, how do you know what the calculator interprets to be under the square root house without using parentheses?

Becaue *you* decide when to press the \(\sqrt{x}\) key.

Pressing it immediately extracts the square root of what's in the X register. The DM42/HP-42S does not work with expressions like an RPL machine does.

[Werner]

Well, you could key in the program I supplied before, the post before Albert's.
Or you could use this, and use it interactively:

Code:

00 { 98-Byte Prgm }
01▸LBL "BORE"
02 MVAR "A"
03 MVAR "B"
04 MVAR "C"
05 1ᴇ-5
06 STO "ACC"
07 CLX
08 STO "LLIM"
09 PGMINT "FX"
10▸LBL 10
11 VARMENU "BORE"
12 STOP
13 RCL "C"
14 STO "ULIM"
15 INTEG "X"
16 4
17 ×
18 GTO 10
19▸LBL "FX"
20 RCL "B"
21 XEQ 01
22 RCL "A"
23 XEQ 01
24 -
25 RCL "C"
26 XEQ 01
27 ×
28 RTN
29▸LBL 01
30 X↑2
31 RCL "X"
32 X↑2
33 -
34 SQRT
35 END

To use, XEQ "BORE".
Enter the values for A, B and C (just enter a value and press the corresponding key), then press R/S. Reply is near instantaneous on a DM42.
You can enter new values - the previous ones are still there so you can just enter the ones that have changed, press R/S for the result.
Press EXIT to.. EXIT ;-)

It would help if you could transfer the file to your DM42, that would save you the trouble of keying it in. Just copy and paste the code into swissmicro's encoder page and upload the .raw file.

Cheers, Werner

[Werner]

Okay, what am I doing wrong?
I tried to derive the formula of the bore hole volume myself.
(Updated - first attempt was wrong)
The volume can be determined as a rotation integral of the area between the inner and outer diameters of the cylinder, and that is simply

\(V = 2\Pi\int _0^c r\left( \sqrt{b^2-r^2} - \sqrt{a^2-r^2} \right)\;dr\)

which, for a=24, b=29 and c=12 is 2393.10478..

Albert, help!

oops I know: it's not a rotational integral! Second attempt also wrong. OK.

[DM48]

It starts out as a triple integral. I’ll find my notes.

[Werner]

No need, I figured it out, but I get my simplified formula straight away, actually.
Werner