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

[Fernando del Rey]

I’ve been working on this Problem 4 since Valentin posted this interesting challenge, but I have not yet reached a solution accurate to 10-12 correct digits as requested. Still, while not abiding to Valentín’s instruction to post only one message featuring the fully refined solution, I am going to share what I have done so far, as I have some doubts that I’ll be able to complete the work within the next few of days.

Problem 4 asks us to calculate the area in the (x,y) plane where the criteria (1) below is met:

1/e^(((x-d)^3-y)^2 ) > y^2/M+1/e^(sin⁡(y)) (1)

To simplify the equation a little bit (and also the code to follow), let’s define:

x’=x-d

Then we can represent equation (1) as:

1/e^((x'^3-y)^2 ) > y^2/M+1/e^(sin⁡(y))

Let’s then name:

1/e^((x'^3-y)^2 ) = f(x',y)

and:

y^2/M+1/e^(sin⁡(y)) = g(y)

We know that f(x’,y) will always be in the range from 0 to 1 as the denominator of the division is positive.

0<f(x',y)≤1

We also know that g(y) will be in the range from y2/M+e^-1 to ∞, therefore in the range from e^-1 to ∞.

e^(-1)<g(y)<∞

In the (f,g) plane, we can represent the area that meets the three criteria:

0<f≤1
g>e^(-1)
f>g

I wish I knew how to include an image I have drawn inside this post, but I don't know how to do it. The area in the (f,g) plane meeting the three criteria listed above is a triangle defined by the points A=(1,1), B=(e^-1,e^-1) and C=(1,e^-1).

From point A we can now derive:

1>g(y)

1>y^2/M+1/e

M(1-1/e)>y^2

|y|<√(M(1-1/e))=4.3598

-4.3598 < y < 4.3598

And from point B we can derive:

f(x’,y)>e-1

1/e^((x'^3-y)^2 ) >1/e → (x'^3-y)^2<1

|x'^3-y|<1

And taking the extremes of y calculated earlier, we get:

-1-4.3598<x'^3<1+4,3598

∛(-5.3598)<x'<∛5.3598

-1.75 < x’ < 1.75

To allow for some margin, we will use the following rectangle area where the original Problem 4 criteria is met:

-1.76 < x’ < 1.76
-4.36 < y < 4.36

Now, to obtain an approximate value of the area where the criteria defined by equation (1) is met, we can use the following program for the HP-71B:

10 X1=-1.76 @ X2=1.76 @ Y1=-4.36 @ Y2=4.36 @ INPUT "N? ";N
20 D1=(X2-X1)/N @ D2=(Y2-Y1)/N @ D=D1*D2 @ S=0 @ T0=TIME
30 FOR Y=Y1+D2/2 TO Y2-D2/2 STEP D2
40 G=Y*Y/30.07+EXP(-SIN(Y))
50 FOR X=X1+D1/2 TO X2-D1/2 STEP D1
60 IF EXP(-(X*X*X-Y)*(X*X*X-Y))>G THEN S=S+D
70 NEXT X
80 NEXT Y
90 DISP "Area=";S;TIME-T0

This program divides the rectangle bounded by the limits of x’ and y as defined above into NxN small rectangles and scans all small rectangles progressing row by row. If the criteria of equation (1) is met at the middle point of a given rectangle, its area is added to the summation S. If the criteria is not met, nothing is added to S.

At the end, the value S gives a rough approximation of the area where criteria (1) is met. Trying with different values of N, I get the followings times are for Emu71/Win running on my PC (which is about 970 times faster than the real 71B):

Grid size; Area; Time (s)
10x10; 2.455552; .02
20X20; 2.148608; .04
50x50; 2.11177472; .24
100x100; 2.09028864; .94
200x200; 2.07647616; 3.72
500x500; 2.0788703232; 23.61
1000x1000; 2.0765989376; 97.96
2000X2000; 2.07705168; 378.71

So, we seem to be converging to a value of around 2.077, but the convergence is too slow to make this brute-force method practical. We need to apply a smarter method, which could take advantage of the INTEGRAL function provided by the 71B Math ROM. But I need to work on it, and I’m sure others will post correct solutions before I have time to do so.

By adding DISP X;Y at the end of line 60 you can get a rough idea of the shape and position of the area meeting criteria (1), by looking at the pairs of values being displayed as the program runs.

60 IF EXP(-(X*X*X-Y)*(X*X*X-Y))>G THEN S=S+D @ DISP X;Y

If you make the grid accurate enough, say 500x500, you can see a small area at the upper left side of the rectangle where criteria (1) is met, and then a much larger area in the middle to lower half of the rectangle.

f you make the grid 2000x2000, the small area at the upper left corner can be seen with better definition.

With a grid of 100x100 you don’t get to see the small area in the upper left part, you will only see the large area starting towards the middle of the rectangle. It would be interesting to plot these areas, but I can't do it with a bare-bones 71B or with Emu71/Win.

It’s very possible that there be more than 2 different areas where the criteria (1) is met, but I have been able to “see” only those two so far.

My proposed way forward would be to express the equation bounding the area where criteria (1) is met:

1/e^((x'^3-y)^2 ) = y^2/M+1/e^(sin⁡(y))

Transforming it to isolate x’ as a function of y, and then using the INTEGRAL function of the Math ROM to do the calculation. But I sill need to figure out what the boundary values for the integration would be, and how to cater for the fact that there seems to be separate areas (at least two) where criteria (1) is met. So still some work to do!