# Looking for a quadrature rule

**URL:** <https://discourse.julialang.org/t/looking-for-a-quadrature-rule/111923>\
**Category:** Numerics\
**Created:** [March 21, 2024, 8:47am UTC](https://discourse.julialang.org/t/looking-for-a-quadrature-rule/111923 "2024-03-21T08:47:02Z")\
**Posts on this page:** 18\
**Page:** 1

<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:** [March 21, 2024, 8:47am UTC](https://discourse.julialang.org/t/looking-for-a-quadrature-rule/111923/1 "2024-03-21T08:47:02Z")

</div>

Hi,

I need compute thousands of integrals of the form

\int\_{0}^{t} \frac{e^{at}}{\sum\_{i=1}^n b\_i e^{c\_i t}} dt.

I know that :

- the time frame [0,t] is typically small (say t \le 1)
- rates a,c\_1,...c\_n are of order 1e-4, 1e-5
- weights b\_1,...,b\_n are of order 1.

For the moment, I use the _very dumb but fast_ forward scheme

\int\_0^t f(s)ds = t \* f(0).

Is there something smarted that i could do without blowing up the computational cost ?

---

<div class="post-metadata">

**Author:** ![IlianPihlajamaa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ilianpihlajamaa/32/28766_2.png) [@IlianPihlajamaa](https://discourse.julialang.org/u/IlianPihlajamaa)\
**Post date:** [March 21, 2024, 8:50am UTC](https://discourse.julialang.org/t/looking-for-a-quadrature-rule/111923/2 "2024-03-21T08:50:40Z")

</div>

The expression in your post you can integrate exactly. I expect there is a missing t in the denominator?

---

<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:** [March 21, 2024, 8:51am UTC](https://discourse.julialang.org/t/looking-for-a-quadrature-rule/111923/3 "2024-03-21T08:51:36Z")

</div>

Yep there was, just fixed it. I am having trouble showing the integral centered i do not know why.

---

<div class="post-metadata">

**Author:** ![IlianPihlajamaa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ilianpihlajamaa/32/28766_2.png) [@IlianPihlajamaa](https://discourse.julialang.org/u/IlianPihlajamaa)\
**Post date:** [March 21, 2024, 8:56am UTC](https://discourse.julialang.org/t/looking-for-a-quadrature-rule/111923/4 "2024-03-21T08:56:03Z")

</div>

How fast and accurate do you want your solution to be? Would a low order [Gauss-Legendre](https://juliaapproximation.github.io/FastGaussQuadrature.jl/dev/gaussquadrature/) rule work?

---

<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:** [March 21, 2024, 8:58am UTC](https://discourse.julialang.org/t/looking-for-a-quadrature-rule/111923/5 "2024-03-21T08:58:53Z")

</div>

Maybe, thanks for the link I’ll try it.

I do not need too much accuracy, I think that an error of order 0.001 would be largely good enough. I am more looking around that certain of what i want here…

---

<div class="post-metadata">

**Author:** ![IlianPihlajamaa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ilianpihlajamaa/32/28766_2.png) [@IlianPihlajamaa](https://discourse.julialang.org/u/IlianPihlajamaa)\
**Post date:** [March 21, 2024, 9:01am UTC](https://discourse.julialang.org/t/looking-for-a-quadrature-rule/111923/6 "2024-03-21T09:01:31Z")

</div>

If the rates are low and the time domain is small, you can also do a Taylor expansion of the exponentials, and integrate the result exactly.

---

<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:** [March 21, 2024, 9:03am UTC](https://discourse.julialang.org/t/looking-for-a-quadrature-rule/111923/7 "2024-03-21T09:03:04Z")

</div>

Taylor expensions of the exponentials would give me a rational function right ? Rational functions are integrated exactly ?

---

<div class="post-metadata">

**Author:** ![IlianPihlajamaa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ilianpihlajamaa/32/28766_2.png) [@IlianPihlajamaa](https://discourse.julialang.org/u/IlianPihlajamaa)\
**Post date:** [March 21, 2024, 9:10am UTC](https://discourse.julialang.org/t/looking-for-a-quadrature-rule/111923/8 "2024-03-21T09:10:21Z")

</div>

Yes. So if I’m not mistaken, using a first order expansion, the integrand  
becomes

I = \int\_0^t\frac{1+at'}{B+Dt'}dt'

where B=\sum\_i b\_i and D = \sum\_i b\_i c\_i.  
This can be integrated exactly, yielding

I = \frac{{D} \left(a t+\log \left(\frac{{D} t}{B}+1\right)\right)+a B \log \left(\frac{B}{B+\text{D} t}\right)}{{D}^2}

---

<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:** [March 21, 2024, 9:17am UTC](https://discourse.julialang.org/t/looking-for-a-quadrature-rule/111923/9 "2024-03-21T09:17:12Z")

</div>

Wow. Indeed this looks a bit stronger than my zeroth order approx. Could the same thing be done for second order too ?

Edit: WolframAlpha gives

 ![image](https://global.discourse-cdn.com/julialang/original/3X/b/7/b76846ec7405c8a799dbfa1cf4e130e07c0ef081.png)

and

 ![image](https://global.discourse-cdn.com/julialang/original/3X/3/4/34a62a9fd5f7f85e4019b9e1cd224c28cc809a54.png)

But this is starting to be a bit messy ^^

I’ll try to implement the first order one to see how it behaves. If it is not sucessfull enough, I’ll try a low-order quadrature as you pointed out. Thanks 🙂

---

<div class="post-metadata">

**Author:** ![IlianPihlajamaa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ilianpihlajamaa/32/28766_2.png) [@IlianPihlajamaa](https://discourse.julialang.org/u/IlianPihlajamaa)\
**Post date:** [March 21, 2024, 9:26am UTC](https://discourse.julialang.org/t/looking-for-a-quadrature-rule/111923/10 "2024-03-21T09:26:26Z")

</div>

EDIT: ah I see you beat me to it. 🙂

Yes, but the integral becomes more hairy. Mathematica says  
I\approx \frac{\frac{2 \left(a^2 \left(D^2-B E\right)-a D E+E^2\right) \tan ^{-1}\left(\frac{D+E \text{t}}{\sqrt{2 B E-D^2}}\right)}{\sqrt{2 B E-D^2}}-\frac{2 \left(a^2 \left(D^2-B E\right)-a D E+E^2\right) \tan ^{-1}\left(\frac{D}{\sqrt{2 B E-D^2}}\right)}{\sqrt{2 B E-D^2}}+a^2 E \text{t}+a (E-a D) \log \left(2 B+2 D \text{t}+E \text{t}^2\right)+a \log (B) (a D-E)+a \log (2) (a D-E)}{E^2},  
where E = \sum\_i b\_ic\_i^2.

you would probably be better off doing a Taylor expansion of the whole integrand at this point.

---

<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:** [March 21, 2024, 9:30am UTC](https://discourse.julialang.org/t/looking-for-a-quadrature-rule/111923/11 "2024-03-21T09:30:03Z")

</div>

Maybe indeed I should taylor expand the whole integrand. Could your Mathematica do this for me ? Wolframalpha does not want too 😛

Or maybe I should go with Symbolics.jl + TaylorDiff.jl

---

<div class="post-metadata">

**Author:** ![vavrines](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vavrines/32/8024_2.png) [@vavrines](https://discourse.julialang.org/u/vavrines)\
**Post date:** [March 21, 2024, 9:32am UTC](https://discourse.julialang.org/t/looking-for-a-quadrature-rule/111923/12 "2024-03-21T09:32:39Z")

</div>

So numerically we are solving the problem with the form f\_t=I(f,t) within t\in [0,t']?  
Then you may check [OrdinaryDiffEq.jl: ODE solvers and utilities · OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/).  
In terms of efficiency, to achieve the same accuracy, higher-order integrators are usually cheaper than the forward Euler scheme.

---

<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:** [March 21, 2024, 9:34am UTC](https://discourse.julialang.org/t/looking-for-a-quadrature-rule/111923/13 "2024-03-21T09:34:45Z")

</div>

Yes but the integrand is piecewise defined and has a different shape on every time step, which is why i ended up making the scheme manually.

I am more looking for a _better approximation_ than t \* f(0) than for a _good approximation_: I do not want the runtime to explode.

---

<div class="post-metadata">

**Author:** ![IlianPihlajamaa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ilianpihlajamaa/32/28766_2.png) [@IlianPihlajamaa](https://discourse.julialang.org/u/IlianPihlajamaa)\
**Post date:** [March 21, 2024, 9:38am UTC](https://discourse.julialang.org/t/looking-for-a-quadrature-rule/111923/14 "2024-03-21T09:38:00Z")

</div>

Symbolics.jl can do all of this as well, but here you go

 ![image](https://global.discourse-cdn.com/julialang/original/3X/8/0/80da90b550c2b4636789f72ca9ff041d6499a08b.png)

---

<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:** [March 21, 2024, 9:44am UTC](https://discourse.julialang.org/t/looking-for-a-quadrature-rule/111923/15 "2024-03-21T09:44:22Z")

</div>

Using ForwardDiff, I could directly do it with the following:

\int\_0^t f(s) ds = \sum\_{k=0}^\infty f^{(k)}(0) \frac{t^{k+1}}{(k+1)!} \approx tf(0) + \frac{t^2}{2} f^{(1)}(0)

which would probably be numerically efficient, and gives me one more order than the previous tf(0) approximation.

---

<div class="post-metadata">

**Author:** ![IlianPihlajamaa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ilianpihlajamaa/32/28766_2.png) [@IlianPihlajamaa](https://discourse.julialang.org/u/IlianPihlajamaa)\
**Post date:** [March 21, 2024, 9:45am UTC](https://discourse.julialang.org/t/looking-for-a-quadrature-rule/111923/16 "2024-03-21T09:45:21Z")

</div>

Note that if you are using this as a way of solving a differential equation, the scheme may become unstable (just like forward Euler does). If so, trying a midpoint or backward rule may work better.

---

<div class="post-metadata">

**Author:** ![Sushrut\_Deshpande](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sushrut_deshpande/32/52754_2.png) [@Sushrut\_Deshpande](https://discourse.julialang.org/u/Sushrut_Deshpande)\
**Post date:** [March 21, 2024, 10:57am UTC](https://discourse.julialang.org/t/looking-for-a-quadrature-rule/111923/17 "2024-03-21T10:57:40Z")

</div>

Hello,  
I think most of the computation time for a such a problem lies in the denominator of the function. If a possible approximation works for the denominator doing thousands of integrals should not be a problem. But I do not know the degree of accuracy needed.

---

<div class="post-metadata">

**Author:** ![PeterSimon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petersimon/32/25193_2.png) [@PeterSimon](https://discourse.julialang.org/u/PeterSimon)\
**Post date:** [March 22, 2024, 4:52pm UTC](https://discourse.julialang.org/t/looking-for-a-quadrature-rule/111923/18 "2024-03-22T16:52:48Z")

</div>

Here is an approximation tested with some typical values of `a`, `b`, and `c`:

```julia
using QuadGK, BenchmarkTools

a = 0.0003
b = [1.0, 1.5, 2.3, 0.6]
c = 1.e-4 * b

function exact_integrand(t, a, b, c)
    return exp(a*t) / sum(b[i] * exp(c[i]*t) for i in eachindex(b,c))
end

approx_exp(x) = 1 + x 

function approx_integrand(t, a, b, c)
    approx_exp(a*t) / sum(b[i] * approx_exp(c[i]*t) for i in eachindex(b,c))
end

```

If we plot the exact and approximate integrands we see the following:  
 ![integrand](https://global.discourse-cdn.com/julialang/original/3X/3/9/3914b43889f2ee4bc64c642e211bffcfcaf519fb.png)  
The approximation looks very good and both traces appear to be simple linear functions. So we can try a one-point quadrature rule that samples the integrand at the center of the integration interval:

```julia
ts = 0:0.01:1
ints_exact = [quadgk(t -> exact_integrand(t,a,b,c), 0, t)[1] for t in ts]
ints_approx = [t * approx_integrand(t/2,a,b,c) for t in ts]
@show maximum(abs, (ints_exact - ints_approx))

```

which produces

```julia
maximum(abs, ints_exact - ints_approx) = 1.4760460631535466e-9

```

The timing for the approximate and exact integrands:

```julia-repl
julia> @btime approx_integrand(0.95, $a, $b, $c)
  7.800 ns (0 allocations: 0 bytes)
0.1852089640409081

julia> @btime exact_integrand(0.95, $a, $b, $c)
  31.194 ns (0 allocations: 0 bytes)
0.1852089689678891

```