# Higher accuracy in numerical cumulative integration

**URL:** <https://discourse.julialang.org/t/higher-accuracy-in-numerical-cumulative-integration/96898>\
**Category:** Numerics\
**Tags:** question, numerics, approxfun\
**Created:** [March 31, 2023, 3:08pm UTC](https://discourse.julialang.org/t/higher-accuracy-in-numerical-cumulative-integration/96898 "2023-03-31T15:08:35Z")\
**Posts on this page:** 16\
**Page:** 1

<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:08pm UTC](https://discourse.julialang.org/t/higher-accuracy-in-numerical-cumulative-integration/96898/1 "2023-03-31T15:08:35Z")

</div>

There are some threads about numerical cumulative integration like [this](https://discourse.julialang.org/t/how-do-i-do-a-fast-cumulative-integration/85574/) and [this](https://discourse.julialang.org/t/solve-for-upper-lower-bound-of-a-numerical-definite-integral/58164). However, I have a more basic question.

I am interested in performing the integral S(x; a) = \int\_a^x f(x') dx' where x \le b. Note that f(x) is always positive so S(x;a) is monotonically increasing. f(x) has the form f(x) = \sqrt{a\_1 - a\_2 x^4 + a\_3 x^2}.

After I have S(x; a), I want to find x(S; a), namely, I want to invert it.  
I understand that I can do the integral numerically for many values of x, and then solve the equation S(x; a) = x\_0 to build x(S; a).

Another complication is that the integration limits a, b actually are a function of some other parameter c, but I am guaranteed that S(b(c); a(c)) = \int\_{a(c)}^{b(c)} f(x') dx' is always the same value. I mean the total integral is a constant.  
That means that I have a bunch of functions S(x; c) that are monotonically increasing from 0 to some fixed value S\_0 = S(b(c); a(c)) .  
And now I want to invert it.

I am not sure what is the best way to get high accuracy.  
I tried the idea by @stevengj that I read [here](https://discourse.julialang.org/t/how-do-i-do-a-fast-cumulative-integration/85574/6) and that is implemented [here](https://discourse.julialang.org/t/how-do-i-do-a-fast-cumulative-integration/85574/5) by @tobydriscoll. I mean, I used `ApproxFun` to approximate the function f(x), and then to integrate it.

```julia
a = 1
b = 2
f_approx = Fun(x->f(x), a..b)
Stmp = integrate(f_approx)
S = Stmp - Stmp(a)

```

Then I used Newton’s Method to invert the function S(x; a) to x(S;a).  
The problem is that for some values of c, I find that the total integral S(b(c); a(c)) = \int\_{a(c)}^{b(c)} f(x') dx' gives a smaller value than my initially calculated S\_0. This means that I cannot find the root.

Are there any ideas on how to improve this calculation?  
Maybe I should use `quadgk` to calculate the integral instead of `ApproxFun`?

---

<div class="post-metadata">

**Author:** ![JM\_Beckers](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jm_beckers/32/22482_2.png) [@JM\_Beckers](https://discourse.julialang.org/u/JM_Beckers)\
**Post date:** [March 31, 2023, 3:42pm UTC](https://discourse.julialang.org/t/higher-accuracy-in-numerical-cumulative-integration/96898/2 "2023-03-31T15:42:58Z")

</div>

So if you know S\_0, why not solve

{dx \over dS}={1 \over f(x)} ;

x(0)=a 

by integrating between 0 and S\_0 and check that at the endpoint you have x(S\_0)=b. If not, increase the accuracy requirement of the integrator.  
Maybe not the most efficient way, but you could in this way enforce your integral constraint and have the solution x(S) .

---

<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, 3:51pm UTC](https://discourse.julialang.org/t/higher-accuracy-in-numerical-cumulative-integration/96898/3 "2023-03-31T15:51:13Z")

</div>

> [@roi.holtzman](#):
>
> f(x) has the form f(x) = \sqrt{a\_1 - a\_2 x^4 + a^3 x^2}

A basic difficulty here may be that your function has a square root singularity, which may make it harder to numerically integrate accurately unless you build this singularity into the quadrature scheme. Does the argument of your square root ever go through zero in your integration interval? (Of course, that would make the integrand complex.)

---

<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, 4:00pm UTC](https://discourse.julialang.org/t/higher-accuracy-in-numerical-cumulative-integration/96898/4 "2023-03-31T16:00:27Z")

</div>

> [@stevengj](#):
>
> A basic difficulty here may be that your function has a square root singularity

Yes, I see that is a difficulty.

> [@stevengj](#):
>
> Does the argument of your square root ever go through zero in your integration interval?

The argument of the sqrt stays positive or zero always. Though I am interested in values of x that make f(x) vanish. I guess that this is a problem, but I would like to still have better accuracy.

---

<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, 4:02pm UTC](https://discourse.julialang.org/t/higher-accuracy-in-numerical-cumulative-integration/96898/5 "2023-03-31T16:02:40Z")

</div>

> [@JM\_Beckers](#):
>
> So if you know S\_0, why not solve

That’s an interesting suggestion! I will try to implement that.

---

<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, 4:07pm UTC](https://discourse.julialang.org/t/higher-accuracy-in-numerical-cumulative-integration/96898/6 "2023-03-31T16:07:58Z")

</div>

> [@roi.holtzman](#):
>
> The argument of the sqrt stays positive or zero always.

Places where it hits zero are still a problem for convergence rates. At the very least, you should break the integral into subintervals where such roots occur. Then an adaptive scheme like QuadGK will have an easier time refining the integral near the singularity. (And I would tend to favor an h-adaptive scheme like QuadGK over a p-adaptive scheme like ApproxFun unless you build the nature of the singularity into the function approximations, e.g. by using Gauss–Jacobi quadrature as I suggested [in another thread](https://discourse.julialang.org/t/domain-error-with-quadgk-to-derive-caputo-derivative/96550/2).)

---

<div class="post-metadata">

**Author:** ![JM\_Beckers](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jm_beckers/32/22482_2.png) [@JM\_Beckers](https://discourse.julialang.org/u/JM_Beckers)\
**Post date:** [March 31, 2023, 4:32pm UTC](https://discourse.julialang.org/t/higher-accuracy-in-numerical-cumulative-integration/96898/7 "2023-03-31T16:32:00Z")

</div>

I think in the OPs case (function f stays positive or zero), the zero can only happen at x=a or x=b, otherwise you would have negative arguments for the square root within the domain (I assume a\_2 is positive) ?

---

<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, 4:33pm UTC](https://discourse.julialang.org/t/higher-accuracy-in-numerical-cumulative-integration/96898/8 "2023-03-31T16:33:31Z")

</div>

> [@JM\_Beckers](#):
>
> the zero can only happen at x=ax=ax=a or x=bx=bx=b

True!

---

<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, 4:38pm UTC](https://discourse.julialang.org/t/higher-accuracy-in-numerical-cumulative-integration/96898/9 "2023-03-31T16:38:34Z")

</div>

Thanks!!

So in my case, due to this singularity, I want to use QuadGK and not the ApproxFun method.

> [@stevengj](#):
>
> unless you build the nature of the singularity into the function approximations, e.g. by using Gauss–Jacobi quadrature as I suggested

So I read your post there – it was very good. Unfortunately, for me, I did not see that my type of singularity

> [@roi.holtzman](#):
>
> f(x) has the form f(x) = \sqrt{a\_1 - a\_2 x^4 + a\_3 x^2} .

is included in these methods, due to the powers of x.

So I guess that the best option is to just use `quadgk` trying the best accuracy I can get.

---

<div class="post-metadata">

**Author:** ![JM\_Beckers](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jm_beckers/32/22482_2.png) [@JM\_Beckers](https://discourse.julialang.org/u/JM_Beckers)\
**Post date:** [March 31, 2023, 4:46pm UTC](https://discourse.julialang.org/t/higher-accuracy-in-numerical-cumulative-integration/96898/10 "2023-03-31T16:46:08Z")

</div>

You might also try symbolic integration, Mathematica finds the expression for the indefinite integral and maybe with your specific constants you can do the work analytically.

---

<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, 4:54pm UTC](https://discourse.julialang.org/t/higher-accuracy-in-numerical-cumulative-integration/96898/11 "2023-03-31T16:54:59Z")

</div>

> [@roi.holtzman](#):
>
> Unfortunately, for me, I did not see that my type of singularity is included in these methods, due to the powers of x.

It is included in these methods. What matters is not the exact form of the integrand, but the form of the singularity. In particular, suppose that the argument of the square root goes to zero at a the lower endpoint x=a. In this case, the integrand is of the form  
\sqrt{p(x) (x - a)} = g(x) \sqrt{x - a} where p(x) is a (cubic) polynomial) and g(x) = \sqrt{p(x)} is non-singular at the endpoint. That is, the singularity is only in the \sqrt{x - a} weight, which is handled by Gauss–Jacobi quadrature.

---

<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, 7:26pm UTC](https://discourse.julialang.org/t/higher-accuracy-in-numerical-cumulative-integration/96898/12 "2023-03-31T19:26:40Z")

</div>

> [@JM\_Beckers](#):
>
> You might also try symbolic integration, Mathematica finds the expression for the indefinite integral and maybe with your specific constants you can do the work analytically.

I tried that. Mathematica indeed gives an analytic expression (with elliptic integral functions), but for some reason, it comes out complex (even when I put all the proper assumptions). Then, in the singular points, namely x\_0 such that f(x\_0)=0, I get that the integral is evaluated to be infinity. So I need to understand first why this comes out as infinity.

---

<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, 7:32pm UTC](https://discourse.julialang.org/t/higher-accuracy-in-numerical-cumulative-integration/96898/13 "2023-03-31T19:32:29Z")

</div>

@stevengj Perfect explanation! thanks!  
So I guess that this means I need to partition my segment to two to handle the right and left sides separately.

I will try that.

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [April 1, 2023, 7:03pm UTC](https://discourse.julialang.org/t/higher-accuracy-in-numerical-cumulative-integration/96898/14 "2023-04-01T19:03:21Z")

</div>

It’s really easy in ApproxFun to factor out a sqrt singularity using `JacobiWeight`

---

<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:** [April 1, 2023, 7:28pm UTC](https://discourse.julialang.org/t/higher-accuracy-in-numerical-cumulative-integration/96898/15 "2023-04-01T19:28:13Z")

</div>

Thanks @dlfivefifty !  
Could you maybe point me to an example that shows how to use the weights in ApproxFun?

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [April 1, 2023, 7:35pm UTC](https://discourse.julialang.org/t/higher-accuracy-in-numerical-cumulative-integration/96898/16 "2023-04-01T19:35:42Z")

</div>

```julia
julia> f = Fun(x -> sqrt(x+1)*cos(x), JacobiWeight(0.5,0));

julia> F = cumsum(f)
Fun((1+x)^1.5[Chebyshev()], [0.527618, 0.104918, -0.0649745, -0.00159742, 0.000894048, 1.02892e-5, -5.55621e-6, -3.72868e-8, 1.97678e-8, 8.70474e-11, -4.56323e-11, -1.42217e-13, 7.40716e-14])

julia> F(0.1)
0.6957793344826042

```