MathGroup Archive 2007

[Date Index] [Thread Index] [Author Index]

Search the Archive

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



  • 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