Statalist


[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: st: mata polyroots issue


From   Hua Peng <hpeng@stata.com>
To   statalist@hsphsun2.harvard.edu
Subject   Re: st: mata polyroots issue
Date   Mon, 11 May 2009 12:39:15 -0500

Lirac<lirac271@yahoo.com> asked a question about mata function polyroot():

>I am running Stata 10.1 on Mac OS 10.4.  The following issue is perplexing me:
>
>-------
>mata:
>mata clear
>polynomial1 = (1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-30,29)
>polyeval(polynomial1,1)
>polyroots(polynomial1)
>-------
>
>The polyeval line makes it clear that x=1 is a root for this polynomial.
>However, 1 is not a reported result for polyroots.

The first two results reported by polyroots(polynomial1) are
(1-1.1729e-08i, 1+1.1729e-08i). The results are correct.

Notice polynomial1 can be factored into

	(x-1)*(29*x^29 - x^28 - x^27 - x^26 -... - 1),

and the second term is 0 at x=1. Hence polynomial has a factor (x-1)^2,
i.e., 1 is a double root of polynomial1. Mata correctly identifies this.

Due to numerical round-off, Mata computes the degree two factor as
(x-1)^2+e, where e is a small positive number, instead of
exact (x-1)^2, and consequently Mata reports two complex conjugate roots
close to 1 instead of two 1.

>What is odd is that
>if I add or subtract a zero term to the polynomial, polyroots() will
>then correctly report x=1 as a root.  i.e., the following works:
>
>polynomial2 = (1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-30,29)
>polynomial3 = (1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-30,>29)
>
>Does anybody know what is going on here?
>

Both polynomial2 and polynomial3 have only (x-1) factor.  Mata correctly
reports 1 as a simple root.

--Hua
(hpeng@stata.com)
*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/



© Copyright 1996–2014 StataCorp LP   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index