Are there ways to compute cubatures of functions that include the Dirac delta function?
How are your functions represented? How would you be able to detect a delta function if you can only evaluate pointwise?
1 Like
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.
2 Likes
Yeah, my bad. I know where the deltas are. 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 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 in statistical physics.)
1 Like
If you are interested in quadrature rules that can handle singularities,
then you can use the Tanh-Sinh rule. Or do you want to use the dirac function properties in integrals?
That’s only for endpoint singularities.
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))
yes indeed, but then you just have to split your domain in subdomains and then sum the respective results from the subdomains.
1 Like
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.
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
Oh, that can get much more tricky, because that sort of thing ultimately boils down to quadrature over implicit surfaces. I don’t know of any Julia package for this offhand, though there are algorithms in the literature like Saye (2015).
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.
2 Likes