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>
- Re: How to eliminate noises? A better way perhaps.
- References:
- How to eliminate noises?
- From: Artur <grafix@csl.pl>
- How to eliminate noises?