# Approximate inverse of a definite integral

**URL:** <https://discourse.julialang.org/t/approximate-inverse-of-a-definite-integral/83181>\
**Category:** Numerics\
**Tags:** question, integral\
**Created:** [June 22, 2022, 12:21pm UTC](https://discourse.julialang.org/t/approximate-inverse-of-a-definite-integral/83181 "2022-06-22T12:21:43Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [June 22, 2022, 12:21pm UTC](https://discourse.julialang.org/t/approximate-inverse-of-a-definite-integral/83181/1 "2022-06-22T12:21:44Z")

</div>

For a given (fixed) function f: [0,1] \to [0, \infty), let

g(x) = \int\_0^x f(x) dx

and

h(x) = \frac{g(x)}{g(1)}

h is then, by construction, [0,1] \to [0,1]. I need a fast approximation of its _inverse_: I am OK with investing a bit in constructing a good approximation to h^{-1}, because it gets called tens of thousands of times. It is “smooth”: no kinks, discontinuities. Orthogonal polynomials would do fine (eg Chebyshev).

The problem is calculating it smartly — eg if I have a grid of Chebyshev nodes mapped to [0,1]. I thought of formulating the integral as an ODE, and then using a [continuous callback](https://diffeq.sciml.ai/stable/features/callback_functions/#ContinuousCallback) when it crosses a node to save the position. But that requires a first pass to calculate g(1).

Maybe I am overthinking this. Any suggestion is appreciated.

---

<div class="post-metadata">

**Author:** ![tbeason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tbeason/32/15898_2.png) [@tbeason](https://discourse.julialang.org/u/tbeason)\
**Post date:** [June 22, 2022, 12:34pm UTC](https://discourse.julialang.org/t/approximate-inverse-of-a-definite-integral/83181/2 "2022-06-22T12:34:46Z")

</div>

My first thought is that if you are indeed calling h or h^{-1} “tens of thousands” of times then computing g(1) _once_ should not be problematic.

In the same way you might use Chebyshev regression to approximate a function, now approximate the inverse instead. So compute `(x,h(x))` where `x` is “dense enough” in [0,1] but then interpolate as `(h(x),x)` – so that the `x` values are now the LHS of your interpolant.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [June 22, 2022, 1:45pm UTC](https://discourse.julialang.org/t/approximate-inverse-of-a-definite-integral/83181/3 "2022-06-22T13:45:10Z")

</div>

I would just use Newton’s method here to compute h^{-1} from a quadrature routine for h(x). Something like:

```julia
using QuadGK
const g₁ = quadgk(f, 0, 1; rtol=1e-12)[1] # compute g(1) accurately
g(x; atol=1e-9 * g₁) = quadgk(f, 0, x; atol=atol)[1]
h(x; atol=1e-9) = g(x; atol = atol * g₁) / g₁
h′(x) = f(x) / g₁ # the derivative, using the fundamental theorem of calculus
function h⁻¹(y; atol=1e-9)
    x = 0.5 # arbitrary initial guess
    while true
        δx = (h(x; atol) - y) / h′(x)
        x -= δx
        abs(δx) ≤ atol && return x
    end
end

```

which gives, e.g. for `f(x) = x^3 + 3`:

```julia
julia> y = h(0.3)
0.2775461538461539

julia> h⁻¹(y)
0.30000000000000004

```

You can then plug `h⁻¹` into e.g. [FastChebInterp](https://github.com/stevengj/FastChebInterp.jl) or [ApproxFun](https://github.com/JuliaApproximation/ApproxFun.jl) to form an efficient Chebyshev interpolant rather than applying Newton’s method every time you want to evaluate it.

(If f is extremely expensive, you might want to first replace it by a Chebyshev interpolant too, before using the interpolant to construct h and h^{-1}. Of course, once you have the interpolant, you don’t need QuadGK, e.g. with ApproxFun using `cumsum` will form the integral function directly. You still probably want Newton’s method to invert it.)

> [@Tamas\_Papp](#):
>
> I thought of formulating the integral as an ODE

In general, you almost never want to evaluate definite integrals by solving ODEs — specialized integration routines are _much_ more efficient (and simpler to use), because they exploit the fact that you have “random access” to your integrand values at any point in the domain (rather than having to build it up from “left to right” as in a general ODE).

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [June 22, 2022, 2:22pm UTC](https://discourse.julialang.org/t/approximate-inverse-of-a-definite-integral/83181/4 "2022-06-22T14:22:19Z")

</div>

If you “precompute” stuff only once, you might want to precompute them in `BigFloat` format to achieve maximum accuracy in later use of the approximator.

---

<div class="post-metadata">

**Author:** ![roi.holtzman](https://avatars.discourse-cdn.com/v4/letter/r/f05b48/32.png) [@roi.holtzman](https://discourse.julialang.org/u/roi.holtzman)\
**Post date:** [March 28, 2023, 10:47pm UTC](https://discourse.julialang.org/t/approximate-inverse-of-a-definite-integral/83181/5 "2023-03-28T22:47:18Z")

</div>

@Tamas_Papp What did you end up doing? I am curious as I am thinking of the same thing.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [March 29, 2023, 8:21am UTC](https://discourse.julialang.org/t/approximate-inverse-of-a-definite-integral/83181/6 "2023-03-29T08:21:42Z")

</div>

I stored the approximate of the inverse using Chebyshev polynomials, used it as a starting point, and then did 2 Newton refinement steps. This gave me machine precision for my problem. YMMV.

---

<div class="post-metadata">

**Author:** ![roi.holtzman](https://avatars.discourse-cdn.com/v4/letter/r/f05b48/32.png) [@roi.holtzman](https://discourse.julialang.org/u/roi.holtzman)\
**Post date:** [March 31, 2023, 2:01pm UTC](https://discourse.julialang.org/t/approximate-inverse-of-a-definite-integral/83181/7 "2023-03-31T14:01:15Z")

</div>

I am still not sure what is the best (i.e. maximal accuracy) way to do that.  
In this [post](https://discourse.julialang.org/t/how-do-i-do-a-fast-cumulative-integration/85574/6) @stevengj explains why it is better to approximate the function you want to integrate (if you are interested in a cumulative integration).

So I am trying that (following the example in the same [thread](https://discourse.julialang.org/t/how-do-i-do-a-fast-cumulative-integration/85574/5)).  
I am also interested in the inverse function h^{-1}. So I use Newton’s Method to invert it for any value that I am interested in. I am not happy with the accuracy I get. Part of the complication, in my case, is that the function f(x) is actually a function that depends on another parameter f(x, a), and then for different a's I do not get consistent results.

I think the main problem comes from the approximation near the boundaries of the integral \int\_{x\_1}^{x\_2} f(x') dx'. Is there any way to increase the approximation of [ApproxFun.jl](https://juliaapproximation.github.io/ApproxFun.jl/latest/), or to account for the boundaries x\_1, x\_2 better?

Alternatively, is there a better way (more accurate) to compute this integral and its inverse?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 31, 2023, 2:56pm UTC](https://discourse.julialang.org/t/approximate-inverse-of-a-definite-integral/83181/8 "2023-03-31T14:56:43Z")

</div>

> [@roi.holtzman](#):
>
> Alternatively, is there a better way (more accurate) to compute this integral and its invers

Really hard to say without knowing what integral you are computing, and how you are computing it now.

---

<div class="post-metadata">

**Author:** ![roi.holtzman](https://avatars.discourse-cdn.com/v4/letter/r/f05b48/32.png) [@roi.holtzman](https://discourse.julialang.org/u/roi.holtzman)\
**Post date:** [March 31, 2023, 3:09pm UTC](https://discourse.julialang.org/t/approximate-inverse-of-a-definite-integral/83181/9 "2023-03-31T15:09:41Z")

</div>

Thanks for the answer! I just posted a new [question](https://discourse.julialang.org/t/higher-accuracy-in-numerical-cumulative-integration/96898) with some more details as I didn’t want to clutter this thread. So the details about the function are there.