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[1]:= F0 = Integrate[t^n (x - t)^m, t]
Out[1]= (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[2]:= 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[3]:= jmp = Limit[F0, t -> x] - Limit[F0, t -> x, Direction -> 1]
Out[3]= ((-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[4]:= int = (F /. t -> max1) - (F /. t -> min1) -
jmp*Boole[min2 < 0 < max2 && min1 < x < max1];
Here's a check:
In[5]:= B[x_, {min_, max_}] := Boole[min < x < max]
In[6]:= 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[6]= 0.000021
Maxim Rytin
m.r at inbox.ru