       Re: Evaluating a convolution integral in Mathematica

• To: mathgroup at smc.vnet.net
• Subject: [mg80284] Re: Evaluating a convolution integral in Mathematica
• From: m.r at inbox.ru
• Date: Thu, 16 Aug 2007 07:20:14 -0400 (EDT)
• References: <f9ud4q\$a5q\$1@smc.vnet.net>

```On Aug 15, 3:22 am, rdb... at indiana.edu wrote:
> I am having trouble evaluating a convolution integral in Mathematica.
> Define B[x,{min,max}] to be a function that takes on the value 1 when
> x is between min and max and 0 otherwise.  I need to find the
> convolution of x^n B[x, {min1, max1}] with x^m B[x, {min2, max2}],
> where x is Real, n and m are nonnegative integers, and the mins and
> maxs are Real with the constraints that min1<=max1 and min2<=max2.
>
> The resulting convolution integeral is Integrate[t^n B[t, {min1,
> max1}] (x-t)^m B[x-t, {min2, max2}], {t, -Infinity, Infinity}].
>
> Mathematica 6.0.1 has no problem evaluating this integral when
> constant values are substituted for n, m, and the mins and maxs.
> However, I need the general value of this integral.  Mathematica also
> claims to be able to evaluate this general integral, returning a
> complicated peicewise expression involving gamma functions and
> hypergeometric functions.  However, when specific values for n, m and
> the mins and maxs are then substituted into this general expression,
> it always returns either 0 or Indeterminate.
>
> Any help with evaluating this integral in general would be greatly
> appreciated.

We need to evaluate

Integrate[t^n (x - t)^m Boole[min2 < x - t < max2], {t, min1, max1}].

Find the antiderivative:

In:= F0 = Integrate[t^n (x - t)^m, t]

Out= (1/(1 + n))t^(1 + n) (1 - t/x)^-m (-t + x)^m
Hypergeometric2F1[1 + n, -m, 2 + n, t/x]

Make it constant outside the interval min2 < x - t < max2 and
continuous at x - t == min2 and x - t == max2:

In:= F = Piecewise[{{F0 /. t -> x - max2, t < x - max2},
{F0 /. t -> x - min2, x - min2 < t},
{F0, True}}];

F can also be discontinuous at t == x:

In:= jmp = Limit[F0, t -> x] - Limit[F0, t -> x, Direction -> 1]

Out= ((-1)^m (-(1/x))^-m x^(1 + n) Gamma[1 + m] Gamma[2 + n])/((1 +
n) Gamma[2 + m + n]) - ((1/x)^-m x^(1 + n) Gamma[1 + m] Gamma[2 + n])/
((1 + n) Gamma[2 + m + n])

The jump should be taken into account only when the point t == x is
inside the interval min2 < x - t < max2 and also inside the interval
of integration min1 < t < max1:

In:= int = (F /. t -> max1) - (F /. t -> min1) -
jmp*Boole[min2 < 0 < max2 && min1 < x < max1];

Here's a check:

In:= B[x_, {min_, max_}] := Boole[min < x < max]

In:= Block[{m, n, x, min1, max1, min2, max2,
prec = WorkingPrecision -> 20},
Table[
{m, n} = RandomReal[{-1, 10}, 2, prec];
x = RandomReal[{-10, 10}, prec];
{min1, max1} = Sort@ RandomReal[{-10, 10}, 2, prec];
{min2, max2} = Sort@ RandomReal[{-10, 10}, 2, prec];
int - NIntegrate[
t^n B[t, {min1, max1}] (x - t)^m B[x - t, {min2, max2}],
{t, -Infinity,
Sequence @@ Sort@ {0, x, min1, max1, x - min2, x - max2},
Infinity} // Evaluate,
AccuracyGoal -> 5, MaxRecursion -> 15, prec // Evaluate],
{100}]] // Abs // Max

Out= 0.000021

Maxim Rytin
m.r at inbox.ru

```

• Prev by Date: Re: Extracting terms of a multivariate polynomial order by order
• Next by Date: Re: Extracting terms of a multivariate polynomial order by order
• Previous by thread: Re: Evaluating a convolution integral in Mathematica
• Next by thread: Re: Question about memory consumption in my codes by using Mathematica