Programming Challenge: a classic trigonometry problem

[Albert Chan]

(01-01-2014 08:44 AM)Thomas Klemm Wrote:  I used additional variables:
\[
u=\sqrt{a^2-x^2} \\
v=\sqrt{b^2-x^2}
\]...
\[
\frac{1}{u}+\frac{1}{v}=\frac{1}{c}
\]

Avoiding infinities, when u or v is zero, let \(f = c\,(u+v) - uv\)

\( \displaystyle \frac{d}{dx}\sqrt{k-x^2} = \frac{-x}{\sqrt{k-x^2}}\)

For both u and v, if we increase x, slope numerator goes bigger, but denominator goes smaller.
To straighten f curve, let \(y = x^2\)

\( \displaystyle \frac{d}{dy}\sqrt{k-y} = \frac{-1/2}{\sqrt{k-y}}\)

Assume a ≥ b > 0, and 1/c > 1/a + 1/b

\(f(y = 0) = c\,(a+b) - ab = abc \left(\frac{1}{a} + \frac{1}{b} - \frac{1}{c} \right) < 0\)

\(f(y = b^2) = c\,(u + 0) - u×0 > 0\)

We have a guaranteed real root between the y extremes.
Solved with Ridder's method (FYI, Free42/Plus42 SOLVE also use Ridder's method)

Note: each updated interval required cost of 2 function calls.
Note: last result is secant root of last interval, fall back to bisection, or f(x)=0

lua> S = require'solver' -- translated from Ridder solver Python code
lua> f = fn'y: local u,v = sqrt(a*a-y), sqrt(b*b-y); c*(u+v)-u*v'
lua> a, b, c = 40, 30, 15
lua> S.ridder(f, 0, b*b, 1e-9, true)
0                      900
0                      268.5639964554787
134.28199822773934     255.620329669562
194.95116394865067     255.6049251825406
225.27804456559562     255.60492082148892
255.6049208211841
lua> sqrt(_) -- = x
15.987649008568587

This work for full range c = 0 .. a*b/(a+b)

lua> c = 1
lua> S.ridder(f, 0, b*b, 1e-9, true)
0                      900
885.4877793349734      900
897.8097220135339      900
898.9165749684195      900
898.9165749684195      898.9214889155777
898.9199601065723      898.9214889155777
898.9199601065723      898.9199601069578
898.9199601069578
lua> sqrt(_) -- = x
29.981993931474232

lua> c = 17
lua> S.ridder(f, 0, b*b, 1e-9, true)
0                      900
0                      19.936187659463883
9.968093829731941      18.368152011527595
14.168122920629768     18.368139023904078
18.368139023879372
lua> sqrt(_) -- = x
4.285806694646805

---

Ridder's method, reduced interval may not translate to less iterations.
Ridder's method convergence rate may degenerate to bisection.

1/v = 1/c - 1/u ≤ 1/c - 1/a
v ≥ (a*c/(a-c))
(v² = b²-y) ≥ (a*c/(a-c))²

y ≤ (hi = b² - (a*c/(a-c))²)

max(c) = max(u*v/(u+v)) = a*b/(a+b)
min(c) = min(u*v/(u+v)) = u*0/(u+0) = 0

When y is small, (|u'|, |v'|) are relatively small, due to dominated constants, (a,b).
When y close to b², |v'| goes infinite, and c goes to 0 like fallen down a cliff.
This explains why f(y) look like a straight line, except when y very close to b²

We expect u*v/(u+v) concave-down, shaped like a quarter ellipse.

c/max(c) ≥ 1 - y/b²      // equality at the edges, y = 0 or b²

y ≥ (lo = b² * (1 - c*(a+b)/(a*b)))

Previously we simply use y = 0 .. b² = 900

lua> c = 15
lua> S.ridder(f, b*b*(1-c*(a+b)/(a*b)), b*b-(a*c/(a-c))^2, 1e-9, true)
112.5                  324
218.25                 255.6187346135275
236.93436730676376     255.6049211845503
246.26964424565705     255.60492082118668
250.93728253342186     255.6049208211843
255.6049208211839      255.6049208211843
255.60492082118404

lua> c = 1
lua> S.ridder(f, b*b*(1-c*(a+b)/(a*b)), b*b-(a*c/(a-c))^2, 1e-9, true)
847.5                  898.9480604865221
898.8848145203102      898.9480604865221
898.9199598288078      898.9480604865221
898.9199598288078      898.9199601069668
898.9199599678873      898.9199601069578
898.9199601069577      898.9199601069578
898.9199601069578

lua> c = 17
lua> S.ridder(f, b*b*(1-c*(a+b)/(a*b)), b*b-(a*c/(a-c))^2, 1e-9, true)
7.499999999999973      25.89792060491493
16.69896030245745      18.368146413229564
17.53355335784351      18.368139023879934
17.95084619086172      18.368139023879515
18.368139023879092     18.368139023879515
18.3681390238793

---

Reduced interval may not worth the trouble, unless we use it for secant's method.
I would ues 2 maximum guess, and iterate down, since f(y<0) is valid.

lua> S.secant(f, b*b, b*b-(a*c/(a-c))^2, 1e-9, true)

a/(a-c) = 1/(1-(c/a)) = 1 + (c/a) + (c/a)² + ... ≥ 1
We could use a simpler, but worse maximum guess.

lua> S.secant(f, b*b, b*b-c*c, 1e-9, true)

Or, we make f(y>b*b) vaild too, as long as we don't change the root.

lua> ssqrt = fn'x: x<0 and -sqrt(-x) or sqrt(x)'
lua> f = fn'y: local u,v = ssqrt(a*a-y), ssqrt(b*b-y); c*(u+v)-u*v'