[VA] SRC #012f - Then and Now: Angle

[Valentin Albillo]

      
Hi, all,

Almost two weeks have elapsed since I posted Problem 6 and, as expected, its perceived difficulty has resulted in very few posts, just a single solution (other than my original one, of course) by Fernando del Rey and worthwhile refinements to it by the indefatigable J-F Garnier. Thank you very much to both of you for your interest and overall outstanding contributions.

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


My sleuthing process

First of all, the relevant formula is Newton's equation for universal gravitation, further simplified as all seven masses and the gravitational constant G are to be taken as 1. Also, as force, acceleration, velocity and position all are two-component values, representing them using complex numbers is ideal and greatly simplifies the coding while also reducing the running time.

Using the formula we obtain the force exerted by the fixed stars upon the moving planet and thus the acceleration, velocity and position at any given time. I fully knew that those values can be obtained by numerically solving the relevant 2nd-order ordinary differential equation by using for instance a favorite of mine, the 4th-order Runge Kutta method (aka RK4).

However, though this would surely be the choice method for best efficiency, just for fun I decided to use a much simpler, more intuitive procedure, which consist of first computing the force (and thus the acceleration) by adding up the contributions from each star, then updating velocity and position by assuming that the acceleration remains constant within the (tiny) time increment.

This procedure is extremely simple to implement and requires a single call per step to the function which computes the force/acceleration while RK4 would require four calls per step, which for the current complicated function is quite time-consuming. On the other hand, RK4 has a local truncation error (LTE) of O(h⁵) while the simpler procedure's LTE is O(h²) so we'd be trading accuracy/runtime for algorithmic simplicity; however, for the modest accuracy goal it will likely be adequate enough.

After selecting the algorithm to use, I tried out assorted angles in order to decide on a range where the optimum angle would surely lie, as angles outside it would obviously result in wasteful, unnecessarily longer or even "retrograde" trajectories. Also, though I did notice, I didn't take advantage of the initial position's symmetry. See Notes below.

Now, I'd need to scan the range using a coarse increment (for speed; see Notes below for details on the various heuristics applied) while keeping track of the minimum distance to the destination and, if within a certain tolerance, a refining procedure would then be applied to improve the accuracy, followed by linearly interpolating both angle and time so as to exactly hit the destination (within a small tolerance).

Finally, once the scan is completed, the optimum angle and resulting minimum time among all suitable trajectories would be subjected to a finer refinement to further improve the accuracy. As seen below, this strategy ultimately results in 4-5 correct digits.


My original solution

My original solution is this 12-line, 703-byte HP-71B program: 

  1  DESTROY ALL @ OPTION BASE 0 @ DIM B,C,D,J,T,U @ COMPLEX X(6),A,P,Q,S,V,Y,Z @ SETTIME 0
  2  DATA .005,.01,.1,.01,1,9,(0,0),(2,-1),(2,0),(2,1),(3,-1),(3,0),(3,1),4,1

  3  FIX 2 @ READ H,F,M,L,V0,G,X,B,C @ Y=(B,C) @ FOR R=-.5 TO .75 STEP F @ IF FND(R)>M THEN 5
  4  W=FNROOT(R,R,FNY) @ DISP W;T;TIME$ @ IF T<G THEN N=W @ G=T
  5  NEXT R @ H=H/10 @ W=FNROOT(N,N,FNY) @ FIX 4 @ DISP "Min T:";T;"for A:";W;"(";TIME$;")" @ END

  6  Q=P @ A=0 @ FOR J=1 TO 6 @ Z=X(J)-P @ A=A+Z/ABS(Z)^3 @ NEXT J @ P=P+H*V+K*A @ V=V+H*A @ RETURN

  7  DEF FND(R) @ P=X(0) @ V=RECT((V0,R)) @ K=H*H/2 @ U=9 @ FOR T=H TO 3 STEP H @ R=REPT(P)
  8  D=U @ GOSUB 6 @ U=ABS(P-Y) @ IF U>D THEN 'N' ELSE IF REPT(P)<R THEN FND=9 @ END
  9  NEXT T @ 'N': T=T-H @ FND=D @ END DEF

 10  DEF FNY @ H=H/2 @ P=X(0) @ V=RECT((V0,FVAR)) @ K=H*H/2 @ T=0
 11  T=T+H @ GOSUB 6 @ IF REPT(P)<B THEN 11 ELSE S=P-Q @ D=(B-REPT(Q))/REPT(S)
 12  T=D*H+T-H @ H=H*2 @ D=D*IMPT(S)+IMPT(Q)-C @ FNY=D*(ABS(D)>L)

Notes:
      
• Lines 1-2 perform the initialization and define the necessary data. There's no need for a RADIANS declaration because no trigonometric functions are used but RECT, which with a complex argument always uses radians regardless of the angular mode.
    
  Also, I noticed the symmetry of the starting position (as J-F Garnier did as well) but despite the obvious speed up I decided not to depend on it so that my program would be able to solve generalized cases (including non-symmetric starting positions,) by simply defining all relevant parameters (initial coordinates, velocity, increments, limits, tolerances, etc.) in a single DATA statement.
    
• Lines 3-5 are the main program, which coarsely scans a plausible angle range and once a suitable one is found, it uses FNROOT to try and find more precisely the angle which makes the Y coordinate of the final position exactly match the Y of the destination (within a small tolerance).
    
  Once the scan is over, both a candidate optimum angle and a minimum time have been selected among the tentative trajectories, and then a last, finer refinement is performed to further increase the final solution's accuracy (both angle and minimum time).
    
• Line 6 implements a subroutine which is called from both FND and FNY (defined below) and essentially computes the force (and thus the acceleration) upon the moving planet due to the six fixed stars, plus it also updates its current speed and position. It was initially implemented as another UDF but converting it to a subroutine speeds up the process.
    
• Lines 7-9 define FND, which computes the resulting trajectory for a given angle and returns an approximation to the minimum distance to the destination, as well as the nearest position's coordinates and the time to reach it. Besides there being a limit on this time, FND further includes heuristics to abort the computation if (a) the current position's X coordinate ever goes backwards or if (b) the distance to the destination eventually increases, to avoid computing a trajectory that can't result in the sought-for minimum time.
    
• Lines 10-12 define FNY, which is used while refining the minimum distance for a given angle with greater accuracy than when using FND, and further it interpolates the final coordinates and time so that the X coordinate exactly matches the X of the destination.
    
• The MATH ROM is used for all complex-number variable definitions and operations, plus FNROOT is used for the angle refinements.
  

Let's run it:

>RUN
Angle   Time   Runtime
-----------------------
-0.46   2.15   00:01:26
-0.27   1.91   00:04:01   {minimum time}
-0.19   1.93   00:05:41
 0.33   1.97   00:09:28
 0.55   1.92   00:12:06   {second best}
 0.73   2.05   00:14:43

Min T: 1.9160 for A: -0.2745 (00:24:30) {go71b: 24'30", Emu71/Win: 3'13", HP-71B: ~ 52h}

The timing for a physical HP-71B may seem a bit steep but it's nothing that an AC adapter and a weekend can't deal with. Just start the program on Friday evening, forget about it and go enjoy the weekend, and by Monday 4 a.m. you'll have the solution displayed, waiting for you to wake up.

The correct 10-digit result is Angle = -0.2748812540, Time = 1.916096918, so we can gauge the accuracy obtained this way:

>FIX 9 @ T;W

      1.916021524  -0.274541755

and the errors are:

      (-0.274881254) - (-0.274541755) = 0.000339499     i.e., we got 4 digits (save 3 ulp)
      ( 1.916096918) - ( 1.916021524) = 0.000075394     i.e., we got 5 digits (save 0.75 ulp)

For the record, the correct answer rounded to 33-34 digits is:

Angle = -0.274881254026276715509629216558426
Time  =  1.916096918332004910307460694725436

Additional comments

The six angles obtained by my original solution result in the following trajectories, counterclockwise:

• As can be seen in the graphic above, the straightest, most direct paths are the ones for the angles A ~ -0.27 and A ~ 0.55 and matter of fact the latter is clearly the one having the shortest length but the former takes the shortest time, as asked, namely 1.9161 vs. 1.9176, probably due to the close encounter with the (2,0) star significantly increasing the planet's velocity.
  
  To better see it, an easy challenge extension would be to compute and display not only the arrival time but also the path length and final velocity at the destination for the various trajectories.
        
• I hereby "officially" declare J-F Garnier as the valedictorian of this SRC#12, as he posted solutions and/or significant contributions to all six problems, with Fernando del Rey a worthy salutatorian, as he participated in most problems and delivered an excellent solution for this last, allegedly most difficult one. Congratulations and thanks a lot to both of you.
        
• Last but not least, I (faintly) hoped that some users of vintage HP graphing calcs would produce RPL solutions, obtaining not just the correct numbers but also graphically depicting the resulting trajectories as shown above, but alas, no such luck. Perhaps it was too difficult, as compared to their usual dealings with this library, or that port, or creating yet another bunch of list utilities and such ...

Well, that's all, at long last SRC#12 has been completed 7 months after I posted Problem 1, and the point has been fully made.

I'm now taking a sabbatical from proposing further elaborate challenges. Thanks for your interest and best regards.

V.