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

[DM48]

(05-28-2020 11:11 AM)Werner Wrote:  No need, I figured it out, but I get my simplified formula straight away, actually.
Werner

When you say simplified formula, do you means the one Albert posted?

Can you get this to the final stage where the solver isn't needed?

Thank you for the code! I can't over state that. Loading it up now.

I am going to write it out on paper and work this on the calculator for practice. I suspect RPN users have a clear advantage when using variables in calcs compared to us RPL users.

[Werner]

This one, not Albert's last one:

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

Which I posted before, but not in pretty print; figured out how to do that in the meantime.
This is the integral solved by the BORE program;
Werner

[DM48]

(05-28-2020 12:47 PM)Werner Wrote:  This one, not Albert's last one:

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

Which I posted before, but not in pretty print; figured out how to do that in the meantime.
This is the integral solved by the BORE program;
Werner

I'm sitting down now going over the Bore program. Using Free42 on my iPhone it is instant. Haven't loaded up the DM42 yet but I am guessing it will be the same result.

I have to thank you again for this because You have pointed out a few commands I didn't know about and saved me some time. PGMINT being one. Found it in the manual on Pages 203-204. The 42s continues to impress but the knowledge on this forum is more impressive.

[Werner]

I just wanted to make it as simple as possible, because (as I said in my first post) you can use the interactive integrator as well, but since 'c' is used as upper limit, it is slightly less user-friendly.
Cheers, Werner

[DM48]

(05-28-2020 07:28 AM)Werner Wrote:  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

Cheers, Werner

I am going to ask some really basic and dumb questions about the above code.

Lines 00 - 08 I understand.

Line 09 This calls which program to integrate from within the program? Integrate "FX" which will be between lines 19 - 28?

Lines 10 - 12 This brings up the VARMENU that populates the softmenu with MVAR data and stops (pauses) the program. This allows me to enter the data (be interactive).

Lines 13 - 14 puts C in the X register and stores it to the ULIM (upper limit) and leaves C in the x register.

Line 15 tells what variable to integrate with respect to. X in this case.

How am I doing so far?

A few notes. a, b, and c are radii. On my 48G program, I would enter the diameters, 48, 58 and 24. The first thing the program does is half each one of them and stores them in a local variable It was because my integration was a quarter rotation.

I bring this up because of Lines 16 and 17.

Hear you multiply the ULIM or C by 4. Why? I feel this is related to my explanation above and the fact I am multiplying the final answer by 4.

I'll stop here.

[Werner]

All correct.
And I thought you needed the final answer *4.. delete the lines if you don't need them.

Cheers, Werner

[Albert Chan]

(05-28-2020 11:00 AM)DM48 Wrote:  It starts out as a triple integral. I’ll find my notes.

Please do.

I asked months ago about bead-volume, with square hole, but still don't know how to solve.
Your notes may help ...

(05-28-2020 12:14 PM)DM48 Wrote:  Can you get this to the final stage where the solver isn't needed?

This version required Hypergeometric2F1, see https://www.hpmuseum.org/forum/thread-14148.html

Mathematica 4.0:

In[1]:= {a,b,c} = {24, 29, 12} * 2;              (* all diameter *)
In[2]:= cv[d_,h_] := Pi (d/2)^2 h                 (* cyclinder volume *)

In[3]:= N[cv[c,b] - cv[c,a], 20]                     (* estimated bore-hole volume *)
Out[3]= 4523.8934211693022634

In[4]:= hv[d_,D_] := cv[d,D] Hypergeometric2F1[-1/2, 1/2, 2, (d/D)^2]

In[5]:= N[hv[c,b] - hv[c,a], 20]                    (* actual bore-hole volume *)
Out[5]= 4651.2464366171070976

Note: above assumed drilling through the pipe, penetrating both side.
For 1 sided bore-hole removed volume, cut above in half.

Update: found an equivalent hv[d,D], using Elliptic functions (#12)

hv2[d_,D_] := D/3 ((D^2+d^2) EllipticE[(d/D)^2] - (D^2-d^2) EllipticK[(d/D)^2])

[DM48]

(05-28-2020 04:06 PM)Werner Wrote:  All correct.
And I thought you needed the final answer *4.. delete the lines if you don't need them.

Cheers, Werner

Werner, I am just trying to understand the program. This is a lot different than RPL. I do need the *4, but wasn't sure why it was included there. I'm still not sure why it is but I am a little slow sometimes.

I will need to adjust the program for diameter entry (divide a, b and c by 2) and wanted to be sure I understood what the program was doing to that point. I am going through the next bit which is giving me trouble understanding.

[DM48]

(05-28-2020 04:40 PM)Albert Chan Wrote:  [quote='DM48' pid='132375' dateline='1590663642']
It starts out as a triple integral. I’ll find my notes.

Please do.

I asked months ago about bead-volume, with square hole, but still don't know how to solve.
Your notes may help ...

(05-28-2020 12:14 PM)DM48 Wrote:  Can you get this to the final stage where the solver isn't needed?

Ok. I'll look for them tonight. I have another problem that solves for the "bore hole" that has a taper to it.

[Werner]

The 4x was there as it is then executed only once, after the integral was calculated.
And there is a very easy way of adapting the program to diameters: just add the lines

Code:

2
/

right after LBL 01 and right before STO “ULIM”

Or, even simpler: just replace the 4 x by 2 /
(A volume with all parameters doubled is 8 times as large ;-)
Cheers, Werner

[DM48]

(05-28-2020 05:30 PM)Werner Wrote:  The 4x was there as it is then executed only once, after the integral was calculated.
And there is a very easy way of adapting the program to diameters: just add the lines

Code:

2
/

right after LBL 01 and right before STO “ULIM”

Or, even simpler: just replace the 4 x by 2 /
(A volume with all parameters doubled is 8 times as large ;-)
Cheers, Werner

I went the 2 / route you mentioned but now I am going to change that based on the halving four. Less lines of code. :-)

PGMINT readies the integration because it is in a running program and then INTEG "X" starts the integration.

Lines 19 - 34 builds the equation to be solved in the X register when INTEG "X" is called. Correct?

I am going all in on the 42s so I hope I won't always be asking obvious questions. RPL programming is linear so that is how I think.

[Albert Chan]

If we assume drilling though the pipe, penetrating both side, we do not need "2 /" correction.
In other words, I = bore-hole removed volume

XCas> [a,b,c] := [24,29,12] .* 2      // convert to diameter values
XCas> gaussquad(sqrt(c^2-x^2) * (sqrt(b^2-x^2) - sqrt(a^2-x^2)), x = 0 .. c) → 4651.24643662

// intersection of 2 cyclinders with radius 1 (diameter 2), volume = 16/3

XCas> b := c := 2
XCas> int(sqrt(c^2-x^2) * sqrt(b^2-x^2), x = 0 .. c) → 16/3

[Albert Chan]

Borrowing my AGM2 method of calculating ellipse perimeter, but factor out AGM.

Code:

local sqrt, pi = math.sqrt, math.pi

function agm(a, b)
    local t = 0.5
    local s = a*a
    while true do
        local k = b-a
        local s1 = s - t*k*k
        if s == s1 then return a+0.5*k, s end
        b = sqrt(a*b)
        a = a + 0.5*k
        t = t + t
        s = s1
    end
end

function perimeter(a, b)    -- ellipse perimeter
    a, b = agm(a+b, 2*sqrt(a*b))
    return pi*b/a           -- b/a = "diameter"
end

We can get *both* K(m), E(m) from 1 call of agm(), thus bore-hole removed volume.
see equation 12: https://mathworld.wolfram.com/SteinmetzSolid.html

Code:

function KE(m)    -- return K(m), E(m), m < 1
    m = sqrt(1-m*m)
    local K, E = agm(1+m, 2*sqrt(m))
    K = pi/K
    return K, K*E/4
end

function hv(d,D)  -- bore-hole removed volume, thru solid pipe
    local K, E = KE(d/D)
    return D/3 * (d*d*(K+E) - D*D*(K-E))
end

lua> a, b, c = 24*2, 29*2, 12*2    -- all diameters
lua> hv(c,b) - hv(c,a)                   -- assumed penetrating pipe both side
4651.24643661709

[John Keith]

(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 the reference, it is very informative. I would highly recommend users of later HP calculators, from the HP-28 through the Prime, read this as well. It explains the operation of these commands and PROOT (polynomial roots) as well as the underlying algorithms.

[Werner]

Well, I had to try that out, Albert!
It's more than 10x faster than the integral solution: (which is noticeable on a DM42):

Code:

00 { 207-Byte Prgm }
01▸LBL "BORE"
02 MVAR "A"
03 MVAR "B"
04 MVAR "C"
05▸LBL 10
06 VARMENU "BORE"
07 RTN
08 RCL "C"
09 RCL "A"
10 XEQ "HV"
11 LSTO "t"
12 RCL "C"
13 RCL "B"
14 XEQ "HV"
15 RCL- "t"
16 GTO 10
17▸LBL "AGM2"
18 LSTO "B"
19 R↓
20 LSTO "A"
21 X↑2
22 LSTO "S"
23 0.5
24 LSTO "T"
25▸LBL 02
26 RCL "A"
27 RCL× "B"
28 SQRT
29 X<> "B"
30 RCL- "A"
31 RCL ST X
32 2
33 ÷
34 STO+ "A"
35 RCL "S"
36 STO ST Y
37 RCL "T"
38 STO+ "T"
39 R↑
40 X↑2
41 ×
42 -
43 STO "S"
44 X≠Y?
45 GTO 02
46 RCL "A"
47 X<>Y
48 RTN
49▸LBL "KE"
50 X↑2
51 1
52 -
53 +/-
54 SQRT
55 1
56 X<>Y
57 STO+ ST Y
58 SQRT
59 STO+ ST X
60 XEQ "AGM2"
61 X<>Y
62 PI
63 X<>Y
64 ÷
65 STO× ST Y
66 X<>Y
67 4
68 ÷
69 RTN
70▸LBL "HV"
71 LSTO "D"
72 X<>Y
73 LSTO "d"
74 X<>Y
75 ÷
76 XEQ "KE"
77 RCL- ST Y
78 X<>Y
79 LASTX
80 +
81 RCL "d"
82 X↑2
83 ×
84 X<>Y
85 RCL "D"
86 X↑2
87 ×
88 +
89 RCL× "D"
90 3
91 ÷
92 END

Cheers, Werner

[DM48]

Albert, thank you for posting that update. That is over my head but it allowed Werner to create a noticeably faster implementation of Bore on my DM42 and for that I am grateful to the both of you.

I have the "coupon" that was "tapped" from ductile iron pipe on my desk. I have always assumed, when flattened, it would be an ellipse. Is that correct?

*Update 1* I ran Werner's original Bore program on an authentic HP42s. It took approximately 30 seconds to complete. Will program in this new program to compare.

*Update 2* Reading closer through the program I found LSTO command. Guessing this is a local variable command not found on the HP 42s?

[Albert Chan]

(06-01-2020 01:02 PM)Werner Wrote:  Well, I had to try that out, Albert!
It's more than 10x faster than the integral solution: (which is noticeable on a DM42):

Very nice !

Suggestion:
Your AGM2 code had AGM(a,b) in Y register, instead of expected X
The fix actually made the code shorter, by deleting Line 47, Line 61

To make the name more suggestive, rename "AGM2" as simply "AGM"
With these changes, on Free42, we calculate Pi via AGM, like this.

0.5 SQRT 1 XEQ "AGM" [X^2] [X<>Y] / 2 *
→ 3.141592653589793238462643383279504

Another example, ellipse_perimeter(50, 10):

50 Enter 10 +
50 Enter 10 * SQRT 2 *
XEQ "AGM" /                             ; "diameter" of ellipse = y/x
PI *
→ 210.100445397

Or,

50 Enter 10 XEQ "AGM"
[X<>Y] 10 [X↑2] + [X<>Y] /     ; "diameter" of ellipse = (y+b²)/x
PI *
→ 210.100445397

For the same reason, I would rename "KE" as "EK", since E is in X register

Again, for ellipse_perimeter(50, 10):

Perimeter = 4 a E(e), where a = 50, e = √(1-0.2^2) = √0.96

0.96 SQRT
XEQ "EK"                     ; E(e)
200 *                          ; 4 a
→ 210.100445397

What do you think ?

[Werner]

Albert: well it was just my first try, and it worked straight away - this doesn't happen often ;-) so I was eager to post it.
And I called it AGM2 because I already have a stack-only version of AGM, that returns only the agm.

DM48: yes, LSTO is one of the Free42 extensions, Local variables that are automatically cleared at RTN. Just replace it by STO for the 42S.
(ah that will not work: I see I have "A" and "B" twice there. So, change the A,B,C in BORE to lowercase). Along with Albert's suggestions and one more small change, here's the latest version:

Code:

00 { 204-Byte Prgm }
01▸LBL "BORE"
02 MVAR "a"
03 MVAR "b"
04 MVAR "c"
05▸LBL 10
06 VARMENU "BORE"
07 RTN
08 RCL "c"
09 RCL "a"
10 XEQ "HV"
11 LSTO "t"
12 RCL "c"
13 RCL "b"
14 XEQ "HV"
15 RCL- "t"
16 GTO 10
17▸LBL "AGM2"
18 LSTO "B"
19 R↓
20 LSTO "A"
21 X↑2
22 LSTO "S"
23 0.5
24 LSTO "T"
25▸LBL 02
26 RCL "A"
27 RCL× "B"
28 SQRT
29 X<> "B"
30 RCL- "A"
31 RCL ST X
32 2
33 ÷
34 STO+ "A"
35 RCL "S"
36 RCL "T"
37 STO+ "T"
38 R↑
39 X↑2
40 ×
41 -
42 ENTER
43 X<> "S"
44 X≠Y?
45 GTO 02
46 RCL "A"
47 RTN
48▸LBL "EK"
49 X↑2
50 1
51 -
52 +/-
53 SQRT
54 1
55 X<>Y
56 STO+ ST Y
57 SQRT
58 STO+ ST X
59 XEQ "AGM2"
60 PI
61 X<>Y
62 ÷
63 STO× ST Y
64 X<>Y
65 4
66 ÷
67 RTN
68▸LBL "HV"
69 LSTO "D"
70 X<>Y
71 LSTO "d"
72 X<>Y
73 ÷
74 XEQ "EK"
75 RCL- ST Y
76 X<>Y
77 LASTX
78 +
79 RCL "d"
80 X↑2
81 ×
82 X<>Y
83 RCL "D"
84 X↑2
85 ×
86 +
87 RCL× "D"
88 3
89 ÷
90 END

Cheers, Werner

[ijabbott]

(06-01-2020 04:46 PM)DM48 Wrote:  Albert, thank you for posting that update. That is over my head but it allowed Werner to create a noticeably faster implementation of Bore on my DM42 and for that I am grateful to the both of you.

I have the "coupon" that was "tapped" from ductile iron pipe on my desk. I have always assumed, when flattened, it would be an ellipse. Is that correct?

I don't think so. If you define variables \(R_p\) for the pipe radius and \(R_h\) for the hole radius, the parametric equations for the "coupon" would be:

\[\begin{cases}
X_c = R_p \arctan\left(\frac{R_h \cos(t)}{\sqrt{{R_p}^2 - (R_h \cos(t))^2}}\right) \\
Y_c = R_h \sin(t)
\end{cases}
t \in [0, 2\pi)\]

EDIT: Corrected it, I think.

But the parametric equations for an ellipse with the same axes would be:

\[\begin{cases}
X_e = R_p \arcsin\left(\frac{R_h}{R_p}\right) \cos(t) \\
Y_e = R_h \sin(t)
\end{cases}
t \in [0, 2\pi)\]

[Werner]

Continuous improvement.. when I rewrite agm as

function agm2(a, b)
local t = -1
local s1 = a*a
do
local s = s1
local k = (b-a)/2
t = t + t
b = sqrt(a*b)
a = a + k
s1 = s + t*k*k
until s == s1
return a, s
end

and I switch the order of the arguments in hv, I get

Code:

00 { 192-Byte Prgm }
01▸LBL "BORE"
02 MVAR "a"
03 MVAR "b"
04 MVAR "c"
05▸LBL 10
06 VARMENU "BORE"
07 RTN
08 RCL "a"
09 XEQ 03
10 LSTO "t"
11 RCL "b"
12 XEQ 03
13 RCL- "t"
14 GTO 10
15▸LBL 03
16 RCL "c"
17▸LBL "HV"
18 LSTO "d"
19 X<>Y
20 LSTO "D"
21 ÷
22 XEQ "EK"
23 RCL- ST Y
24 X<>Y
25 LASTX
26 +
27 RCL "d"
28 X↑2
29 ×
30 X<>Y
31 RCL "D"
32 X↑2
33 ×
34 +
35 RCL× "D"
36 3
37 ÷
38 RTN
39▸LBL "EK"
40 X↑2
41 1
42 -
43 +/-
44 SQRT
45 1
46 X<>Y
47 STO+ ST Y
48 SQRT
49 STO+ ST X
50 XEQ "AGM2"
51 PI
52 X<>Y
53 ÷
54 STO× ST Y
55 X<>Y
56 4
57 ÷
58 RTN
59▸LBL "AGM2"
60 LSTO "B"
61 R↓
62 LSTO "A"
63 X↑2
64 LSTO "S"
65 SIGN
66 +/-
67 LSTO "T"
68▸LBL 02
69 RCL "A"
70 RCL× "B"
71 SQRT
72 X<> "B"
73 RCL- "A"
74 2
75 STO× "T"
76 ÷
77 STO+ "A"
78 X↑2
79 RCL× "T"
80 RCL+ "S"
81 ENTER
82 X<> "S"
83 X≠Y?
84 GTO 02
85 RCL "A"
86 END

Of course, if you don't need the alpa labels AGM2, EK and HV, you can replace them with numeric ones ( leave out LBL "HV" entirely), saving another 22 bytes, I think.

Cheers, Werner