Programming Challenge: a classic trigonometry problem
|
03-10-2023, 03:53 PM
(This post was last modified: 03-16-2023 11:11 AM by Albert Chan.)
Post: #35
|
|||
|
|||
RE: Programming Challenge: a classic trigonometry problem
(03-08-2023 04:09 PM)Albert Chan Wrote: (1) and (2): (1+1/p)² * (1-p²) = (a²-b²)/c² = k .... (3) 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 |
|||
« Next Oldest | Next Newest »
|
User(s) browsing this thread: 1 Guest(s)