# Quadratures and cubatures with Dirac deltas

**URL:** <https://discourse.julialang.org/t/quadratures-and-cubatures-with-dirac-deltas/108683>\
**Category:** Numerics\
**Tags:** integral, integrate, quadrature\
**Created:** [January 11, 2024, 1:20pm UTC](https://discourse.julialang.org/t/quadratures-and-cubatures-with-dirac-deltas/108683 "2024-01-11T13:20:10Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![homocomputeris](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/homocomputeris/32/8933_2.png) [@homocomputeris](https://discourse.julialang.org/u/homocomputeris)\
**Post date:** [January 11, 2024, 1:20pm UTC](https://discourse.julialang.org/t/quadratures-and-cubatures-with-dirac-deltas/108683/1 "2024-01-11T13:20:10Z")

</div>

Are there ways to compute cubatures of functions that include the Dirac delta function?

---

<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:** [January 11, 2024, 2:03pm UTC](https://discourse.julialang.org/t/quadratures-and-cubatures-with-dirac-deltas/108683/2 "2024-01-11T14:03:30Z")

</div>

How are your functions represented? How would you be able to detect a delta function if you can only evaluate pointwise?

---

<div class="post-metadata">

**Author:** ![tom-plaa](https://avatars.discourse-cdn.com/v4/letter/t/d26b3c/32.png) [@tom-plaa](https://discourse.julialang.org/u/tom-plaa)\
**Post date:** [January 11, 2024, 2:55pm UTC](https://discourse.julialang.org/t/quadratures-and-cubatures-with-dirac-deltas/108683/3 "2024-01-11T14:55:11Z")

</div>

I have done this in the past by rewriting the target integral in a way that separates the dirac delta from the rest, analytically computing the integral w.r.t. the delta (it’s just the integrand’s value at the point) and then evaluating the leftover integral. I did this in a case where I had a probability density that had a dirac delta at zero (for a stochastic jump process).

As @dlfivefifty implies, you need to know a priori where you’re going to get the delta. If you do, just define the rule to evaluate it and then compute the non-delta components numerically.

If you have a concrete example, we can try to help you think about the problem more precisely.

---

<div class="post-metadata">

**Author:** ![homocomputeris](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/homocomputeris/32/8933_2.png) [@homocomputeris](https://discourse.julialang.org/u/homocomputeris)\
**Post date:** [January 11, 2024, 4:47pm UTC](https://discourse.julialang.org/t/quadratures-and-cubatures-with-dirac-deltas/108683/4 "2024-01-11T16:47:28Z")

</div>

Yeah, my bad. I know where the deltas are. Can I make cubatures step into these points somehow?

---

<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:** [January 11, 2024, 5:18pm UTC](https://discourse.julialang.org/t/quadratures-and-cubatures-with-dirac-deltas/108683/5 "2024-01-11T17:18:31Z")

</div>

> [@homocomputeris](#):
>
> Can I make cubatures step into these points somehow?

I don’t know of any integration packages that handle delta functions _automatically_, but as @tom-plaa wrote it should be easy to include the delta contributions separately from the rest of your integral.

(A separate issue arises when you have _lots_ of delta functions, or even an infinite summation. In that case it’s possible to develop a quadrature scheme for a [discrete measure](https://en.wikipedia.org/wiki/Discrete_measure) that evaluates your integrand at _fewer_ points than the number of delta functions, essentially by deriving a specialized interpolation scheme. For example, a student of mine did this to derive specialized [“quadrature” schemes for infinite Matsubara sums](https://arxiv.org/abs/2101.01655) in statistical physics.)

---

<div class="post-metadata">

**Author:** ![svretina](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/svretina/32/206157_2.png) [@svretina](https://discourse.julialang.org/u/svretina)\
**Post date:** [January 11, 2024, 5:43pm UTC](https://discourse.julialang.org/t/quadratures-and-cubatures-with-dirac-deltas/108683/6 "2024-01-11T17:43:33Z")

</div>

If you are interested in quadrature rules that can handle singularities,  
then you can use the [Tanh-Sinh](https://en.wikipedia.org/wiki/Tanh-sinh_quadrature) rule. Or do you want to use the dirac function properties in integrals?

---

<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:** [January 11, 2024, 5:59pm UTC](https://discourse.julialang.org/t/quadratures-and-cubatures-with-dirac-deltas/108683/7 "2024-01-11T17:59:51Z")

</div>

> [@svretina](#):
>
> If you are interested in quadrature rules that can handle singularities,  
> then you can use the [Tanh-Sinh](https://en.wikipedia.org/wiki/Tanh-sinh_quadrature) rule.

That’s only for endpoint singularities.

---

<div class="post-metadata">

**Author:** ![homocomputeris](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/homocomputeris/32/8933_2.png) [@homocomputeris](https://discourse.julialang.org/u/homocomputeris)\
**Post date:** [January 11, 2024, 6:08pm UTC](https://discourse.julialang.org/t/quadratures-and-cubatures-with-dirac-deltas/108683/8 "2024-01-11T18:08:14Z")

</div>

> [@dlfivefifty](#):
>
> How are your functions represented?

> [@svretina](#):
>
> Or do you want to use the dirac function properties in integrals?

In general, I want to handle functions like where there can be several deltas in different variables, e.g. for a 3D case:

f(x,y,z) = h(x)g(y)\delta((a(z))

---

<div class="post-metadata">

**Author:** ![svretina](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/svretina/32/206157_2.png) [@svretina](https://discourse.julialang.org/u/svretina)\
**Post date:** [January 11, 2024, 6:13pm UTC](https://discourse.julialang.org/t/quadratures-and-cubatures-with-dirac-deltas/108683/9 "2024-01-11T18:13:48Z")

</div>

yes indeed, but then you just have to split your domain in subdomains and then sum the respective results from the subdomains.

---

<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:** [January 11, 2024, 6:37pm UTC](https://discourse.julialang.org/t/quadratures-and-cubatures-with-dirac-deltas/108683/10 "2024-01-11T18:37:17Z")

</div>

> [@svretina](#):
>
> yes indeed, but then you just have to split your domain in subdomains and then sum the respective results from the subdomains.

Still has nothing to do with delta functions, which can’t be represented as ordinary functions at all (because they aren’t — they are distributions). Tanh–sinh methods are really offtopic for the current thread.

---

<div class="post-metadata">

**Author:** ![svretina](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/svretina/32/206157_2.png) [@svretina](https://discourse.julialang.org/u/svretina)\
**Post date:** [January 11, 2024, 6:40pm UTC](https://discourse.julialang.org/t/quadratures-and-cubatures-with-dirac-deltas/108683/11 "2024-01-11T18:40:44Z")

</div>

I agree, that’s why I asked if the OP needs something with regard to singularities or he/she wants something about Dirac functions solely

---

<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:** [January 11, 2024, 6:46pm UTC](https://discourse.julialang.org/t/quadratures-and-cubatures-with-dirac-deltas/108683/12 "2024-01-11T18:46:46Z")

</div>

> [@homocomputeris](#):
>
> In general, I want to handle functions like where there can be several deltas in different variables, e.g. for a 3D case: f(x,y,z) = h(x)g(y)\delta((a(z))

Oh, that can get much more tricky, because that sort of thing ultimately boils down to quadrature over [implicit surfaces](https://en.wikipedia.org/wiki/Implicit_surface). I don’t know of any Julia package for this offhand, though there are algorithms in the literature like [Saye (2015)](https://epubs.siam.org/doi/10.1137/140966290).

For that particular example (assuming integration over a simple xyz box) it simplifies because it is separable, giving:

\left(\sum\_k \frac{1}{|a'(z\_k)|}\right) \iint h(x) g(y) dx dy

where |a'| is a Jacobian factor and z\_k is the k-th root of a(z) in your domain (which you can find numerically by one of the many root-finding methods in Julia).

But in more general cases it won’t be separable and you will be back to implicit surfaces. Depending on how complicated your problem is, you’ll have to decide if you can simplify it analytically like this or if you need some more general method.

---

<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:** [January 11, 2024, 8:53pm UTC](https://discourse.julialang.org/t/quadratures-and-cubatures-with-dirac-deltas/108683/13 "2024-01-11T20:53:35Z")

</div>

A post was split to a new topic: [Guide to numerical integration methods?](https://discourse.julialang.org/t/guide-to-numerical-integration-methods/108704)