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


18

If you're using Matlab, you may be interested in the Chebfun project. Chebfun takes a function, samples it and tries to represent it as a polynomial interpolant. If your function has discontinuities, Chebfun should be able to detect them with the splitting on command. You can find some examples here. If you're interested in the underlying algorithms, a good ...


14

For the first part of my question, I found this very useful comparison for performance of different linear interpolation methods using python libraries: http://nbviewer.ipython.org/github/pierre-haessig/stodynprog/blob/master/stodynprog/linear_interp_benchmark.ipynb Below is list of methods collected so far. Standart interpolation, structured grid: http:/...


14

Why not simply show the absolute value function? Approximation with e.g. Legendre-polynomial expansion works, but pretty badly: Taylor expansion is of course completely useless here, always giving only a linear function, either always decreasing or always increasing (depending on whether the point you expand around is negative or positive).


11

I suspect that chebfun algorithm must appear to be more practical, but it is neccesary to mention one more way to detect discontinuities, namely the discrete wavelet transform. You can get an idea how it works by looking at this Mathematica documentation page, see section > Applications > Detect Discontinuities and Edges. Briefly speaking, you can take DWT ...


11

If you need higher-order derivatives, you won't get good results using equidistant data points. If you can sample your function at arbitrary nodes, I would recommend using Chebyshev points, i.e. $$x_k = \cos\left( \pi \frac{k}{n} \right), \quad k=0\dots n$$ for a polynomial of degree $n$. You can evaluate the polynomial stably using Barycentric ...


11

Have you thought of using Barycentric Interpolation? A detailed description on how to do it efficiently for Chebyshev nodes is given in Section 5 of this paper. This is actually an exact evaluation of the Chebyshev interpolant. If you're evaluating a polynomial of degree $n$ at $m$ nodes, the cost is in $\mathcal O(nm)$. Update Another alternative, if you ...


11

Reducing the kernel width $\sigma_m$ will usually reduce the condition number. However, kernel matrices can become singular, or close to singular, for any basis function or point distribution, provided the basis functions overlap. The reason for this is actually quite simple: The kernel matrix $K$ is singular when its determinant $\det(K)$ is zero. ...


10

Given interpolation $I_H^h$ and restriction $I_h^H$ (where restriction is typically $(I_H^h)^T$ for symmetric problems), with fine grid discretized operator $A^h$, there are two common approaches for constructing the coarse grid operator $A^H$. (Petrov-)Galerkin coarse operators This explicitly computes the matrix triple product $$ A^H = I_h^H A^h I_H^h .$$...


10

Examining the documentation for the UnivariateSpline object, it seems as though you're constructing a 4th order spline, which is why you can't get values for derivatives of greater than 4th order. (All of those derivatives, as you're aware, would be zero.) SciPy limits the polynomial degree of splines to be 5th order or less (i.e., k <= 5). If you need ...


10

This is the classic colorization using optimization problem. Optimization and Linear Algebra To see how this can be expressed as a linear system, it's helpful to use a slightly different notation (and a slightly different objective function). Think of your image as graph $G$ with a node for each pixel in the image. There is an edge $(i,j)$ between two ...


10

It's a pathological case, but you can always resort to the Weierstrass monster function. It illustrates a broader point, namely that functions that are not smooth -- e.g., that have a kink -- are difficult to approximate because the interpolation error estimates require the function being interpolated to be differentiable a number of times. In other words, ...


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

I particularly like the bivariate spline class for what I think you are describing. You can use it to make a function (i.e. it is callable at any point) which interpolates the data using a spline. If you want just an interpolation then you simply set the kx and ky values to 1. If you want a smoother function then increasing the order of the spline (arguably ...


8

Weighted Essentially Non-Oscillatory (WENO) methods use "smoothness indicators" to detect discontinuities in finite volume and difference methods. From the description of Chebfun that Pedro gave, it seems as though the general idea is the same: construct a set of interpolating polynomials and use them to compute some measure of smoothness. See G.S. Jiang,...


8

You're asking how to produce dense output from your Runge-Kutta method. There are a number of ways to do this (see e.g. Hairer/Nørsett/Wanner). As noted in that reference, if you don't want to do more function evaluations aside from those already done by your fourth-order Runge-Kutta method, the best you can hope for is a third-order interpolant. This is ...


8

There are standard methods for these types of quadrature in Python, in NumPy and SciPy: Gauss-Laguerre quadrature Gauss-Legendre quadrature Gauss-Hermite quadrature (as noted in your post) Gauss-Chebyshev quadrature QUADPACK adaptive quadrature adaptive Gaussian quadrature and other routines. Given a degree of quadrature (essentially, the degree of the ...


8

There is no built-in Fortran functionality to do linear interpolation. You could either use a library or write your own routine. I haven't tried compiling or testing and my fortran may be a bit rusty, but something like the following should work. subroutine interp1( xData, yData, xVal, yVal ) ! Inputs: xData = a vector of the x-values of the data to be ...


8

This problem is generally called the Minimax Problem. Unfortunately the step function is not continuous and therefore the Weierstrass approximation theorem does not apply. Any continuous approximation will have $\epsilon \ge 0.5$ since there is a jump of size 1 so the best you can do at that point is split the difference. In fact, $y = 1/2$ is as good as you ...


8

That depends on how well you know the coordinates and velocities. If you have exact values, you can get a reasonable answer using Hermite interpolation. This will give you a degree-3 polynomial in each window that matches the coordinates and velocities at all the endpoints. Alternatively, if you do not know the coordinates and velocities exactly, but you do ...


7

Links to diverse software packages for scattered data interpolation are on my web page http://www.mat.univie.ac.at/~neum/stat.html#fit The book G.E. Fasshauer, Meshfree Approximation Methods using MATLAB, World Scientific 2007. gives a comprehensive state of the art (as of 2006). A few more recent papers on scattered data interpolation: http://www.stanford....


7

You can find a good overview of methods and vocabulary on interpolation in two dimensions at http://en.wikipedia.org/wiki/Multivariate_interpolation


7

If I understand you correctly, you want to fill in the values of the finer grid by interpolating on the coarser grid. One way to do linear interpolation on an unstructured grid is with Delaunay triangulations (this is how Matlab's griddata and TriScatteredInterp commands are implemented). After constructing a triangulation of your grid points, ...


7

Bernstein polynomials and Lagrange polynomials both span the same spaces. So in terms of the possible functions one can represent, using one or the other makes no difference. However, if you are thinking of using these as basis functions in either a finite element method or an interpolation problem, the spectral properties of the linear operator you create ...


7

I use Bernstein polynomials in a collocation method to solve boundary value problems for ODEs and PDEs. They are quite interesting. Convergence was exponential for some linear BVPs, but little slower compared to Chebyshev collocation, Legendre Galerkin, and Tau. Here's the figure comparing convergence rates with some Chebyshev spectral methods. The example ...


7

The RK4 method implicitly constructs a degree 3 polynomial interpolant, using the data $f(x_i)$, $f(x_{i+1})$, $f'(x_i)$, and $f'(x_{i+1})$ in each interval. This interpolant can be constructed rather easily and efficiently using a linear combination of shifted Hermite basis functions in each interval.


7

For a polynomial interpolant we have the formula for the pointwise error $$ E_n(x) = f(x) - p_n(x) = \frac{f^{(n+1)}}{(n+1)!}\prod_{i=0}^n(x-x_i),$$ where the $x_i$ are the interpolation knots. In general, we want to work with the right hand side of the error formula (a derivation appears in the book by Quarteroni et. al.) as the function we interpolate is ...


7

The frequency content of the interpolated signal is significantly influenced by the interpolation basis. If you have a band-limited function that you have adequately sampled (i.e. satisfying Nyquist criterion), interpolating with any function that is not band-limited to the same frequency will indeed introduce high-frequency noise. Unfortunately, exact band-...


7

Your choice of parameterization is creating problems. Instead of spanning one in $t$ between points, span an amount proportional to the line segment between the two points in $(x,y)$ space. I've created a Python example that demonstrates the issue. It compares a uniform parameterization $(x_t(t),y_t(t))$ (like yours, but mine is scaled so that $t \in [0,1]$...


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