The inverted integral of LinearInterpolation

[SouthEndMusic]

What kind of problems is it mostly used for? Please describe.

For my application I need to compute the water level in a bucket given a piecewise linear area(level) relationship and a volume:

It is assumed that the area is always positive.

Describe the algorithm you’d like

volume(level) relationship:

$$ V = V_i + \int_{\ell_i}^{\ell} A_i + \frac{\Delta A_i}{\Delta \ell_i}(\ell-\ell_i)\text{d}t = \frac{1}{2}\frac{\Delta A_i}{\Delta \ell_i}(\ell-\ell_i)^2 + A_i(\ell-\ell_i) + V_i $$

So this quadratic equation must be solved for $\ell$:

$$ \ell = \ell_i + \Delta \ell_i\frac{-A_i + \text{sign}(\Delta A_i)\sqrt{A_i^2 + 2\frac{\Delta A_i}{\Delta \ell_i}(V-V_i)}}{\Delta A_i} $$

This does lead to numerical problems when $\frac{\Delta A_i}{\Delta \ell_i}$ is small, so I could use some help with that.

Other implementations to know about

I implemented this as CLinearInterpolationInvInv here: https://github.com/SouthEndMusic/CachedInterpolations.jl.

Is this too niche for DataInterpolations.jl?

[ChrisRackauckas]

Is this doable for all interpolations?

[SouthEndMusic]

In favour of doability:

• We could cache the cumulative integral at the data points, which could be used to speed up the existing integral calculations and can serve as the t for the integral inverse

Against doability:

• First of all: we need a not too expensive way to determine whether the integral is invertible, for which we can check whether the original interpolation is strictly positive everywhere. This is trivial for linear interpolation, but not for more complicated interpolation types.
• Then there is the matter of computing the inverse analytically, which generally means solving a polynomial equation. This already gave me some problems with edge cases in the linear interpolation case as mentioned above, which only gets worse with more complicated methods. Or we have to use a numerical solver for this, but that is probably more expensive to compute

[ChrisRackauckas]

It could be added to the interface with a trait for which algorithms support it, and just not supported by most algorithms to start. That could be fine.

[SouthEndMusic]

I think the following expression is sufficiently stable (where $s$ is the slope of the linear interpolation):

$$ \ell = \ell_i + \frac{2(V - V_i)}{A_i + \sqrt{A_i^2 + 2s(V - V_i)}}, $$

as we require $A_i > 0$, and for $s = 0$ it nicely reduces to the linear

$$ \ell = \ell_i + \frac{V - V_i}{A_i}. $$