[Date Index]
[Thread Index]
[Author Index]
Re: Approximate/asymptotic factorization
*To*: mathgroup at smc.vnet.net
*Subject*: [mg73630] Re: [mg73509] Approximate/asymptotic factorization
*From*: Daniel Lichtblau <danl at wolfram.com>
*Date*: Fri, 23 Feb 2007 04:40:18 -0500 (EST)
*References*: <200702210635.BAA17453@smc.vnet.net>
Paul Abbott wrote:
> Consider the monomial H[4,4]
>
> 256*(23625 + 126000*z^4 - 7200*z^8 + 3840*z^12 + 256*z^16)
>
> where H[m,n] is a Generalized Hermite Polynomial (Noumi & Yamada 1998),
> a rational solution to the Painlevï¿½ equation P_IV. A plot of the roots
> of H[4,4] displays an interesting quasi-rectangular structure.
>
> Then consider the monomial Q[4,4]
>
> 1528899609315374375625 + 6795109374734997225000*z^4 -
> 560866170613047390000*z^8 + 153399294526645440000*z^12 +
> 2734590598399296000*z^16 - 167891551278796800*z^20 +
> 2948686820352000*z^24 - 40649991782400*z^28 + 277762867200*z^32 -
> 920125440*z^36 + 1048576*z^40
>
> where Q[m,n] is a Generalized Okamoto Polynomial, another rational
> solution to P_IV. A plot of the roots displays an even more interesting
> structure, with a quasi-rectangular core and 4 triangular "lobes"
> (related, in some way, to the the Yablonskiiï¿½Vorobï¿½ev Polynomials, Q[n],
> which are Rational Solutions of P_II)
>
> If you overlay plots of the roots of H[4,4] and Q[4,4] what is most
> striking is that the "common" roots are very similar. In other words,
> H[4,4] is "approximately" a factor of Q[4,4].
>
> Is there a "standard" way of quantifying this approximate factorization?
>
> This approximate factorization holds for any m and n and appears to
> improve as m and n increases. Perhaps there is a limit in which the
> factorization holds asymptotically?
>
> Cheers,
> Paul
>
> _______________________________________________________________________
> Paul Abbott Phone: 61 8 6488 2734
> School of Physics, M013 Fax: +61 8 6488 1014
> The University of Western Australia (CRICOS Provider No 00126G)
> AUSTRALIA http://physics.uwa.edu.au/~paul
I'm not sure how to tackle the asymptotics you have in mind. For
particular examples this is very much the province of univariate
approximate gcd computation. There has been a good deal of work in this
area since around 1995, with some work even earlier.
A good place to start would be the "4M" article of Corless, Gianni,
Trager, and Watt, "The singular value decomposition for polynomial
systems", ISSAC 1995. The idea of their approach is to use the SVD of
the Sylvester matrix in order to determine the degree of a "good" approx
gcd, then derive that gcd also using the SVD.
For your example I show some quick-and-dirty code that, I think, does
more or less what they have in mind. For the Sylvester matrix I cribbed
and slightly modified code from MathWorld:
http://mathworld.wolfram.com/SylvesterMatrix.html
SylvesterMatrix[poly1_, poly2_, var_] := Module[
{coeffs1, coeffs2, l1, l2},
coeffs1 = Reverse[CoefficientList[poly1, var]];
coeffs2 = Reverse[CoefficientList[poly2, var]];
l1 = Length[coeffs1];
l2 = Length[coeffs2];
Join[Table[
Join[Table[0,{j}], coeffs1,
Table[0,{l2-2-j}]], {j,0,l2-2}],
Table[Join[Table[0,{j}], coeffs2,
Table[0,{l1-2-j}]], {j,0,l1-2}]]
]
h44 = 256*(23625 + 126000*z^4 - 7200*z^8 + 3840*z^12 + 256*z^16);
q44 = 1528899609315374375625 + 6795109374734997225000*z^4 -
560866170613047390000*z^8 + 153399294526645440000*z^12 +
2734590598399296000*z^16 - 167891551278796800*z^20 +
2948686820352000*z^24 - 40649991782400*z^28 + 277762867200*z^32 -
920125440*z^36 + 1048576*z^40;
I'll normalize so largest coefficient is unity.
normPoly[poly_, var_] := Expand[poly /
Max[Abs[CoefficientList[poly,var]]]]
p1 = normPoly[h44,z];
p2 = normPoly[q44,z];
I'll use an epsilon of 10^(-3) as tolerance for what singular values to
remove.
eps = 10^(-3);
We count the number of singular values that were zero.
nulls = Select[S,Max[Abs[#]]<eps&];
We now get the degree of our approximate gcd as below (for the epsilon
we chose, it is 12).
k = Length[sylmat] - Length[nulls]
To form that gcd we'll row reduce a certain matrix and use the final row.
newmat = Take[newS.Transpose[V], k];
coeffrow = Last[RowReduce[newmat]];
gcdapprox = coeffrow . z^Reverse[Range[0,Length[coeffrow]-1]]
So here is what I get.
In[150]:= InputForm[gcdapprox]
Out[150]//InputForm=
6.980499621431014 + 0.000585041745605589*z + 0.0035109293315991083*z^2 +
0.002602929711965681*z^3 + 33.10928622044708*z^4 +
0.0019204205746833373*z^5 + 0.01152459048089943*z^6 +
0.008544087788492293*z^7 - 3.49084586854123*z^8 +
0.00011025782928825875*z^9 + 0.0006624721373616564*z^10 +
0.0004911544753497421*z^11 + z^12
I'll note that an epsilon of around 10^(-1) would have given a degree 16
approx gcd (not equal to a scalar multiple of p1, though). One can get
an idea of how close are the various roots by looking at them as below.
s1 = Sort[z /. NSolve[p1==0, z]];
s2 = Sort[z /. NSolve[p2==0, z]];
sg = Sort[z /. NSolve[gcdapprox==0, z]];
Playing tricks with Intersection and a comparison criterion might
automate this process.
Note that there are other approaches to approximate univariate gcds
using QR decomposition, remainder sequences, approximate factoring, and
probably some others I am forgetting.
Daniel Lichtblau
Wolfram Research
Prev by Date:
**FullSimplify question**
Next by Date:
**Re: split**
Previous by thread:
**Approximate/asymptotic factorization**
Next by thread:
**Re: Approximate/asymptotic factorization**
| |