[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
Prev by Date:
**Re: Bug in 5.1,Linux, GUIKit?**
Next by Date:
**Re: Solve Feature?**
Previous by thread:
**Graphics Conversion Bug with DensityGraphics**
Next by thread:
**Re: Solve Feature?**
| |