Comments and discussion on Valentin's 4th "Then and Now" - Area

[EdS2]

(In these Comments and Discussion threads, of which this is the first, LENGTHY MATHS SESSIONS are welcome - provided of course the intent is to spread knowledge and have an enjoyable exchange of ideas. It's especially useful, I think, to adopt the convention of showing long decimal results with the correct and incorrect digits emphasised in some way: part of the interest with calculators, I think, is how much accuracy you really have, even if you always see lots of digits.)

[Albert Chan]

(01-22-2023 06:54 AM)C.Ret Wrote:  Simply add the counter variable C to the list of call indices: SUB SIMPSON(X1,X3,X5,F1,F3,F5,E,S,C) and increment C in Simpson's procedure each time it calls the function F.

Thanks! I was hoping for something similar to Python "global c", but this tip work!
Here is full listing, integrating main area.

10 X1=0 @ X3=2.82740261413 @ E=10^(-8)
20 T=TIME @ X2=(X3+X1)/2 @ CALL F(X1,F1) @ CALL F(X2,F2) @ CALL F(X3,F3)
30 S=(F1+4*F2+F3)*(X3-X1)/6 @ C=3 ! F calls counter
40 CALL SIMPSON(X1,X2,X3,F1,F2,F3,32*E,S,C) @ DISP S,C,TIME-T @ END

50 SUB F(X,Y) @ S=SQR(-LN(X*X/30.07+EXP(-SIN(X))))
60 Y=SGN(X+S)*ABS(X+S)^(1/3)-SGN(X-S)*ABS(X-S)^(1/3)

100 SUB SIMPSON(X1,X3,X5,F1,F3,F5,E,S,C)
110 H=(X5-X1)/4 @ X2=X1+H @ X4=X5-H
120 CALL F(X2,F2) @ CALL F(X4,F4) @ C=C+2
130 H=H/3 @ S1=(F1+4*F2+F3)*H @ S2=(F3+4*F4+F5)*H
140 S3=S1+S2 @ D=S3-S @ H=E/2
150 IF H<1.E-15 OR ABS(D)<H THEN S=S3+D/15 @ STOP
160 CALL SIMPSON(X1,X2,X3,F1,F2,F3,H,S1,C)
170 CALL SIMPSON(X3,X4,X5,F3,F4,F5,H,S2,C)
180 S=S1+S2

>run
 2.07662636771      2601      9.13 ! @200x, HP71B about 30 minutes

Overhead of adding counter = 9.13/8.96 - 1 ≈ 2%, not too bad
SIMPSON calls = ((F calls) - 3)/2 = 1299

Above accuracy is a fluke. Without extended precision, best we could hope for is about 10 digits.
E (max absolute error) should setup to match this. Smaller E may make the integral estimate worse.

[Albert Chan]

(01-21-2023 08:06 PM)Albert Chan Wrote:  I had also considered breaking up main area into more natural 2 pieces, and apply u-transformation.

\(\displaystyle \int_a^{a+b} = \int_0^a + \int_a^{a+b} \)

Code is more readable, because this does not involve "back folding".
But, I may need math to explain how u-substitution work, and why I need it here.
INTEGRAL built-in u-sub is not enough. (that's why JFG sum main area from 4 pieces)

To keep post short, I gave up on that idea.

This was my 2-splits main area solution, as opposed to 3-splits: \( \displaystyle
\int_0^{a+b} = \int_0^{a/2} - \int_a^{a/2} + \int_a^b\)

10 M=30.07 @ A=.831971149978 @ B=2.82740261413-A @ P=10^(-9)
20 T=1/3 @ DEF FND(Y,S)=SGN(Y+S)*ABS(Y+S)^T-SGN(Y-S)*ABS(Y-S)^T
30 DEF FNF(Y)=FND(Y,SQR(-LN(Y*Y/M+EXP(-SIN(Y)))))
40 DEF FNG(Z)=FNF(A*Z)*A+FNF(A+B*Z)*B
45 DEF FNU(U,U2)=FNG(U2*U*(4-3*U))*U2*(1-U)*12
50 SETTIME 0 @ DISP INTEGRAL(0,1,P,FNU(IVAR,IVAR*IVAR)),TIME

>run
 2.07662636771        .71      ! @200x, HP71B = 142 seconds
>p=1e-10 @ run 50
 2.07662636775      1.42      ! @200x, HP71B = 284 seconds

Both 3-splits and 2-splits are equally good.

3-splits use 3*127 = 381 FNF() calls to get 12 accurate digits.
2-splits use 2*127 = 254 calls to get 11+ digits, 2*255 = 510 to get full 12.

---

We wanted u-substitution to produce polynomial-like integrand, easy for INTEGRAL to handle.
Let x = U, where U is a function of u

\(\displaystyle \int _0^1 g(x)\;dx = \int _0^1 g(U)\,(U' \; du) \)

Based from plots, I had guessed g(x) both ends curve like cube root, left side more extreme.
To transform to polynomial like integrand, I cube left edge, and square the other.
We start from derivatives: (u^3)' = 3u^2, ((1-u)^2)' = -2(1-u)

\(→ U' = K\,u^2\,(1-u)\)

For u = 0 .. 1, U' has sign of K, a non-zero constant --> x and u are one-to-one.

\(→ U = \int U' du = K \left(\frac{u^3}{3} - \frac{u^4}{4} \right) + C \)

Matching integral limits, x=0 → u=0, x=1 → u=1, we have K = 12, C = 0

\(→ U = u^3\,(4-3u) \quad,\quad U' = 12u^2\,(1-u)\)

NOTE: U does not have to be a polynomial of u. see Kahan's Q(x) example

[Gjermund Skailand]

I found that double exponential integration is also well suited to evaluate the integrals.
With pre-calculated boundaries, and the DM42 running on batteries I can evaluate the integrals in about 3 seconds to 12 digits accuracy using only 169 function evaluations.
Reported relative accuracy of 6.8E-13 compared to actual relative accuracy of 3.7E-13 is also good.

However, when I increase the required accuracy to e.g 1E-20, then in 7.2 seconds, using only 403 function evaluations, I get a reported relative accuracy of 2.7E-26, but comparing with
2.07669834394 7601111175686517522841 (double exponential)
2.07669834394 76059651.. (Valentin Free42, standard integration?)

the result only agree up to 14-15 digits.

Pushing the limits to e.g 1E-30 does results in a reported accuracy of 4E-32, but the results does not really change.
2.07669834394 7601111175686517522 91 (16.7 seconds, 941 function evaluations)

So it is possible that the transformations cause a lot of cancellation but we got 32 digits to work with, and the requirement was 12 digits, which is quickly calculated in a few seconds without excessive evaluations.
Use is similar to using standard integrator, the DE-integration routine takes all inputs from standard 4-level stack
T: lower limit
Z: upper limit
Y: "function"
X: "X", variable to integrate
and returns
Y: error estimate
X: result
I have posted the DE program in .raw format with some description in the "General software library" https://www.hpmuseum.org/forum/thread-19547.html
Thanks to Didier for the zip tip.
best regards
Gjermund

[Albert Chan]

(02-12-2023 05:00 PM)Gjermund Skailand Wrote:  I found that double exponential integration is also well suited to evaluate the integrals.
...
Pushing the limits to e.g 1E-30 does results in a reported accuracy of 4E-32, but the results does not really change.
2.07669834394 7601111175686517522 91 (16.7 seconds, 941 function evaluations)

Thanks for DE integrate code!

A weakness with DE algorithm is we cannot have fuzzy integrand.
Above calculated area only had about half-precision accuracy.

Let's compare the same calculations, fuzzy vs "sharp" integrand:

>>> from mpmath import *
>>> mp.pretty = 1
>>> mp.dps = 34
>>> fnc = lambda x: sign(x)*abs(x)**(1/3)
>>> fnd = lambda y,s: fnc(y+s) - fnc(y-s)
>>> fne = lambda y: y*y*100/3007 + expm1(-sin(y))
>>> fnf = lambda y: fnd(y, sqrt(-log1p(fne(y))))
>>> main = [0, findroot(lambda y: fne(y)-expm1(-y*y),[0.5,1]), findroot(fne,[2.5,3])]
>>> speck = [findroot(fne, [-4.09,-4.08]), findroot(fne, [-4.05,-4.04])]

>>> main
[0, 0.8319711499790767979993060165734035, 2.82740261412956009186848044583741]
>>> speck
[-4.085146747638309747433733595149126, -4.049212264396950476554495084637107]

Note: mpmath quad use DE algorithm by default. Also, main had 2 parts, treated as 2 integrals.

>>> quad(fnf, main) + quad(fnf, speck)
2.076698343947601079521669308070185

The fuzziness comes from catastrophic cancellation of fnd
We redefine fnd, using identity: a^3 - b^3 = (a-b) * (a^2 + a*b + b^2)

>>> fnd = lambda y,s: 2*s/(fnc((y+s)**2) + fnc((y+s)*(y-s)) + fnc((y-s)**2))

>>> quad(fnf, main) + quad(fnf, speck)
2.076698343947601116938636757001192

[Gjermund Skailand]

Thanks Albert,

I was hoping for insight on the problem, and wishing for someone to calculate a result to high precision.
PS I did split the large area into two, so I actually did calculate 3 integrals to get the results.
br Gjermund

[Albert Chan]

Hi, Gjermund Skailand

My numbers were wrong!

I was using sign(x) * abs(x) ** (1/3) for cube root, but (1/3) evaluated as Python float.
I should have used mpf(1)/3 exponent, or simply use mpmath cbrt(x) = exp(log(x)/3)

For total area, fnd cancellation errors does not matter.

>>> fnc = lambda x: sign(x)*cbrt(abs(x))
>>> fnd = lambda y,s: fnc(y+s) - fnc(y-s)
>>> quad(fnf, main) + quad(fnf, speck)
2.076698343947601111175686517522914

Now, it match your numbers, almost exactly. Sorry for the noise.

[pier4r]

(01-20-2023 01:56 AM)Valentin Albillo Wrote:  which with the delay and the suggestion (not mandatory rule) that people would do best refining their initial solutions before posting them, to reduce the clutter in the thread, essentially prevented that learning process for him.

This being so, I'm not sure if the benefits outdo the disadvantages, so perhaps I will reinstate the 2-3 day delay for future problems, or perhaps I won't. What do you think ?

Since the forum is flexible and EdS2 (or anyone else) can create a sibling thread for "things that do not fit the rules", they could post the iterative refinement in such sibling topics, and only the final solution in the original one.
I agree that reading solutions that grow, like learning the historical process of some scientific/mathematical discoveries , may be easier to understand things.

Though one can have both, the process with the scaffolding, in one thread, and "final solution only, without the scaffolding" in the other (quoting "no self-respecting architect leaves the scaffolding in place after completing the building" of Gauss).

But that's my suggestion.

PS: I disagree with Gauss for the motives mentioned, from the scaffolding one can learn, because not everyone can make large comprehension steps at once, so seeing more steps could be beneficial.

[pier4r]

(01-20-2023 12:30 PM)Albert Chan Wrote:  I vote for no delay, and less rules, to encourage brainstorming among members.
Many times, OT ideas are more interesting than original puzzle!

The delay does not avoid spoiling the game ... only delay it.
Whoever post solution first, already spoil the game ... but someone has to do it!

But for that there are the sibling threads, it is perfectly fine (and I'd say even more ordered). Opening a new thread costs like, dunno, 5 minutes? After all the effort that Valentin (or whoever would do the same, but so far mostly Valentin does it) puts in the posts, I find opening a sibling thread very easy in comparison.

[EdS2]

It is true, and fortunate, that both Valentin and Albert are (by their responses) happy with the sibling threads idea.

It's also true that in the latest challenge offering, Valentin has chosen not to request a delay.

So, I think, and hope, there is no particular need to re-open the question. Perhaps I could respectfully suggest another thread, a new one, for discussing the merits and optimisation of how to corral or constrain the contributions to challenge threads. I'd say it won't be of interest to everyone who is interested in discussing the challenges themselves - it's a speciality subject, on the design of challenges.