0
$\begingroup$

I am given with a function $f(x)$ on a grid $X_{old}\in \{x_{min},x_{max}\}$ with a non uniform spacing. I need to interpolate that function on a new log-spaced grid $X_{new}\in\{x^{\prime}_{min},x^{\prime}_{max}\}$ with same numbers of discretized point, $\mathcal{N}$, such that the following constraint is satisfied, $$\int_{X_{old}}f_{old}(x)\,dx=\int_{X_{new}}f_{new}(x)\,dx$$.

Is there any algorithm for doing such kind of interpolation. I tried to find online but could not succeed. Therefore any help would be highly appreciated.

$\endgroup$
6
  • $\begingroup$ The $f$ in the two integrals are not the same. You might want to use an index on $f$ to make clear which mesh it is defined on. $\endgroup$ Apr 7, 2023 at 21:40
  • $\begingroup$ What do we mean by the integral? The algorithm with depend on what numerical procedure is used for the integration. But once a definition is made, .e.g., the Simpson rule, then one can probably come up with a conservative interpolation algorithm. $\endgroup$ Apr 7, 2023 at 23:02
  • $\begingroup$ @MaximUmansky, thanks for the comment. I am using trapezoidal rule for computing the integral. $\endgroup$
    – Sayan
    Apr 8, 2023 at 3:50
  • $\begingroup$ @WolfgangBangerth, sure. Now modified. Thanks $\endgroup$
    – Sayan
    Apr 8, 2023 at 3:51
  • $\begingroup$ @Sayan Can we make an assumption that the end points for the two grids are the same? Otherwise, if we are allowed to arbitrarily change the limits of integration, but the integral must stay the same, it is hard to imagine something sensible can be done there. $\endgroup$ Apr 8, 2023 at 4:19

2 Answers 2

2
$\begingroup$

I'll add some thoughts and terminology.

  1. First, your question doesn't make that much sense or is a bit underspecified. You are given a function values on a grid $\{x_i\}$, and you want the interpolation to a new grid $\{x^\prime_i\}$. However, this is usually not what interpolation is usually for, because interpolation takes a input a function and not function values.

  2. So the first thing you should do is to interpolate the original function values on the original grid to obtain the base function $f(x)$. For example, you can use Lagrange interpolation for that, $f(x) = \sum_i f_i L_i(x)$, where $L_i(x) = \prod_{i\neq j} \frac{x-x_i}{x_j - x_i}$.

  3. Next you'll need the integrated value, let's call that $I = \int_{x_0}^{x_{N+1}} f(x) dx$. This is usually evaluated by some sort of quadrature on the original gridpoints/function values.

  4. Now you can face the actual problem, which is then basically a resampling: $f^\prime(x^\prime) = \sum_i f(x^\prime_i) L^\prime(x^\prime)$. The result is the identical function as $f(x)$, but sampled on a new grid. It's just that you define the constructing polynomials on the new grid. Note tat this resampling might be very badly conditioned -- for example, try to map the Chebyshav grid with $N=128$ to an equally spaced grid ... this won't work due to an exponentially bad conditioning.

As long as you are using a quadrature rule that integrates the Lagrange-representation exactly, or close to machine precision, the integral of the resampled function should be identical -- as it's the same polynomial.

In any other case, the new representation will give a somewhat different integral. You can then try to solve the least squares problem and include the constraint that the integral is constant. This will then yield, however, $N+1$ conditions, which can't be solved exactly but only approximately.

$\endgroup$
5
  • $\begingroup$ Instead of using least squares, you may rather shift the polynomial functions by the difference of both integrals. $\endgroup$
    – ConvexHull
    Apr 8, 2023 at 19:24
  • 1
    $\begingroup$ @ConvexHull: ok, that's clever and in fact a more easy approach, but could give suboptimal results as it prioritizes the integral-condition over the interpolation condition of all other points. Still, it's a bit random to reason without seeing this in action. $\endgroup$
    – davidhigh
    Apr 8, 2023 at 20:05
  • $\begingroup$ Maybe. The new nodal values will definitely not match the old polynomial. This is generally not possible if the polynomials are different, as @wolfgangbangerth already said. $\endgroup$
    – ConvexHull
    Apr 9, 2023 at 3:23
  • $\begingroup$ But then it's not interpolation. $\endgroup$ Apr 10, 2023 at 2:32
  • $\begingroup$ @Bangerth: one could reduce the number of interpolated nodes to $N-1$ and apply the quadrature constraint as the $N$th. Not interpolation either, but probably suitable for the current use case. $\endgroup$
    – davidhigh
    Apr 10, 2023 at 18:03
1
$\begingroup$

You can't do this with interpolation. Imagine for example that you are representing the function $f(x)=x^2$ on the interval $[-1,1]$ on a mesh $X_{old}$ with a very small mesh width. This representation will be rather accurate due to the small mesh, and consequently $$ \int_{-1}^{1} f_{old}(x) dx \approx \int_{-1}^{1} f(x) dx = \frac 13 x^3|_{-1}^1 = \frac 23. $$ But then re-interpolate $f_{old}$ on a new mesh that has only one cell $[-1,1]$. The function is defined by two node points at $\pm 1$ where $f_{old}=1$ and so $f_{new}(x)=1$, with $$ \int_{-1}^{1} f_{new}(x) dx =2. $$ The point is that interpolation simply doesn't preserve the mean value.

$\endgroup$
2
  • $\begingroup$ One extra point I forgot to add in the question (now added) is that both the meshes $X_{old}$ and $X_{new}$ are discetized with same numbers of points (128 typically for my case). $\endgroup$
    – Sayan
    Apr 8, 2023 at 3:50
  • $\begingroup$ @Sayan But you get the idea. You can do the same kind of argument having all 128 points shifted somewhere close to the left end of the interval in the one case, and shifted close to the right hand in the other case. The point is simply that interpolation does not preserve the mean value of a function -- that's just how it is. $\endgroup$ Apr 10, 2023 at 2:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.