![]() |
[QUOTE=Stargate38;375241]By exact, I mean something like this:
[TEX]\sum\limits_{i=a}^\infty f(i)[/TEX] where f(i) is a function of i such that the above sum converges to the root. If you need to, use nested sums.[/QUOTE] In what manner is that "more exact" (or even "more concise") than giving an iterative formula which encodes a sequence of iterates that converges to the root? The only difference I see is that you appear to desire an explicit expression for your summand, whereas the iterative-formula approach defines its argument inductively. |
Should you not be asking this question on Stack Exchange?
|
[QUOTE=ewmayer;375255]In what manner is that "more exact" (or even "more concise") than giving an iterative formula which encodes a sequence of iterates that converges to the root? The only difference I see is that you appear to desire an explicit expression for your summand, whereas the iterative-formula approach defines its argument inductively.[/QUOTE]
It is not "exact" in any meaningful sense, especially since ALL real-analytic functions can be put in this form. It is also meaningless without saying anything about f(x). One could define f() to be "the smallest real root of the following septic polynomial.........." |
[QUOTE=R.D. Silverman;375287]It is not "exact" in any meaningful sense, especially since ALL
real-analytic functions can be put in this form. It is also meaningless without saying anything about f(x). One could define f() to be "the smallest real root of the following septic polynomial.........."[/QUOTE] Exactly - pardon the pun. It is an interesting epistemological question to try to delineate the (fuzzy) boundary between exact and approximate here. For instance, one could require that any candidate for exactness should encode a constructive method for producing the answer to any desired precision, thus the self-referential quote you give above would be ruled out on those grounds. ================ I spent a bit of time considering higher-order analogs of the N-R scheme which will not only converge quickly from a good initial guess at the root, but can also cope with the deliberately-problematic initial iterate x0 = 0, at which f = -1 and f`-f```` all = 0. Define a set of (k)th-order-convergent iterative schemes x[sub]n+1[/sub] = x[sub]n[/sub] - f*G[sub]k-1[/sub](x)/G[sub]k[/sub](x) , [1] where the G[sub]k[/sub](x) are functions-of-the-function-f(x) defined via the recurrence G[sub]k[/sub](x) = f`*G[sub]k[/sub] - f*G`[sub]k[/sub]/(k-1) , starting with G[sub]1[/sub] = 1 . [2] We can see that G[sub]2[/sub] = f`, which is just the 2nd-order N-R scheme. k = 3 gives G[sub]3[/sub](x) = f`[sup]2[/sup] - f*f``/2, yielding the 3rd-order Halley scheme, &c. The lowest-order iteration safe to start at x0 = 0 needs a 5th derivative (f*f`````) or higher term in denominator, which means the 6th-order scheme of the family, but that gives dx = 0 since the numerator = 0 at x = 0, thus we need thw 7th-order scheme. I'll spare you the gory algebraic details. but we end up with the following 2 G-terms in the 7th-order version of [1]: G[sub]6[/sub] = ((f`)^2 - 2.f.f``).(f`)^3 + f^2.f`.(2.f`.f``` + 3.(f``)^2)/4 - f^3.((f`.f```` + 2.f``.f```) + f.f`````/10)/12 G[sub]7[/sub] = (f`)^4.((f`)^2 - (5/2).f.f``) + f^2.(f`)^2.(4.f`.f``` + 9.(f``)^2)/6 - f^3.(f``.((f``)^2 + 4.f`.f```) + (f`)^2.f````)/8 + f^4.(12.f`.f````` + 30.f``.f```` + 20.(f```)^2 - f.f``````)/720 Now try the resulting iterative scheme starting from x0 = 0; here is simple Pari code: default(realprecision, 200) x = 0; [i][Repeat the following sequence][/i] f = 243*x^7 + 16*x^6 + x^5 - 1; f1 = 1701*x^6 + 96*x^5 + 5*x^4; f2 = 1701*6*x^5 + 480*x^4 + 20*x^3; f3 = 1701*30*x^4 + 1920*x^3 + 60*x^2; f4 = 1701*120*x^3 + 5760*x^2 + 120*x; f5 = 1701*360*x^2 + 11520*x + 120; f6 = 1701*720*x + 11520; print(f) G_6 = (f1^2 - 2*f*f2)*f1^3 + f^2*f1*(2*f1*f3 + 3*f2^2)/4 - f^3*((f1*f4 + 2*f2*f3) - f*f5/10)/12; G_7 = f1^4*(f1^2 - 5*f*f2/2) + f^2*f1^2*(4*f1*f3 + 9*f2^2)/6 - f^3*(f2*(f2^2 + 4*f1*f3) + f1^2*f4)/8 + f^4*(12*f1*f5 + 30*f2*f4 + 20*f3^2 - f*f6)/720; dx = -f*G_6/G_7; x += dx; x Note that because our function is a polynomial one could further substitute the derivatives f-f6 into the expressions for G[sub]6[/sub] and G[sub]7[/sub] and collect common powers of x in each, but I'm only interested in checking the convergence properties, so didn't bother doing that. Starting from x0 = 0 we obtain the following sequence of iterates (and f(x) value at each) - Notice that the approximant barely budges away from x = 0 on the first iteration, but that is enough for the above functions to start getting a nonzero signal for the intermediate derivatives, so on step 2 we get more than halfway to the root, and it's all downhill from there: [code]j f(x) x_j -- ------------ -------------------------------- 0 -1 1 -0.999997... 0.0625 2 -0.955180... 0.28249... 3 -1.75..E-002 0.44507... 4 4.73..E-018 0.4462235223189588700... [good to ~19 digits] 5 -5.10..E-127 0.4462235223189588697797792058806101653432216192296867823821087097882260386219976960158801016445364065226392916613728558306678486... [good to 127 digits][/code] Now from a computational-efficiency perspective this is of course the wrong way to go - The problem at x = 0 here illustrates the first thing one should do whenever deploying such a scheme, namely to produce a bracketing interval a < x < b containing one or more roots (1 is the preferred count). For the OP's 7th-order poly it is easy to verify that 0 < x < 1 is a bracketing interval since f(0) < 0 and f(1) > 0. The initial guess for the root is drawn from this interval - typically x = (a+b)/2 is as good a starting value as any, and indeed x0 = 1/2 works perfectly well as a starting value for the basic 2nd-order N-R scheme here. During the ensuing iteration loop if the updated iterate x lands outside (a,b) one should try to narrow the bracketing interval and retry, and/or retry with a higher-order iterative scheme. |
In the special case of polynomial roots, one could also directly compute approximations to the numerical values via a variety of methods: computing the eigenvalues of a matrix derived from the polynomial coefficients, or using sophisticated algorithms like Jenkins-Traub that are nearly bulletproof. Then use your favorite iterative method to 'polish' the computed roots. You get more robustness than when using iterative methods alone, and only need high-precision arithmetic after the bulletproof methods finish.
|
[QUOTE=jasonp;375377]In the special case of polynomial roots, one could also directly compute approximations to the numerical values via a variety of methods: computing the eigenvalues of a matrix derived from the polynomial coefficients, or using sophisticated algorithms like Jenkins-Traub that are nearly bulletproof. Then use your favorite iterative method to 'polish' the computed roots. You get more robustness than when using iterative methods alone, and only need high-precision arithmetic after the bulletproof methods finish.[/QUOTE]Another approach, hinted at in Knuth, is to use exact rational arithmetic. That is treat a floating point number as p/(2^n) where p and n are integers. This approach can give greatly reduced issues with rounding and cancellation errors together with the easily availability and simplicity of use of multiprecision arithmetic compared with the typical floating point counterpart.
|
[QUOTE=xilman;375409]Another approach, hinted at in Knuth, is to use exact rational arithmetic. That is treat a floating point number as p/(2^n) where p and n are integers. This approach can give greatly reduced issues with rounding and cancellation errors together with the easily availability and simplicity of use of multiprecision arithmetic compared with the typical floating point counterpart.[/QUOTE]
Yep! Why not?? Let's solve it over the 2-adics, then lift the result to R. :smile: :smile::smile::smile: for the humor impaired.// |
[QUOTE=R.D. Silverman;375426]Yep! Why not??
Let's solve it over the 2-adics, then lift the result to R. :smile: :smile::smile::smile: for the humor impaired.//[/QUOTE]Nice to see someone is paying attention and spots a gentle troll. |
Did that mean "An approximation is an approximation, in any base" ?
(Otherwise, wouldn't computing in multiprecision integer plus "large" binary part (say > 10^5 bits) produce as valid a "solution" ? ) |
Actually, it is an interesting exercise to perform an iterative approximation of one's choice over Q in a package like Pari which auto-reduces the resulting quotients to LCD form and compare how fast the "bitness" of the approximants grows relative to (error vs exact) compared to a multiprecision-float approach.
|
Post #48 helps, but how do I see the whole number in Pari? I keep seeing "[+++]" and the number is cut off.
|
| All times are UTC. The time now is 23:28. |
Powered by vBulletin® Version 3.8.11
Copyright ©2000 - 2021, Jelsoft Enterprises Ltd.