# Tag Info

11

I posted some benchmarks here: http://www.cecm.sfu.ca/~rpearcea/mgb.html (archived copy) 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 ...

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

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

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

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

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

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

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

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

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

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

Philips in 1988 proved the following relationship: If $f(x)$ is an infinitely differentiable function defined on the interval $[-1,+1]$ and its Legendre expansion is given by $f(x) = \sum_{n=0}^{\infty}a_{n}P_{n}(x)$, then the Legendre coefficients $a_{n}^{(q)}$ of the $q$-th derivative of $f(x)$ are given by $$a_{n}^{(q)} = \frac{(2n+1)}{2^{q-2}(q-1)!} \... 4 Yes. You can compute a running error bound, i.e, a number \mu such that the difference between the exact value of y = p(x) and the computed value satisfies \hat{y} satisfies$$|y - \hat{y}| \leq \mu u.$$Here u is the unit roundoff. You can trust the sign computed sign of y, when |\hat{y}| > \mu u. Let p(x) = \sum_{j=0}^n a_j x^j, then ... 4 One technique is to use a library with arbitrary large integers and use fixed point arithmetic. Just scale up x by some integer factor, then compute your polynomial in fixed-point, keeping only the sign. From two values, that you know positive and negative, you can then perform a binary search. At the very end, divide back by the fixed point factor and get ... 4 Taking your questions in order, the definition of D_r as the operator converting the nodal values of u_h into derivatives follows immediately from the definition, although the notation is a little inexact,:$$ u_h(r) = \sum_i u(r_i)l_i(r)$$hence$$\frac{du_h}{dr} = \sum_j u(r_j)\frac{dl_j}{dr} $$so that$$\left.\frac{du_h}{dr}\right|_{r_i} = \sum_j u(...

3

What you have is most likely a mesh for a Spectral Element Method or something similar). The solution is represented using a high order polynomial each section of the cube (referred to as an element). For convenience, these polynomials are represented as Lagrange nodal functions, which are associated with a set of points on each element (which are denoted by ...

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