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

Factor 2 error in Inverse Laplace Transform

• To: mathgroup at smc.vnet.net
• Subject: [mg52797] Factor 2 error in Inverse Laplace Transform
• From: ab_def at prontomail.com (Maxim)
• Date: Mon, 13 Dec 2004 04:23:11 -0500 (EST)
• Sender: owner-wri-mathgroup at wolfram.com

There was a discussion of a curious LaplaceTransform glitch a while
ago in this group ( http://forums.wolfram.com/mathgroup/archive/2004/Oct/msg00248.html
). Here is a similar example in version 5.1:

In[1]:=
InverseLaplaceTransform[p*(PolyGamma[p + 2] - PolyGamma[p]), p, t]

Out[1]=
2 - 1/(2*E^t)

Since PolyGamma[p + 2] - PolyGamma[p] == 1/p + 1/(p + 1), the correct
answer is 2*DiracDelta[t] - E^(-t). It is interesting to trace where
the error comes from. Let

f(t) = 1 + E^(-t), F(p) = L[f] = 1/p + 1/(p + 1).

Mathematica evaluates the inverse transform of p*F(p) by using the
identity

(L^-1)[p*F(p)] = convolution((L^-1)[p], (L^-1)[F(p)]) =

convolution(DiracDelta'(t), f(t)).

Then it writes the convolution as

In[2]:=
Integrate[(1 + E^(-tau))*DiracDelta'[t - tau], {tau, low, t},
  Assumptions -> low < 0]

Out[2]=
-1/(2*E^t) + DiracDelta[-low + t] + DiracDelta[-low + t]/E^low

where the parameter low is used to 'include' the zero in the
integration range, after which Mathematica takes the limit as low
tends to 0 from the left (all this can be seen by doing Trace).

The first error is that in fact it is needed to include the point
tau=t, not tau=0; otherwise Mathematica multiplies the value of the
function at tau=t by 1/2: in Mathematica Integrate[DiracDelta[t], {t,
0, Infinity}] == 1/2, and similarly in the above case we get -E^(-t)/2
instead of -E^(-t).

The next error is encountered when taking the limit:

In[3]:=
Limit[-1/(2*E^t) + DiracDelta[-low + t] + DiracDelta[-low + t]/E^low,
  low -> 0, Direction -> 1]

Out[3]=
2 - 1/(2*E^t)

which is the answer returned by InverseLaplaceTransform; the
DiracDelta terms somehow evaluated to 2. At the same time, taking the
limit term by term gives

In[4]:=
Limit[#, low -> 0, Direction -> 1]& /@
  (-1/(2*E^t) + DiracDelta[-low + t] + DiracDelta[-low + t]/E^low)

Out[4]=
-1/(2*E^t) + 2*DiracDelta[t, t]

which is somewhat closer to the correct answer, but also incorrect:
DiracDelta[t, t] is the same as DiracDelta[t]^2, which is
indeterminate -- it contains a product of distributions with
intersecting supports.

In fact, rather than performing all those dubious transformations such
as taking pointwise limits of generalized functions, we could evaluate
the convolution in a straightforward way: it is simply

In[5]:=
Integrate[(1 + E^(-tau))*UnitStep[tau]*DiracDelta'[t - tau],
  {tau, -Infinity, Infinity}]

Out[5]=
-(t + Sqrt[t^2])/(2*E^t*Sqrt[t^2])

When evaluated correctly, with UnitStep viewed as a generalized
function, this should give the required answer 2*DiracDelta[t] -
E^(-t). Unfortunately, Mathematica is unable to handle this integral
as well (although, curiously, there is a way to make it work -- by
rewriting the integral as Integrate[(1 + E^(-tau))*DiracDelta'[t -
tau], {tau, 0, Infinity}]).

The second example from Peter's post still gives an incorrect result:

In[6]:=
InverseLaplaceTransform[p*ArcTan[1/p], p, t]

Out[6]=
(2*t^2 + t*Cos[t] - Sin[t])/(2*t^2)

The correct value is DiracDelta[t] + (t*Cos[t] - Sin[t])/t^2. Besides,
if we try to verify Out[6] by taking the direct transform, we find out
that LaplaceTransform of (t*Cos[t] - Sin[t])/t^2 is incorrect too:

In[7]:=
LaplaceTransform[(t*Cos[t]-Sin[t])/t^2, t, p] // FullSimplify

Out[7]=
p*ArcCot[p]

The correct value is p*ArcCot[p] - 1.

The only conclusion that comes to mind is something along the lines
"if it were a horse..."

Maxim Rytin
m.r at inbox.ru