Programming Challenge: a classic trigonometry problem
|
03-07-2023, 06:43 PM
(This post was last modified: 03-18-2023 11:17 AM by Albert Chan.)
Post: #23
|
|||
|
|||
RE: Programming Challenge: a classic trigonometry problem
(01-01-2014 08:44 AM)Thomas Klemm Wrote: I used additional variables: 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' |
|||
« Next Oldest | Next Newest »
|
User(s) browsing this thread: 4 Guest(s)