Programming Challenge: a classic trigonometry problem

[Albert Chan]

(03-08-2023 04:09 PM)Albert Chan Wrote:  (1) and (2): (1+1/p)² * (1-p²) = (a²-b²)/c² = k      .... (3)

Equation (3), p and c have positive relationship.

(c → min(c) = 0) ⇒ (p → 0)
(c → max(c) = 1/(1/a+1/b) = (a*b)/(a+b)) ⇒ (p → b/a)      .... (4)

We can estimate max(p), assuming (3) term (1-p²) = 1

From (3) and (4): p ≤ min(1/(√k - 1), b/a)

With maximum guess, it is possible secant line over-shoot passed y-axis.
It may blow up, or converge to the "wrong" root.

It would be better if we have a minimum guess instead.
This way, iterations converge away from y-axis, making it more robust.

Let q = 1/p = u/v      // p max guess = q min guess

Q = (1+q)^2*(1-q^-2) - k

Q = 0       ⇒ |q| and k have positive relationship.

k = 16 ε³  ⇒ revert taylor series, q ≈ -horner([2,2,1],ε), +horner([2,2,1],ε³)
k → ∞      ⇒ Q ≈ q^2 - k = 0,     q ≈ ±√k

I have modified secant verbose option to display iterations, 1 per line.

I also added option for negative eps, to swap (x1,x2) if |f(x1)| < |f(x2)|
"Worse" point is moved as the first, thrown-away ASAP.

lua> S = require'solver'
lua> a, b, c = 40, 30, 15
lua> k = (a*a-b*b)/(c*c)
lua> Q = fn'q: (1+q)^2*(1-q^-2) - k'
lua> eps = -1e-9 -- swap worse point as first

lua> lo = max(sqrt(k)-1, a/b) -- ladder crossed solution
lua> S.secant(Q, lo, 1.1*lo, eps, true)
1.3333333333333333 -- = lo
1.4666666666666668
1.4445731516722746
1.444399371124978
1.444399566520114
1.4443995665182763

lua> q = _
lua> u = c*(1+q)
lua> x = sqrt((a-u)*(a+u))
lua> x, u, a -- right triangle
15.987649008568585      36.665993497774146      40

Extended-line crossed solution (negative q), we may just flip the sign of lo, and solve.
It take more iterations than with good guess, but it can be done. Over-shoot is harmless.

lua> lo = -lo -- bad guess, but OK
lua> S.secant(Q, lo, 1.1*lo, eps, true)
-1.3333333333333333 -- = lo
-1.4666666666666668
-7.344688904660508
-1.9135094120851166
-2.263168070942035
-3.2061293483685995
-2.8172286303461243
-2.875248024422489
-2.8808890380615035
-2.8807851830598548
-2.8807853490119513
-2.880785349016917
-2.880785349016917

lua> lo = -max(sqrt(k)+1, a/b) -- good guess
lua> S.secant(Q, lo, 1.1*lo, eps, true)
-3.040217628114033
-2.7638342073763935 -- = lo, swapped here
-2.875471300898932
-2.880971458358896
-2.880785063180847
-2.8807853490015667
-2.880785349016917
-2.880785349016917

lua> q = _
lua> u = c*(1+q)
lua> x = sqrt((a-u)*(a+u))
lua> x, u, a -- right triangle
28.3565769435901        -28.211780235253755      40

solve (1+q)^2 * (1-q^-2) * UnitStep[q^2-1] = 28/9, for q

Comment: since k ≥ 0, plot only show |q| = |u/v| ≥ a/b ≥ 1

---

From my next post, I suggested a correction, inside sqrt: (v/c)² = (1+1/q)²

Correction = 1 if q=±∞, but k ≈ q², even bigger. Correction doesn't matter.
Correction = 4 if q=1, it is worthwhile to use. Bonus: now hi ≥ 1

lua> hi = -1 + sqrt(k+4)
lua> S.secant(Q, hi, hi*0.9, eps, true)
1.6666666666666665 -- = hi
1.5
1.4442155853072434
1.4444000308015605
1.4443995665229092
1.4443995665182763
1.4443995665182763

Correction = 0 if q=-1., lo = -1 - sqrt(k) seems OK.
But when k is tiny, lo guess might be "too good".
Q'(q=-1) = 0. Iterations may bounce around wildly.
Last 2 points are extrapolated again for secant root, which may be garbage.
Going from hi side, steadily reach the true root is safer. (Q is concave-up)

lua> hi = -1 - sqrt(k+1)
lua> S.secant(Q, hi, hi*0.9, eps, true)
-3.0275875100994063 -- = hi
-2.7248287590894655
-2.874159866479985
-2.881098675537131
-2.880784748823651
-2.8807853489626534
-2.880785349016917
-2.880785349016917