Re: How to eliminate noises? A better way perhaps.

*To*: mathgroup at smc.vnet.net*Subject*: [mg122734] Re: How to eliminate noises? A better way perhaps.*From*: Noqsi <noqsiaerospace at gmail.com>*Date*: Tue, 8 Nov 2011 07:15:44 -0500 (EST)*Delivered-to*: l-mathgroup@mail-archive0.wolfram.com*References*: <201111021121.GAA03503@smc.vnet.net> <j8tl5t$f3t$1@smc.vnet.net> <j90gmp$shi$1@smc.vnet.net>

On Nov 4, 4:01 am, Richard Fateman <fate... at cs.berkeley.edu> wrote: > On 11/3/2011 1:59 AM, Bob Hanlon wrote:> Use higher precision. > > .. so say you all... > Well, that's the brute force method, and one reason to recommend it is > it doesn't require much thought. > > On the other hand, the person asking the question has an apparently > large expression which he knows has a double zero at n=1/2 and he wants > to know the behavior of the expression from 0.35 to 0.53. Namely, > around that zero. > > Let us call that expression p. A non-brute force, but mostly automatic > method is to note that (since p is a polynomial of degree 29) it is > EXACTLY equal to s, where > > s = Normal[Series[p, {n, 1/2, 29}]] > > A plot of s can be done in ordinary float arithmetic and looks to the > eye just like the plot of p, done with high working precision. Unique to > computer algebra systems, it is also possible to look at the structure > of s, and note that expanded around n=1/2 it is mighty sparse. > > Let k = n-1/2 . The expression being plotted is then > > 14 k^2 - 1296 k^5 + 1451188224 k^11 - 2002639749120 k^14 + > 598075300380672 k^17 - 83414577743659008 k^20 + > 3977900382788517888 k^23 - 113721152119718805504 k^26 + > 1516282028262917406720 k^29 > > Sometimes a better way to format expressions for evaluation is to use > HornerForm. In this case it is > > k^2 (14 + > k^3 (-1296 + .... > > Evaluating this expression at a value for k requires computing > k^2,k^3,k^6,and some other arithmetic for a total of about 18 arithmetic > operations, which can be done in double-float. > The real win here is that you can show how much more insight you might > get from using computer algebra. I find it a great pleasure to be in agreement with Richard for a change. But let me offer another perspective, less insightful perhaps, but more susceptible to automation. First, let me note the nice behavior of N[Root[]]. For this problem, it is useful to have all of the roots in numerical form, and I think the following approach is the best way to get them in Mathematica: nRoots[p_, x_] := With[{f = p /. x -> #}, Table[N[Root[f, n]], {n, 1, Exponent[p, x]}]] Now, note the the documentation for Root makes the following promise: "The ordering used by Root[f,k] takes real roots to come before complex ones, and takes complex conjugate pairs of roots to be adjacent. " One difficulty with this promise is that it doesn't tell you how to find the break between the real part of the vector and the complex part. The following code assumes that N[Root[f,n]] will reliably have head Real for real roots. Oh Wolfram experts, is this true? If you want to blindly make a numerical polynomial for evaluation in floating point, it is useful to put it in factored form. Factored form isn't as quick to evaluate as "streamlined" forms like Horner, but it is more stable. The roots are the raw material you need to create this form. The real roots correspond to linear terms. The following function constructs them and collects their product: linearTerms[r_, x_] := Times @@ (x - # & /@ Cases[r, _Real]) The complex root pairs correspond to quadratic terms. There are various ways to construct these. This function constructs the form that Hart et al. ("Computer Approximations", Kreiger 1978, p. 69) recommend for stability: quadraticTerm[r_, x_] := With[{z = r[[1]]}, (x - Re[z])^2 + Im[z]^2] This collects the product of those quadratic terms: quadraticTerms[r_, x_] := Times @@ (quadraticTerm[#, x] & /@ Partition[Cases[r, _Complex], 2]) nFactor[ p, x ] yields a factored numerical representation of polynomial p in variable x: nFactor[p_, x_] := With[{r = nRoots[p, x]}, N[Coefficient[p, x, Exponent[p, x]]] linearTerms[r, x] quadraticTerms[r, x]] Now test it: nfn = nFactor[(-1 + 2 n)^2 (-5077311495040 + 256771815840909 n - 6285087811295946 n^2 + 99108751736880828 n^3 - 1130572436430176064 n^4 + 9932051583279890496 n^5 - 69845429986320295296 n^6 + 403489321408322272512 n^7 - 1949815250878086761472 n^8 + 7985109816575584026624 n^9 - 27976513439837284724736 n^10 + 84423827219928959287296 n^11 - 220448430557212109488128 n^12 + 499528909151513065046016 n^13 - 983522463480512590086144 n^14 + 1682263859668487045185536 n^15 - 2495534893540205641334784 n^16 + 3200240859272788875411456 n^17 - 3529950418080000008257536 n^18 + 3325177849291516859645952 n^19 - 2648782194444972790382592 n^20 + 1760454182292395183308800 n^21 - 958285503480825479430144 n^22 + 416134125894081039040512 n^23 - 138626084433937223909376 n^24 + 33263436995017750609920 n^25 - 5117451845387346247680 n^26 + 379070507065729351680 n^27), n] Plot[nfn, {n, 0.35, 0.53}] This produces the now-familiar plot. Both this approach and Richard's yield an expression that can be used for a variety of purposes, including export via CForm[] or FortranForm[]. Richard's has the advantage that in this case it yields more human insight, and a more efficient expression. The disadvantage of Richard's approach is that it *requires* human insight, while this approach may be applied in a blind, automated fashion. John Doty Noqsi Aerospace, Ltd. http://www.noqsi.com/

**Follow-Ups**:**Re: How to eliminate noises? A better way perhaps.***From:*Andrzej Kozlowski <akoz@mimuw.edu.pl>

**References**:**How to eliminate noises?***From:*Artur <grafix@csl.pl>