Programming Challenge: a classic trigonometry problem

[Albert Chan]

(03-10-2023 03:53 PM)Albert Chan Wrote:  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.

With Secant's method, we want to avoid zero slope situation.
And, if q ≈ -1, u = c*(1+q) suffer subtraction cancellation.

For Q negative root, perhaps we should let r = -(1+q) ≥ 0

(1+q)^2 * (1 - q^-2) = k
r^2 - r^2/(1+r)^2 = k

Let s = r/(1+r) ≤ min(r, 1)

(r-s) * (r+s) = k
(r-s) = r - r/(1+r) = r*s = k/(r+s)

r*r ≥ r*s = k/(r+s) ≥ k/(r+r)
2 r³ ≥ k
r ≥ ³√(k/2)

From previous post, we have another: r ≤ √(k+1)      (*)

Solve for r is more robust (no zero slope issue), faster (less iterations), and better accuracy.
This setup is stable. We can even go for full convergence!

lua> S = require'solver'
lua> eps = 0 -- full convergence
lua> R = fn'r,s: s=r/(1+r); r*s - k/(r+s)'

lua> k = 1e9
lua> r = S.secant(R, cbrt(k/2), sqrt(k+1), -eps, true)
793.7005259840997
31622.776617495183
31622.776617495136
31622.776617494186
31622.776617494186

lua> k = 1e-9
lua> r = S.secant(R, cbrt(k/2), sqrt(k+1), -eps, true)
0.0007937005259840997
1.0000000005
-- points swapped
0.0007937020238029177
0.0007940155478145599
0.0007940155478967004

For k = 1e-9 --> R = 5 calls

Solve Q(k = 1e-9) = 0 for negative q, for comparison:
hi = -1-sqrt(k+1), 0.9*hi .. hi --> Q = 32 calls
lo = -1-sqrt(k),     lo .. 1.1*lo  --> Q = 795 calls!

(*) we also have r ≥ √k, so technically r ≥ max(√k, ³√(k/2))
But, I don't want too complicated limit, especially since calculating R is so cheap.
We are only interested in asymptote estimate when k in tiny, thus cube root only.

(01-07-2014 06:03 AM)Thomas Klemm Wrote:  Wikipedia provides yet another solution: Crossed ladders problem
\[ x^3(x-c)=1 \]

Another way is to solve above, for negative X root.
Let Y = -1/X > 0, we have:

Y^4 = C*Y + 1      → Y = ⁴√(C*Y+1)

If C is big, we can drop +1, and get Y ≈ ³√C
The formula is so stable, that same guess is good enough for small C
We can even use it as iteration formula, without using solver!

Again, we solve for r = -(1+q), indirectly from y, also full convergence.
To simplify, we assume c=1. In other words, all variables are dimensionless.

lua> Y = fn'y: y - (C*y+1)^0.25'

lua> k = 1e9
lua> C = 4/sqrt(k); lo = cbrt(C)
lua> y = S.secant(Y, lo, lo*1.01, -eps, true)
0.05019802884366822
0.0507000091321049
1.000031623628705
1.0000316222766017
lua> (y+1/y) * (2/C) -- = r
31622.776617494186

lua> k = 1e-9
lua> C = 4/sqrt(k); lo = cbrt(C)
lua> y = S.secant(Y, lo, lo*1.01, -eps, true)
50.19802884366822
50.7000091321049
-- points swapped
50.198031475628454
50.19803147889932
50.198031478899324
lua> (y+1/y) * (2/C) -- = r
0.0007940155478967004