26

The problem with equispaced points is that the interpolation error polynomial, i.e. $$ f(x) - P_n(x) = \frac{f^{(n+1)}(\xi)}{(n+1)!} \prod_{i=0}^n (x - x_i),\quad \xi\in[x_0,x_n] $$ behaves differently for different sets of nodes $x_i$. In the case of equispaced points, this polynomial blows up at the edges. If you use Gauss-Legendre points, the error ...


22

This is a really interesting question, and there are a lot of possible explanations. If we are attempting to use a polynomial interpolation, then note that polynomial satisfy the following annoying inequality Given a polynomial $P$ of degree not exceeding $N$ we have $$ |P^{\prime}(x)| \leq \frac{N}{\sqrt{1-x^2}}\max _x |P(x) | $$ for every $x \in (-1,1)$...


20

I would strongly advice against using closed form solutions since they tend to be numerically very unstable. You need to take extreme care in the way and order of your evaluations of the discriminant and other parameters. The classical example is the one for the quadratic equation $ax^2+bx+c=0$. Calculating the roots as $$x_{1,2} = \frac{-b \pm \sqrt{b^2-...


10

I posted some benchmarks here: http://www.cecm.sfu.ca/~rpearcea/mgb.html These are for total degree orders. To solve systems you typically need to do more work. Timings are for a typical midrange desktop as of 2015 (Haswell Core i5 quad core). The fastest system on one core is Magma, which uses floating point arithmetic and SSE/AVX. Magma is the ...


10

According to Wolfram Alpha, $x^5+3(x-1)=0$ has no closed-form solution, so you can forget about a nice closed-form expression. :) I see nothing wrong with Newton's method; it should be quick and accurate, and with some analysis like the one you are sketching I think you can identify a safe starting point and prove global convergence. It might even be faster ...


9

If you can stand to use complex arithmetic, simultaneous iteration methods might be preferable for computing all the roots of your polynomial. The simplest simultaneous iteration method, the (Weierstrass-)Durand-Kerner method, is effectively equivalent to applying Newton-Raphson to the Vieta relations relating the coefficients and roots of a polynomial, ...


9

Short answer You are missing the Jacobian of the transformation for the derivatives. Long answer The conditions that you propose for your interpolator translate into the following system of equations $$ \begin{bmatrix} 1 &x_1 &x_1^2 &x_1^3\\ 1 &x_2 &x_2^2 &x_2^3\\ 0 &1 &2x_1 &3x_1^2\\ 0 &1 &2x_2 &3x_2^2 \...


8

$\mathbf{A}$ is an $(n+1) \times (n+1)$ matrix. It can be obtained as follows: $\textbf{A} = \left[ \matrix{1 & 1 & 1 & \cdots & 1 \cr x_0 & x_1 & x_2 & \cdots & x_{n} \cr \vdots & \vdots & \vdots & \cdots & \vdots \cr x_0^n & x_1^n & x_2^n & \cdots & ...


8

Googling benchmark polynomial systems leads to a few hits, including the University of Mannheim's Computer Algebra Benchmark Initiative. Sadly, most of these are out of date or defunct. The most active seems to be the SymbolicData Wiki, but as far as I can tell, it only collects benchmark problems, not benchmark results. Some comparisons (dating back to ...


8

A possible solution (suggested by @gammatester) is to use Jacobi polynomials. This circumvents the problem of catastrophic cancellation in adding the large polynomial coefficients by 'naive' polynomial evaluation. The radial Zernike polynomial can be expressed by Jacobi polynomials as follows (see equation (6)) $$ R^m_n(\rho) = (-1)^{(n-m)/2}\rho^m \cdot P^...


7

As pointed out by David Ketcheson, one method is to use the companion matrix and find its eigenvalues (this is what Matlab does for the roots function). However, if you want to code everything by yourself, you can try to use the Sturm sequences, which you use as a first step to find an interval with only one zero. Then, you can apply one standard methods ...


7

Horner is indeed the most stable way to evaluate a polynomial (and you get the bonus of evaluating its derivatives with not too much extra cost). Higham presents a nice error analysis of the algorithm (Accuracy and stability of numerical algorithms, 2nd edition, p.94). He also presents an algorithm that includes a running error bound so you have an idea on ...


7

Note that, if $D$ is invertible, the eigenvalues of $A$ and $DAD^{-1}$ are the same. You can avoid floating-point underflow when forming the matrix by scaling the companion matrix by a diagonal matrix $D$ such that $D_{ii} = 1/\sqrt{(n-i)!}$. For the polynomial $x^4 + a x^3/\sqrt{1!} + b x^2/\sqrt{2!} + c x/\sqrt{3!} + d/\sqrt{4!}$, this gives a modified ...


7

In this paper, Honarvar and Paramesran derive an interesting method to compute the radial Zernike polynomials in a very nice recursive way. The recursion formula is surprisingly straightforward, without division or multiplication by large integers: $$ R^m_n(\rho) = \rho \left(R^{|m-1|}_{n-1}(\rho)+R^{m+1}_{n-1}(\rho)\right) - R^{m}_{n-2}(\rho)$$ I'd ...


6

This may seem like a somewhat exotic suggestion, but if your polynomials are of fixed, low degree, why not use the eigenvalues of the companion/colleague-matrix as an initial value for the roots? This may be expensive, but if you have to repeat the operation for slightly modified polynomials, you could use one or two steps of inverse iteration to update the ...


6

One way to get the $d_k$ is to expand each Chebyshev polynomial in powers, and then take linear combinations. However, using the Chebyshev representation itself (as detailed below) is numerically far more stable than using the power series expansion. You can evaluate a Chebyshev sum at any pareticular $x$ in a Horner-like fashion by using the 3-term ...


6

I don't know about the theoretical complexity. But the practically fastest algorithm for finding all zeros of a polynomial given as a finite power series is the algorithm by Traub [1]. It factors a degree $n$ polynomial in typically $O(n^2)$ time. See also [2]. For factoring a polynomial accessible only through its function values, see, e.g., my paper [3]....


6

You need to apply the Jacobian of the variable transformation. In other words, since you know how $x$ and $t$ are related, you can figure out how $dx$ and $dt$ are related. Obviously, they should differ by a factor of $(b-a)/2$ since that's the only substantive difference between the two expressions. In terms of more basic calculus, this is simply a $u$-...


6

Always keep in mind that the results of any benchmark will depend, in addition to the problem's size, on the base field over which the polynomial ring is defined (rational numbers or integers modulo some power of a prime number). The FGb library is an actively developed and high-performance implementation of the F5 algorithm. A benchmark comparing FGb to ...


6

Generally speaking, $P(A)$ would not be close to the zero matrix: if you compute the polynomial $P$ to a relative accuracy of $\epsilon$ (which is $10^{-16}$ in double-precision), then $P(A)$ would be something on the order of $\epsilon \|A^n\|$ ($n$—the size of the matrix), which would be huge. So you can only ask that $P(A)$ be as small as the typical ...


5

Division-free algorithms for the determinant are the Samuelson-Berkowitz algorithm (see german wikipedia) and the Leverrier-Faddeev algorithm (see german wikipedia). The latter requires division by integers. There is also a variant of the Gauß algorithm, the Gauß-Bareiss algorithm, that does book-keeping on numerators, i.e., keeping it in a factored form ...


5

I think you are looking for capabilities not available in Python itself nor Numpy, but that are available in SymPy. Sympy is a computer algebra package developed in Python. You can try it "live" in its on-line shell. Documentation on series development can be found here.


5

Certified homotopy continuation methods are used both for finding roots and for proving that they indeed exist (inside a certain interval). A quick web search turned out this paper: Reliable homotopy continuation by Joris van der Hoeven.


5

Unfortunately the 1D reasoning will not directly carry into 2D and beyond. Depending on the domain for your integral, you will have considerable latitude in the node choice - but of course the quadrature weights can be generated from polynomials as usual. The closest you could achieve to a classic Gaussian-type rule is to coordinate-transform your domain to ...


5

Yes, one can extend the Gaussian quadrature (in one dimension) to multiple dimensions. Often this is called cubature. Much work on this has been done by Ronald Cools and collaborators (see his Encyclopedia of Cubature Formulae). A nice paper for your application (integration on the unit disk) is A survey of known and new cubature formulas for the unit disk....


4

You could look into asymptotic behavior of discrete Dynamical Systems. There is both a rich theoretical literature on this topic in mathematics and more applied literature in physics and computer science.


4

You probably know this already, but it sounds like something from the world of chaos theory and fractals? (hence it is computationally "difficult") To your question, have you looked at the world of planetary mechanics and N-body problems? These are also forced to use iterative solutions, and the fundamental underlying physics is N^2, although the force ...


4

If you have a reasonable starting guess and the function you're trying to find a root of is not overly complex, then you can't beat Newton's method. It's also easy to implement -- a 10 line exercise in most languages.


4

As I was writing the question, I realized that in fact both of the given schemes achieve $O(m(d+1))$. This is because we always have $m \ge d$. In the expression tree scheme, the cost of quadratic polynomial multiplication at higher nodes in the tree can be amortized to lower nodes. In the polynomial interpolation scheme, at worst $h(d) = O(d^2)$ by ...


4

Suppose that you know the orthogonal polynomial basis in a single dimension $(x)$ of each degree $i$ up to some desired order $K$. That is, we know $$p_0(x),p_1(x),...p_i(x),...,p_K(x)$$ To extend this into a two dimensions $(x,y)$, we need only consider the product between 1D polynomials in (x) and (y) and collect only the products whose total degree is ...


Only top voted, non community-wiki answers of a minimum length are eligible