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