[VA] SRC #012d - Then and Now: Area

[Valentin Albillo]

      
Hi, all,

Well, almost two weeks have elapsed since I posted my OP, which has already passed the ~1,400 views mark (alas, once again when the difficulty increases the number of posts and/or views decreases,) and we've got a number of solutions and comments (though no RPL ones, I wonder why ...), namely by PeterP, Martin Hepperle, Werner, J-F Garnier, Fernando del Rey, and Albert Chan. Thank you very much to all of you for your interest and valuable contributions.

Now, this is my detailed sleuthing process and resulting original solution, plus additional comments:

Note: though a number of graphics are featured in this solution, as suggested in the problem's statement, rest assured that none of them are needed, all intermediate and final results can be automatically computed without relying on graphics at all. The ultimate reason behind using them is that I feel they greatly help interested readers to better understand the steps to the solution, and they make for a more enjoyable presentation as compared to a wall of dry, terse math expressions alone.

My sleuthing process

First of all, we rearrange the equation into this equivalent, simpler form:
      
Now, as suggested in my OP, we use a simple function plotter to get a visual idea of the region's shape and determine a rough (d.dd) approximate value for its ymin and ymax limits:

Note: There's no need to get now the xmin and xmax limits because we'll compute them as a function of the y values, as will be seen below.

and at first sight it seems that the region is a single island between ymin = 0 and ymax ~ 2.83, but there's the nagging question: How can we be sure that there's only just one island ? Perhaps it might be the case that the region consists of several disconnected subregions. How can we know for sure ?

One obvious way would be to obtain surefire limits for the valid ranges for x and y and carefully zoom and pan over the corresponding enclosing rectangle or else write scanning code to try and detect additional subregions, but this would be extremely time-consuming if the hypothetical subregions ("islands") happen to be very small (or worse, nonexistent,) which might well be the case as we can't see them outright.

At any rate, even if our zooming and/or scanning did detect some small islands, we still wouldn't be sure that other even smaller ones didn't exist, so a more analytical approach is required and as the variable x appears just once, we isolate it, obtaining:
      
and it's clear that x will be real only if the expression under the square root is ≥ 0. Two cases:

● If > 0, we'll have two distinct values, x₁ and x₂, and their difference will be the width of the island for a given value of y.

● If = 0, then x₁ = x₂ and the width of the island will be 0 (i.e. the extreme top and bottom of the shape), so the corresponding values of y are the limits ymin and ymax for each existing island and we'll compute them by simply finding all the roots of the equation:
      

Graphing the expression will give us the number of islands and a little zooming will provide initial approximations (d.dd) to the roots, i.e. the ymin, ymax limits for each island:


so we see that there's definitely one big island with ymin = 0 and ymax ~ 2.83, and possibly a much smaller island near y ~ -4. Zooming a little more we have:



and indeed there's a very small island with ymin ~ -4.08 and ymax ~ -4.05 and no more islands are possible as the expression goes increasingly negative for y > ~ 2.83 and y < ~ -4.08. Now that we know the small island's location (the x can be obtained from the y using the formula above), we can visualize it:



It's worth noting that the formulas above for x₁ and x₂, if plotted independently, allows us to see how the perimeter of each island is composed of two separate "branches" which are joined seamlessly at their extremes:



Now, knowing the y-range [ymin, ymax] for an island, its area is given by the expression:
      
and the total area is the sum of the areas of both islands. The small island's integral is quite amenable to numerical integration save for the fact that it might lose 2-3 digits to cancellation, but as the resulting area is then added to the much larger one corresponding to the big island, those digits don't affect the sum to 12-digit accuracy.

The integral for the big island, however, has two singularities, one at y = 0 (which doesn't cause any problems as this is the ymin extreme and the Width function being integrated is never evaluated there by INTEGRAL,) and another inside the range [ymin, ymax], which causes INTEGRAL to be at least 400% slower and even so it can only produce a value accurate to at most 7-8 digits, not 12.

The singularities can be recognized as the y values where the graph's tangent is horizontal, i.e., the graph goes parallel to the X axis. The y value of the inside singularity can be obtained by solving the equation SQR(-LN(Y*Y/30.07+EXP(-SIN(Y)))) = Y (which gives ysng ~ 0.83, as can be seen in the zoomed graphic below:



and we can then split the single problematic integral into two parts, the singularities being now at the extremes of the ranges of integration where they do no harm, like this:
      

My original solution

My original solution is this little 5-line, 12-statement, 244-byte HP-71B program (it can be made faster by using an additional line of code, but I like it this way):   

1  DESTROY ALL @ DIM U @ H=1/3 @ D=FNROOT(1,1,SQR(FNG(FVAR))-FVAR)
2  DISP "Area:";FNI(0,D)+FNI(D,FNY(2.83))+FNI(FNY(-4.08),FNY(-4.05))

3  DEF FNG(Y)=-LN(Y*Y/30.07+EXP(-SIN(Y))) @ DEF FNY(Y)=FNROOT(Y,Y,FNG(FVAR))
4  DEF FNI(A,B)=INTEGRAL(A,B,1/10^10,FNW(IVAR)) @ DEF FNR(X)=SGN(X)*ABS(X)^H
5  DEF FNW(Y) @ U=SQR(FNG(Y)) @ FNW=FNR(Y+U)-FNR(Y-U)

Note: Notice that the parameter d = 1.598 does not appear in the solution's code, as the way it's featured in the inequality means that varying it only translates the whole region R along the X axis, which obviously doesn't alter R's area at all. Also, FNROOT/FVAR and INTEGRAL/IVAR are keywords from the Math ROM).

Line 1 performs some initialization and also computes ysng (DIM U is necessary to create variable U here and not inside DEF FNW, where we would fall prey to a well-known system bug).

Line 2 computes and displays the total area by evaluating and adding together the two integrals for the big island and the integral for the small island, while refining to full accuracy on the fly the d.dd approximations we got for the various ymin, ymax.

Line 3 defines FNG for the expression under the square root as a function of y, as well as FNY, which refines to full accuracy the passed d.dd y approximation.

Line 4 defines FNI, which returns the integral of the Width function between specified limits to full accuracy, as well as FNR, which returns the real cubic root (with the correct sign) of the real argument passed to it.

Line 5 defines FNW, which returns the Width of an island at a given y value.

Let's run it:

>STD @ RADIANS @ RUN

            Area: 2.07669834394

The 20-digit value is 2.07669834394 76059651..., so we've got a full 12 correct digits.


Note:  all 20-digit results given here were obtained using my original solution converted to RPN and run on Free42 Decimal. See Additional comments below.

We can obtain additional results if immediately after running the unmodified program above we execute these statements directly from the command prompt:

>RES-IVALUE ►  2.07662636775     (big island's area   :   2.07662636775 45636287...)
>IVALUE     ►  7.19761930517E-5  (little island's area:   7.19761930423 36384588...E-5)
>D          ►   .831971149978    (big island's ysng   :   0.831971149979 07679799...)
>FNY(2.83)  ►  2.82740261413     (big island's ymax   :   2.82740261412 95600918...)
>FNY(-4.08) ► -4.08514674778     (little island's ymin:  -4.08514674763 83097474...)
>FNY(-4.05) ► -4.04921226429     (little island's ymax:  -4.04921226439 69504765...)

For timing, execute instead:

>SETTIME 0 @ CALL @ TIME

Additional comments

● It's quite trivial to convert my original BASIC solution for the HP-71B to HP-42S/Free42/DM42's RPN code, matter of fact it can be done on the fly while typing the RPN instructions on the calculator/simulator, without ever needing to previously write down the RPN code on paper or a notepad file. I did it precisely that way, thanks to the utmost simplicity of classic 4-level RPN code.

This allows for greatly enhanced accuracy if using Free42 Decimal/DM42, as seen in this image, where the total Area has been computed and displayed to 25 correct digits. A bit of the resulting RPN code is shown as well but I'm not including the complete conversion to RPN here, that's left as an easy exercise for the reader.



● As some of you may have guessed (some did not), the problem's statement (gas turbine's efficiency, ethane fuel, etc.) makes absolutely no sense from a "dimensional point of view" and it's but a made up story I concocted in the spirit of those the great late Martin Gardner used to create as a pleasant wrapper around many of his puzzles and mathematical recreations, to make them all the more interesting (e.g., see "The Magic Numbers of Dr [I.J.] Matrix").

For this particular problem, I first carefully selected the two parameters (d = 1.598 and M = 30.07) to achieve the desired effect, and only afterwards did I look for real-life subject matters which would include those values. I quickly found that some ethane properties actually did, and the problem was thus stated in terms of ethane, further deciding on its use as fuel for some hypothetical gas turbine. Should I have found first that certain properties of herrings featured those numbers, the problem's description would've been built around herrings instead, possibly red ones.

And what's the "desired effect" ? Well, to make the little island as inconspicuous as possible without overdoing it. To that effect, each parameter controls a particular aspect of the inconspicuousness, namely:

○  The parameter M controls the size of the region R, including the area of the main island and the number and areas of any existing additional islands, if any. Small variations of M result in a very noticeable effect, see for instance the very large increase in the size of the small island when changing from the problem's M = 30.07 to M = 30.5:



You may want to verify that by decreasing M the area of the small island tends to 0, and for M < ~30.06480779965 it disappears altogether so there's only the main island left. On the other hand, if we keep increasing M the small "microscopic" island grows to "macroscopic" size, becoming quite conspicuous, and as M tends to infinity further islands of various sizes do appear.

○  The parameter d controls the position of the region R along the X axis, and the problem's value d = 1.598 places the small island just under the Y axis, which will tend to greatly obscure and hide it if the user plots the Y axis as usual, even if it would be otherwise quite visible, as seen in this image where d = 1.818 clearly shows the small island while the problem's d = 1.598 almost completely hides it under the Y axis, though its size is the same.

Well, that's all for now, I hope you enjoyed it. Awesome-but-slightly-more-difficult Problem 5 will be posted next February.

P.S.: If you want to discuss other approaches, show any mathematical analysis (math "lectures"), use any tools (Mathematica, etc.), languages (Python, etc.), calculators (HP Prime, etc.) or computers (SHARP Pocket computers, laptops, etc.) please post your comments and/or code to this parallel thread kindly created by EdS2.

V.