Tag Info

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)$...

11

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

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 ...

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

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

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

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

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

What you are doing here is projecting a function onto a basis of a finite dimensional space of polynomials. All you use in the formulas you show is that the basis $T_k$ is orthogonal, i.e., that $\left<T_i,T_j\right>=0$ if $i\neq j$ with regard to some scalar product $\left<\cdot,\cdot\right>$. All of this naturally generalizes to other ...

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....

5

OK, you have a very nice problem, I tried to run some benchmarks. First, I don't have your parameters so I used your small example. Second, since you do not specify the language, I used C + GSL (since I'm not that familiar with C++) #if __STDC_VERSION__ >= 199901L #define _XOPEN_SOURCE 700 #else #define _XOPEN_SOURCE 600 #endif /* __STDC_VERSION__ */ #...

5

The problem that you are seeing is a well known problem of these basis functions: $$\phi_{m,n}(x,y) = x^{m}y^{n}$$ I quote from here: The reason why the coefficient matrix is nearly singular and ill-conditioned is that our basis functions $\Psi_{i}(x) = x^{i}$ are nearly linearly dependent for large i. That is, $x^{i}$ and $x^{i+1}$ are very close for i ...

5

Quadrature points $x_1,x_2$ are also unknown. You need two more equations corresponding to $x^2$ and $x^3$. Then you solve for $a_1, a_2, x_1, x_2$.

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 ...

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

For this kind of large scale problem, using method like Gröbner basis, or other generally used to solve polynomial system, require lots of calculation time and many times your "solution" is so big you can't do anything with it. Also to solve the system, those method required that your system is zero-dimensional which is not all the time true. Probably a ...

4

You can easily do this with for loops. Just use n to define your loop limits. Here's a fully functional solution for MATLAB (just define n and x first): %pre-allocate A A = zeros(n+1); %first row: for j=0:n A(1,j+1)=sum(x.^j); end %rows 2 through n for i=1:n A(i+1,1:n)=A(i,2:n+1); %copy from previous row A(i+1,n+1)=sum(x.^(n+i)); %compute last ...

4

Let's use the example (approximate the square root function on $[0.25,1.0]$ with a quartic polynomial) to step through your calculations. I suspect that the code is going to work with only modest changes. Having chosen an initial set of six interpolation points: $$x_k = 0.25, 0.4, 0.55, 0.7, 0.85, 1.0\;\; (k=0,\ldots,5)$$ we proceed to interpolate the ...

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