Coverage for /usr/lib/python3/dist-packages/sympy/integrals/laplace.py: 9%
881 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1"""Laplace Transforms"""
2from sympy.core import S, pi, I
3from sympy.core.add import Add
4from sympy.core.cache import cacheit
5from sympy.core.function import (
6 AppliedUndef, Derivative, expand, expand_complex, expand_mul, expand_trig,
7 Lambda, WildFunction, diff)
8from sympy.core.mul import Mul, prod
9from sympy.core.relational import _canonical, Ge, Gt, Lt, Unequality, Eq
10from sympy.core.sorting import ordered
11from sympy.core.symbol import Dummy, symbols, Wild
12from sympy.functions.elementary.complexes import (
13 re, im, arg, Abs, polar_lift, periodic_argument)
14from sympy.functions.elementary.exponential import exp, log
15from sympy.functions.elementary.hyperbolic import cosh, coth, sinh, asinh
16from sympy.functions.elementary.miscellaneous import Max, Min, sqrt
17from sympy.functions.elementary.piecewise import Piecewise
18from sympy.functions.elementary.trigonometric import cos, sin, atan
19from sympy.functions.special.bessel import besseli, besselj, besselk, bessely
20from sympy.functions.special.delta_functions import DiracDelta, Heaviside
21from sympy.functions.special.error_functions import erf, erfc, Ei
22from sympy.functions.special.gamma_functions import digamma, gamma, lowergamma
23from sympy.integrals import integrate, Integral
24from sympy.integrals.transforms import (
25 _simplify, IntegralTransform, IntegralTransformError)
26from sympy.logic.boolalg import to_cnf, conjuncts, disjuncts, Or, And
27from sympy.matrices.matrices import MatrixBase
28from sympy.polys.matrices.linsolve import _lin_eq2dict
29from sympy.polys.polyerrors import PolynomialError
30from sympy.polys.polyroots import roots
31from sympy.polys.polytools import Poly
32from sympy.polys.rationaltools import together
33from sympy.polys.rootoftools import RootSum
34from sympy.utilities.exceptions import (
35 sympy_deprecation_warning, SymPyDeprecationWarning, ignore_warnings)
36from sympy.utilities.misc import debug, debugf
39def _simplifyconds(expr, s, a):
40 r"""
41 Naively simplify some conditions occurring in ``expr``,
42 given that `\operatorname{Re}(s) > a`.
44 Examples
45 ========
47 >>> from sympy.integrals.laplace import _simplifyconds
48 >>> from sympy.abc import x
49 >>> from sympy import sympify as S
50 >>> _simplifyconds(abs(x**2) < 1, x, 1)
51 False
52 >>> _simplifyconds(abs(x**2) < 1, x, 2)
53 False
54 >>> _simplifyconds(abs(x**2) < 1, x, 0)
55 Abs(x**2) < 1
56 >>> _simplifyconds(abs(1/x**2) < 1, x, 1)
57 True
58 >>> _simplifyconds(S(1) < abs(x), x, 1)
59 True
60 >>> _simplifyconds(S(1) < abs(1/x), x, 1)
61 False
63 >>> from sympy import Ne
64 >>> _simplifyconds(Ne(1, x**3), x, 1)
65 True
66 >>> _simplifyconds(Ne(1, x**3), x, 2)
67 True
68 >>> _simplifyconds(Ne(1, x**3), x, 0)
69 Ne(1, x**3)
70 """
72 def power(ex):
73 if ex == s:
74 return 1
75 if ex.is_Pow and ex.base == s:
76 return ex.exp
77 return None
79 def bigger(ex1, ex2):
80 """ Return True only if |ex1| > |ex2|, False only if |ex1| < |ex2|.
81 Else return None. """
82 if ex1.has(s) and ex2.has(s):
83 return None
84 if isinstance(ex1, Abs):
85 ex1 = ex1.args[0]
86 if isinstance(ex2, Abs):
87 ex2 = ex2.args[0]
88 if ex1.has(s):
89 return bigger(1/ex2, 1/ex1)
90 n = power(ex2)
91 if n is None:
92 return None
93 try:
94 if n > 0 and (Abs(ex1) <= Abs(a)**n) == S.true:
95 return False
96 if n < 0 and (Abs(ex1) >= Abs(a)**n) == S.true:
97 return True
98 except TypeError:
99 pass
101 def replie(x, y):
102 """ simplify x < y """
103 if (not (x.is_positive or isinstance(x, Abs))
104 or not (y.is_positive or isinstance(y, Abs))):
105 return (x < y)
106 r = bigger(x, y)
107 if r is not None:
108 return not r
109 return (x < y)
111 def replue(x, y):
112 b = bigger(x, y)
113 if b in (True, False):
114 return True
115 return Unequality(x, y)
117 def repl(ex, *args):
118 if ex in (True, False):
119 return bool(ex)
120 return ex.replace(*args)
121 from sympy.simplify.radsimp import collect_abs
122 expr = collect_abs(expr)
123 expr = repl(expr, Lt, replie)
124 expr = repl(expr, Gt, lambda x, y: replie(y, x))
125 expr = repl(expr, Unequality, replue)
126 return S(expr)
129def expand_dirac_delta(expr):
130 """
131 Expand an expression involving DiractDelta to get it as a linear
132 combination of DiracDelta functions.
133 """
134 return _lin_eq2dict(expr, expr.atoms(DiracDelta))
137def _laplace_transform_integration(f, t, s_, simplify=True):
138 """ The backend function for doing Laplace transforms by integration.
140 This backend assumes that the frontend has already split sums
141 such that `f` is to an addition anymore.
142 """
143 s = Dummy('s')
144 debugf('[LT _l_t_i ] started with (%s, %s, %s)', (f, t, s))
145 debugf('[LT _l_t_i ] and simplify=%s', (simplify, ))
147 if f.has(DiracDelta):
148 return None
150 F = integrate(f*exp(-s*t), (t, S.Zero, S.Infinity))
151 debugf('[LT _l_t_i ] integrated: %s', (F, ))
153 if not F.has(Integral):
154 return _simplify(F.subs(s, s_), simplify), S.NegativeInfinity, S.true
156 if not F.is_Piecewise:
157 debug('[LT _l_t_i ] not piecewise.')
158 return None
160 F, cond = F.args[0]
161 if F.has(Integral):
162 debug('[LT _l_t_i ] integral in unexpected form.')
163 return None
165 def process_conds(conds):
166 """ Turn ``conds`` into a strip and auxiliary conditions. """
167 from sympy.solvers.inequalities import _solve_inequality
168 a = S.NegativeInfinity
169 aux = S.true
170 conds = conjuncts(to_cnf(conds))
171 p, q, w1, w2, w3, w4, w5 = symbols(
172 'p q w1 w2 w3 w4 w5', cls=Wild, exclude=[s])
173 patterns = (
174 p*Abs(arg((s + w3)*q)) < w2,
175 p*Abs(arg((s + w3)*q)) <= w2,
176 Abs(periodic_argument((s + w3)**p*q, w1)) < w2,
177 Abs(periodic_argument((s + w3)**p*q, w1)) <= w2,
178 Abs(periodic_argument((polar_lift(s + w3))**p*q, w1)) < w2,
179 Abs(periodic_argument((polar_lift(s + w3))**p*q, w1)) <= w2)
180 for c in conds:
181 a_ = S.Infinity
182 aux_ = []
183 for d in disjuncts(c):
184 if d.is_Relational and s in d.rhs.free_symbols:
185 d = d.reversed
186 if d.is_Relational and isinstance(d, (Ge, Gt)):
187 d = d.reversedsign
188 for pat in patterns:
189 m = d.match(pat)
190 if m:
191 break
192 if m and m[q].is_positive and m[w2]/m[p] == pi/2:
193 d = -re(s + m[w3]) < 0
194 m = d.match(p - cos(w1*Abs(arg(s*w5))*w2)*Abs(s**w3)**w4 < 0)
195 if not m:
196 m = d.match(
197 cos(p - Abs(periodic_argument(s**w1*w5, q))*w2) *
198 Abs(s**w3)**w4 < 0)
199 if not m:
200 m = d.match(
201 p - cos(
202 Abs(periodic_argument(polar_lift(s)**w1*w5, q))*w2
203 )*Abs(s**w3)**w4 < 0)
204 if m and all(m[wild].is_positive for wild in [
205 w1, w2, w3, w4, w5]):
206 d = re(s) > m[p]
207 d_ = d.replace(
208 re, lambda x: x.expand().as_real_imag()[0]).subs(re(s), t)
209 if (
210 not d.is_Relational or d.rel_op in ('==', '!=')
211 or d_.has(s) or not d_.has(t)):
212 aux_ += [d]
213 continue
214 soln = _solve_inequality(d_, t)
215 if not soln.is_Relational or soln.rel_op in ('==', '!='):
216 aux_ += [d]
217 continue
218 if soln.lts == t:
219 debug('[LT _l_t_i ] convergence not in half-plane.')
220 return None
221 else:
222 a_ = Min(soln.lts, a_)
223 if a_ is not S.Infinity:
224 a = Max(a_, a)
225 else:
226 aux = And(aux, Or(*aux_))
227 return a, aux.canonical if aux.is_Relational else aux
229 conds = [process_conds(c) for c in disjuncts(cond)]
230 conds2 = [x for x in conds if x[1] !=
231 S.false and x[0] is not S.NegativeInfinity]
232 if not conds2:
233 conds2 = [x for x in conds if x[1] != S.false]
234 conds = list(ordered(conds2))
236 def cnt(expr):
237 if expr in (True, False):
238 return 0
239 return expr.count_ops()
240 conds.sort(key=lambda x: (-x[0], cnt(x[1])))
242 if not conds:
243 debug('[LT _l_t_i ] no convergence found.')
244 return None
245 a, aux = conds[0] # XXX is [0] always the right one?
247 def sbs(expr):
248 return expr.subs(s, s_)
249 if simplify:
250 F = _simplifyconds(F, s, a)
251 aux = _simplifyconds(aux, s, a)
252 return _simplify(F.subs(s, s_), simplify), sbs(a), _canonical(sbs(aux))
255def _laplace_deep_collect(f, t):
256 """
257 This is an internal helper function that traverses through the epression
258 tree of `f(t)` and collects arguments. The purpose of it is that
259 anything like `f(w*t-1*t-c)` will be written as `f((w-1)*t-c)` such that
260 it can match `f(a*t+b)`.
261 """
262 func = f.func
263 args = list(f.args)
264 if len(f.args) == 0:
265 return f
266 else:
267 args = [_laplace_deep_collect(arg, t) for arg in args]
268 if func.is_Add:
269 return func(*args).collect(t)
270 else:
271 return func(*args)
274@cacheit
275def _laplace_build_rules():
276 """
277 This is an internal helper function that returns the table of Laplace
278 transform rules in terms of the time variable `t` and the frequency
279 variable `s`. It is used by ``_laplace_apply_rules``. Each entry is a
280 tuple containing:
282 (time domain pattern,
283 frequency-domain replacement,
284 condition for the rule to be applied,
285 convergence plane,
286 preparation function)
288 The preparation function is a function with one argument that is applied
289 to the expression before matching. For most rules it should be
290 ``_laplace_deep_collect``.
291 """
292 t = Dummy('t')
293 s = Dummy('s')
294 a = Wild('a', exclude=[t])
295 b = Wild('b', exclude=[t])
296 n = Wild('n', exclude=[t])
297 tau = Wild('tau', exclude=[t])
298 omega = Wild('omega', exclude=[t])
299 def dco(f): return _laplace_deep_collect(f, t)
300 debug('_laplace_build_rules is building rules')
302 laplace_transform_rules = [
303 (a, a/s,
304 S.true, S.Zero, dco), # 4.2.1
305 (DiracDelta(a*t-b), exp(-s*b/a)/Abs(a),
306 Or(And(a > 0, b >= 0), And(a < 0, b <= 0)),
307 S.NegativeInfinity, dco), # Not in Bateman54
308 (DiracDelta(a*t-b), S(0),
309 Or(And(a < 0, b >= 0), And(a > 0, b <= 0)),
310 S.NegativeInfinity, dco), # Not in Bateman54
311 (Heaviside(a*t-b), exp(-s*b/a)/s,
312 And(a > 0, b > 0), S.Zero, dco), # 4.4.1
313 (Heaviside(a*t-b), (1-exp(-s*b/a))/s,
314 And(a < 0, b < 0), S.Zero, dco), # 4.4.1
315 (Heaviside(a*t-b), 1/s,
316 And(a > 0, b <= 0), S.Zero, dco), # 4.4.1
317 (Heaviside(a*t-b), 0,
318 And(a < 0, b > 0), S.Zero, dco), # 4.4.1
319 (t, 1/s**2,
320 S.true, S.Zero, dco), # 4.2.3
321 (1/(a*t+b), -exp(-b/a*s)*Ei(-b/a*s)/a,
322 Abs(arg(b/a)) < pi, S.Zero, dco), # 4.2.6
323 (1/sqrt(a*t+b), sqrt(a*pi/s)*exp(b/a*s)*erfc(sqrt(b/a*s))/a,
324 Abs(arg(b/a)) < pi, S.Zero, dco), # 4.2.18
325 ((a*t+b)**(-S(3)/2),
326 2*b**(-S(1)/2)-2*(pi*s/a)**(S(1)/2)*exp(b/a*s) * erfc(sqrt(b/a*s))/a,
327 Abs(arg(b/a)) < pi, S.Zero, dco), # 4.2.20
328 (sqrt(t)/(t+b), sqrt(pi/s)-pi*sqrt(b)*exp(b*s)*erfc(sqrt(b*s)),
329 Abs(arg(b)) < pi, S.Zero, dco), # 4.2.22
330 (1/(a*sqrt(t) + t**(3/2)), pi*a**(S(1)/2)*exp(a*s)*erfc(sqrt(a*s)),
331 S.true, S.Zero, dco), # Not in Bateman54
332 (t**n, gamma(n+1)/s**(n+1),
333 n > -1, S.Zero, dco), # 4.3.1
334 ((a*t+b)**n, lowergamma(n+1, b/a*s)*exp(-b/a*s)/s**(n+1)/a,
335 And(n > -1, Abs(arg(b/a)) < pi), S.Zero, dco), # 4.3.4
336 (t**n/(t+a), a**n*gamma(n+1)*lowergamma(-n, a*s),
337 And(n > -1, Abs(arg(a)) < pi), S.Zero, dco), # 4.3.7
338 (exp(a*t-tau), exp(-tau)/(s-a),
339 S.true, re(a), dco), # 4.5.1
340 (t*exp(a*t-tau), exp(-tau)/(s-a)**2,
341 S.true, re(a), dco), # 4.5.2
342 (t**n*exp(a*t), gamma(n+1)/(s-a)**(n+1),
343 re(n) > -1, re(a), dco), # 4.5.3
344 (exp(-a*t**2), sqrt(pi/4/a)*exp(s**2/4/a)*erfc(s/sqrt(4*a)),
345 re(a) > 0, S.Zero, dco), # 4.5.21
346 (t*exp(-a*t**2),
347 1/(2*a)-2/sqrt(pi)/(4*a)**(S(3)/2)*s*erfc(s/sqrt(4*a)),
348 re(a) > 0, S.Zero, dco), # 4.5.22
349 (exp(-a/t), 2*sqrt(a/s)*besselk(1, 2*sqrt(a*s)),
350 re(a) >= 0, S.Zero, dco), # 4.5.25
351 (sqrt(t)*exp(-a/t),
352 S(1)/2*sqrt(pi/s**3)*(1+2*sqrt(a*s))*exp(-2*sqrt(a*s)),
353 re(a) >= 0, S.Zero, dco), # 4.5.26
354 (exp(-a/t)/sqrt(t), sqrt(pi/s)*exp(-2*sqrt(a*s)),
355 re(a) >= 0, S.Zero, dco), # 4.5.27
356 (exp(-a/t)/(t*sqrt(t)), sqrt(pi/a)*exp(-2*sqrt(a*s)),
357 re(a) > 0, S.Zero, dco), # 4.5.28
358 (t**n*exp(-a/t), 2*(a/s)**((n+1)/2)*besselk(n+1, 2*sqrt(a*s)),
359 re(a) > 0, S.Zero, dco), # 4.5.29
360 (exp(-2*sqrt(a*t)),
361 s**(-1)-sqrt(pi*a)*s**(-S(3)/2)*exp(a/s) * erfc(sqrt(a/s)),
362 Abs(arg(a)) < pi, S.Zero, dco), # 4.5.31
363 (exp(-2*sqrt(a*t))/sqrt(t), (pi/s)**(S(1)/2)*exp(a/s)*erfc(sqrt(a/s)),
364 Abs(arg(a)) < pi, S.Zero, dco), # 4.5.33
365 (log(a*t), -log(exp(S.EulerGamma)*s/a)/s,
366 a > 0, S.Zero, dco), # 4.6.1
367 (log(1+a*t), -exp(s/a)/s*Ei(-s/a),
368 Abs(arg(a)) < pi, S.Zero, dco), # 4.6.4
369 (log(a*t+b), (log(b)-exp(s/b/a)/s*a*Ei(-s/b))/s*a,
370 And(a > 0, Abs(arg(b)) < pi), S.Zero, dco), # 4.6.5
371 (log(t)/sqrt(t), -sqrt(pi/s)*log(4*s*exp(S.EulerGamma)),
372 S.true, S.Zero, dco), # 4.6.9
373 (t**n*log(t), gamma(n+1)*s**(-n-1)*(digamma(n+1)-log(s)),
374 re(n) > -1, S.Zero, dco), # 4.6.11
375 (log(a*t)**2, (log(exp(S.EulerGamma)*s/a)**2+pi**2/6)/s,
376 a > 0, S.Zero, dco), # 4.6.13
377 (sin(omega*t), omega/(s**2+omega**2),
378 S.true, Abs(im(omega)), dco), # 4,7,1
379 (Abs(sin(omega*t)), omega/(s**2+omega**2)*coth(pi*s/2/omega),
380 omega > 0, S.Zero, dco), # 4.7.2
381 (sin(omega*t)/t, atan(omega/s),
382 S.true, Abs(im(omega)), dco), # 4.7.16
383 (sin(omega*t)**2/t, log(1+4*omega**2/s**2)/4,
384 S.true, 2*Abs(im(omega)), dco), # 4.7.17
385 (sin(omega*t)**2/t**2,
386 omega*atan(2*omega/s)-s*log(1+4*omega**2/s**2)/4,
387 S.true, 2*Abs(im(omega)), dco), # 4.7.20
388 (sin(2*sqrt(a*t)), sqrt(pi*a)/s/sqrt(s)*exp(-a/s),
389 S.true, S.Zero, dco), # 4.7.32
390 (sin(2*sqrt(a*t))/t, pi*erf(sqrt(a/s)),
391 S.true, S.Zero, dco), # 4.7.34
392 (cos(omega*t), s/(s**2+omega**2),
393 S.true, Abs(im(omega)), dco), # 4.7.43
394 (cos(omega*t)**2, (s**2+2*omega**2)/(s**2+4*omega**2)/s,
395 S.true, 2*Abs(im(omega)), dco), # 4.7.45
396 (sqrt(t)*cos(2*sqrt(a*t)), sqrt(pi)/2*s**(-S(5)/2)*(s-2*a)*exp(-a/s),
397 S.true, S.Zero, dco), # 4.7.66
398 (cos(2*sqrt(a*t))/sqrt(t), sqrt(pi/s)*exp(-a/s),
399 S.true, S.Zero, dco), # 4.7.67
400 (sin(a*t)*sin(b*t), 2*a*b*s/(s**2+(a+b)**2)/(s**2+(a-b)**2),
401 S.true, Abs(im(a))+Abs(im(b)), dco), # 4.7.78
402 (cos(a*t)*sin(b*t), b*(s**2-a**2+b**2)/(s**2+(a+b)**2)/(s**2+(a-b)**2),
403 S.true, Abs(im(a))+Abs(im(b)), dco), # 4.7.79
404 (cos(a*t)*cos(b*t), s*(s**2+a**2+b**2)/(s**2+(a+b)**2)/(s**2+(a-b)**2),
405 S.true, Abs(im(a))+Abs(im(b)), dco), # 4.7.80
406 (sinh(a*t), a/(s**2-a**2),
407 S.true, Abs(re(a)), dco), # 4.9.1
408 (cosh(a*t), s/(s**2-a**2),
409 S.true, Abs(re(a)), dco), # 4.9.2
410 (sinh(a*t)**2, 2*a**2/(s**3-4*a**2*s),
411 S.true, 2*Abs(re(a)), dco), # 4.9.3
412 (cosh(a*t)**2, (s**2-2*a**2)/(s**3-4*a**2*s),
413 S.true, 2*Abs(re(a)), dco), # 4.9.4
414 (sinh(a*t)/t, log((s+a)/(s-a))/2,
415 S.true, Abs(re(a)), dco), # 4.9.12
416 (t**n*sinh(a*t), gamma(n+1)/2*((s-a)**(-n-1)-(s+a)**(-n-1)),
417 n > -2, Abs(a), dco), # 4.9.18
418 (t**n*cosh(a*t), gamma(n+1)/2*((s-a)**(-n-1)+(s+a)**(-n-1)),
419 n > -1, Abs(a), dco), # 4.9.19
420 (sinh(2*sqrt(a*t)), sqrt(pi*a)/s/sqrt(s)*exp(a/s),
421 S.true, S.Zero, dco), # 4.9.34
422 (cosh(2*sqrt(a*t)), 1/s+sqrt(pi*a)/s/sqrt(s)*exp(a/s)*erf(sqrt(a/s)),
423 S.true, S.Zero, dco), # 4.9.35
424 (
425 sqrt(t)*sinh(2*sqrt(a*t)),
426 pi**(S(1)/2)*s**(-S(5)/2)*(s/2+a) *
427 exp(a/s)*erf(sqrt(a/s))-a**(S(1)/2)*s**(-2),
428 S.true, S.Zero, dco), # 4.9.36
429 (sqrt(t)*cosh(2*sqrt(a*t)), pi**(S(1)/2)*s**(-S(5)/2)*(s/2+a)*exp(a/s),
430 S.true, S.Zero, dco), # 4.9.37
431 (sinh(2*sqrt(a*t))/sqrt(t),
432 pi**(S(1)/2)*s**(-S(1)/2)*exp(a/s) * erf(sqrt(a/s)),
433 S.true, S.Zero, dco), # 4.9.38
434 (cosh(2*sqrt(a*t))/sqrt(t), pi**(S(1)/2)*s**(-S(1)/2)*exp(a/s),
435 S.true, S.Zero, dco), # 4.9.39
436 (sinh(sqrt(a*t))**2/sqrt(t), pi**(S(1)/2)/2*s**(-S(1)/2)*(exp(a/s)-1),
437 S.true, S.Zero, dco), # 4.9.40
438 (cosh(sqrt(a*t))**2/sqrt(t), pi**(S(1)/2)/2*s**(-S(1)/2)*(exp(a/s)+1),
439 S.true, S.Zero, dco), # 4.9.41
440 (erf(a*t), exp(s**2/(2*a)**2)*erfc(s/(2*a))/s,
441 4*Abs(arg(a)) < pi, S.Zero, dco), # 4.12.2
442 (erf(sqrt(a*t)), sqrt(a)/sqrt(s+a)/s,
443 S.true, Max(S.Zero, -re(a)), dco), # 4.12.4
444 (exp(a*t)*erf(sqrt(a*t)), sqrt(a)/sqrt(s)/(s-a),
445 S.true, Max(S.Zero, re(a)), dco), # 4.12.5
446 (erf(sqrt(a/t)/2), (1-exp(-sqrt(a*s)))/s,
447 re(a) > 0, S.Zero, dco), # 4.12.6
448 (erfc(sqrt(a*t)), (sqrt(s+a)-sqrt(a))/sqrt(s+a)/s,
449 S.true, -re(a), dco), # 4.12.9
450 (exp(a*t)*erfc(sqrt(a*t)), 1/(s+sqrt(a*s)),
451 S.true, S.Zero, dco), # 4.12.10
452 (erfc(sqrt(a/t)/2), exp(-sqrt(a*s))/s,
453 re(a) > 0, S.Zero, dco), # 4.2.11
454 (besselj(n, a*t), a**n/(sqrt(s**2+a**2)*(s+sqrt(s**2+a**2))**n),
455 re(n) > -1, Abs(im(a)), dco), # 4.14.1
456 (t**b*besselj(n, a*t),
457 2**n/sqrt(pi)*gamma(n+S.Half)*a**n*(s**2+a**2)**(-n-S.Half),
458 And(re(n) > -S.Half, Eq(b, n)), Abs(im(a)), dco), # 4.14.7
459 (t**b*besselj(n, a*t),
460 2**(n+1)/sqrt(pi)*gamma(n+S(3)/2)*a**n*s*(s**2+a**2)**(-n-S(3)/2),
461 And(re(n) > -1, Eq(b, n+1)), Abs(im(a)), dco), # 4.14.8
462 (besselj(0, 2*sqrt(a*t)), exp(-a/s)/s,
463 S.true, S.Zero, dco), # 4.14.25
464 (t**(b)*besselj(n, 2*sqrt(a*t)), a**(n/2)*s**(-n-1)*exp(-a/s),
465 And(re(n) > -1, Eq(b, n*S.Half)), S.Zero, dco), # 4.14.30
466 (besselj(0, a*sqrt(t**2+b*t)),
467 exp(b*s-b*sqrt(s**2+a**2))/sqrt(s**2+a**2),
468 Abs(arg(b)) < pi, Abs(im(a)), dco), # 4.15.19
469 (besseli(n, a*t), a**n/(sqrt(s**2-a**2)*(s+sqrt(s**2-a**2))**n),
470 re(n) > -1, Abs(re(a)), dco), # 4.16.1
471 (t**b*besseli(n, a*t),
472 2**n/sqrt(pi)*gamma(n+S.Half)*a**n*(s**2-a**2)**(-n-S.Half),
473 And(re(n) > -S.Half, Eq(b, n)), Abs(re(a)), dco), # 4.16.6
474 (t**b*besseli(n, a*t),
475 2**(n+1)/sqrt(pi)*gamma(n+S(3)/2)*a**n*s*(s**2-a**2)**(-n-S(3)/2),
476 And(re(n) > -1, Eq(b, n+1)), Abs(re(a)), dco), # 4.16.7
477 (t**(b)*besseli(n, 2*sqrt(a*t)), a**(n/2)*s**(-n-1)*exp(a/s),
478 And(re(n) > -1, Eq(b, n*S.Half)), S.Zero, dco), # 4.16.18
479 (bessely(0, a*t), -2/pi*asinh(s/a)/sqrt(s**2+a**2),
480 S.true, Abs(im(a)), dco), # 4.15.44
481 (besselk(0, a*t), log((s + sqrt(s**2-a**2))/a)/(sqrt(s**2-a**2)),
482 S.true, -re(a), dco) # 4.16.23
483 ]
484 return laplace_transform_rules, t, s
487def _laplace_rule_timescale(f, t, s):
488 """
489 This function applies the time-scaling rule of the Laplace transform in
490 a straight-forward way. For example, if it gets ``(f(a*t), t, s)``, it will
491 compute ``LaplaceTransform(f(t)/a, t, s/a)`` if ``a>0``.
492 """
494 a = Wild('a', exclude=[t])
495 g = WildFunction('g', nargs=1)
496 ma1 = f.match(g)
497 if ma1:
498 arg = ma1[g].args[0].collect(t)
499 ma2 = arg.match(a*t)
500 if ma2 and ma2[a].is_positive and ma2[a] != 1:
501 debug('_laplace_apply_prog rules match:')
502 debugf(' f: %s _ %s, %s )', (f, ma1, ma2))
503 debug(' rule: time scaling (4.1.4)')
504 r, pr, cr = _laplace_transform(1/ma2[a]*ma1[g].func(t),
505 t, s/ma2[a], simplify=False)
506 return (r, pr, cr)
507 return None
510def _laplace_rule_heaviside(f, t, s):
511 """
512 This function deals with time-shifted Heaviside step functions. If the time
513 shift is positive, it applies the time-shift rule of the Laplace transform.
514 For example, if it gets ``(Heaviside(t-a)*f(t), t, s)``, it will compute
515 ``exp(-a*s)*LaplaceTransform(f(t+a), t, s)``.
517 If the time shift is negative, the Heaviside function is simply removed
518 as it means nothing to the Laplace transform.
520 The function does not remove a factor ``Heaviside(t)``; this is done by
521 the simple rules.
522 """
524 a = Wild('a', exclude=[t])
525 y = Wild('y')
526 g = Wild('g')
527 ma1 = f.match(Heaviside(y)*g)
528 if ma1:
529 ma2 = ma1[y].match(t-a)
530 if ma2 and ma2[a].is_positive:
531 debug('_laplace_apply_prog_rules match:')
532 debugf(' f: %s ( %s, %s )', (f, ma1, ma2))
533 debug(' rule: time shift (4.1.4)')
534 r, pr, cr = _laplace_transform(ma1[g].subs(t, t+ma2[a]), t, s,
535 simplify=False)
536 return (exp(-ma2[a]*s)*r, pr, cr)
537 if ma2 and ma2[a].is_negative:
538 debug('_laplace_apply_prog_rules match:')
539 debugf(' f: %s ( %s, %s )', (f, ma1, ma2))
540 debug(' rule: Heaviside factor, negative time shift (4.1.4)')
541 r, pr, cr = _laplace_transform(ma1[g], t, s, simplify=False)
542 return (r, pr, cr)
543 return None
546def _laplace_rule_exp(f, t, s):
547 """
548 If this function finds a factor ``exp(a*t)``, it applies the
549 frequency-shift rule of the Laplace transform and adjusts the convergence
550 plane accordingly. For example, if it gets ``(exp(-a*t)*f(t), t, s)``, it
551 will compute ``LaplaceTransform(f(t), t, s+a)``.
552 """
554 a = Wild('a', exclude=[t])
555 y = Wild('y')
556 z = Wild('z')
557 ma1 = f.match(exp(y)*z)
558 if ma1:
559 ma2 = ma1[y].collect(t).match(a*t)
560 if ma2:
561 debug('_laplace_apply_prog_rules match:')
562 debugf(' f: %s ( %s, %s )', (f, ma1, ma2))
563 debug(' rule: multiply with exp (4.1.5)')
564 r, pr, cr = _laplace_transform(ma1[z], t, s-ma2[a],
565 simplify=False)
566 return (r, pr+re(ma2[a]), cr)
567 return None
570def _laplace_rule_delta(f, t, s):
571 """
572 If this function finds a factor ``DiracDelta(b*t-a)``, it applies the
573 masking property of the delta distribution. For example, if it gets
574 ``(DiracDelta(t-a)*f(t), t, s)``, it will return
575 ``(f(a)*exp(-a*s), -a, True)``.
576 """
577 # This rule is not in Bateman54
579 a = Wild('a', exclude=[t])
580 b = Wild('b', exclude=[t])
582 y = Wild('y')
583 z = Wild('z')
584 ma1 = f.match(DiracDelta(y)*z)
585 if ma1 and not ma1[z].has(DiracDelta):
586 ma2 = ma1[y].collect(t).match(b*t-a)
587 if ma2:
588 debug('_laplace_apply_prog_rules match:')
589 debugf(' f: %s ( %s, %s )', (f, ma1, ma2))
590 debug(' rule: multiply with DiracDelta')
591 loc = ma2[a]/ma2[b]
592 if re(loc) >= 0 and im(loc) == 0:
593 r = exp(-ma2[a]/ma2[b]*s)*ma1[z].subs(t, ma2[a]/ma2[b])/ma2[b]
594 return (r, S.NegativeInfinity, S.true)
595 else:
596 return (0, S.NegativeInfinity, S.true)
597 if ma1[y].is_polynomial(t):
598 ro = roots(ma1[y], t)
599 if roots is not {} and set(ro.values()) == {1}:
600 slope = diff(ma1[y], t)
601 r = Add(
602 *[exp(-x*s)*ma1[z].subs(t, s)/slope.subs(t, x)
603 for x in list(ro.keys()) if im(x) == 0 and re(x) >= 0])
604 return (r, S.NegativeInfinity, S.true)
605 return None
608def _laplace_trig_split(fn):
609 """
610 Helper function for `_laplace_rule_trig`. This function returns two terms
611 `f` and `g`. `f` contains all product terms with sin, cos, sinh, cosh in
612 them; `g` contains everything else.
613 """
614 trigs = [S.One]
615 other = [S.One]
616 for term in Mul.make_args(fn):
617 if term.has(sin, cos, sinh, cosh, exp):
618 trigs.append(term)
619 else:
620 other.append(term)
621 f = Mul(*trigs)
622 g = Mul(*other)
623 return f, g
626def _laplace_trig_expsum(f, t):
627 """
628 Helper function for `_laplace_rule_trig`. This function expects the `f`
629 from `_laplace_trig_split`. It returns two lists `xm` and `xn`. `xm` is
630 a list of dictionaries with keys `k` and `a` representing a function
631 `k*exp(a*t)`. `xn` is a list of all terms that cannot be brought into
632 that form, which may happen, e.g., when a trigonometric function has
633 another function in its argument.
634 """
635 m = Wild('m')
636 p = Wild('p', exclude=[t])
637 xm = []
638 xn = []
640 x1 = f.rewrite(exp).expand()
642 for term in Add.make_args(x1):
643 if not term.has(t):
644 xm.append({'k': term, 'a': 0, re: 0, im: 0})
645 continue
646 term = term.powsimp(combine='exp')
647 if (r := term.match(p*exp(m))) is not None:
648 if (mp := r[m].as_poly(t)) is not None:
649 mc = mp.all_coeffs()
650 if len(mc) == 2:
651 xm.append({
652 'k': r[p]*exp(mc[1]), 'a': mc[0],
653 re: re(mc[0]), im: im(mc[0])})
654 else:
655 xn.append(term)
656 else:
657 xn.append(term)
658 else:
659 xn.append(term)
660 return xm, xn
663def _laplace_trig_ltex(xm, t, s):
664 """
665 Helper function for `_laplace_rule_trig`. This function takes the list of
666 exponentials `xm` from `_laplace_trig_expsum` and simplifies complex
667 conjugate and real symmetric poles. It returns the result as a sum and
668 the convergence plane.
669 """
670 results = []
671 planes = []
673 def _simpc(coeffs):
674 nc = coeffs.copy()
675 for k in range(len(nc)):
676 ri = nc[k].as_real_imag()
677 if ri[0].has(im):
678 nc[k] = nc[k].rewrite(cos)
679 else:
680 nc[k] = (ri[0] + I*ri[1]).rewrite(cos)
681 return nc
683 def _quadpole(t1, k1, k2, k3, s):
684 a, k0, a_r, a_i = t1['a'], t1['k'], t1[re], t1[im]
685 nc = [
686 k0 + k1 + k2 + k3,
687 a*(k0 + k1 - k2 - k3) - 2*I*a_i*k1 + 2*I*a_i*k2,
688 (
689 a**2*(-k0 - k1 - k2 - k3) +
690 a*(4*I*a_i*k0 + 4*I*a_i*k3) +
691 4*a_i**2*k0 + 4*a_i**2*k3),
692 (
693 a**3*(-k0 - k1 + k2 + k3) +
694 a**2*(4*I*a_i*k0 + 2*I*a_i*k1 - 2*I*a_i*k2 - 4*I*a_i*k3) +
695 a*(4*a_i**2*k0 - 4*a_i**2*k3))
696 ]
697 dc = [
698 S.One, S.Zero, 2*a_i**2 - 2*a_r**2,
699 S.Zero, a_i**4 + 2*a_i**2*a_r**2 + a_r**4]
700 n = Add(
701 *[x*s**y for x, y in zip(_simpc(nc), range(len(nc))[::-1])])
702 d = Add(
703 *[x*s**y for x, y in zip(dc, range(len(dc))[::-1])])
704 debugf(' quadpole: (%s) / (%s)', (n, d))
705 return n/d
707 def _ccpole(t1, k1, s):
708 a, k0, a_r, a_i = t1['a'], t1['k'], t1[re], t1[im]
709 nc = [k0 + k1, -a*k0 - a*k1 + 2*I*a_i*k0]
710 dc = [S.One, -2*a_r, a_i**2 + a_r**2]
711 n = Add(
712 *[x*s**y for x, y in zip(_simpc(nc), range(len(nc))[::-1])])
713 d = Add(
714 *[x*s**y for x, y in zip(dc, range(len(dc))[::-1])])
715 debugf(' ccpole: (%s) / (%s)', (n, d))
716 return n/d
718 def _rspole(t1, k2, s):
719 a, k0, a_r, a_i = t1['a'], t1['k'], t1[re], t1[im]
720 nc = [k0 + k2, a*k0 - a*k2 - 2*I*a_i*k0]
721 dc = [S.One, -2*I*a_i, -a_i**2 - a_r**2]
722 n = Add(
723 *[x*s**y for x, y in zip(_simpc(nc), range(len(nc))[::-1])])
724 d = Add(
725 *[x*s**y for x, y in zip(dc, range(len(dc))[::-1])])
726 debugf(' rspole: (%s) / (%s)', (n, d))
727 return n/d
729 def _sypole(t1, k3, s):
730 a, k0 = t1['a'], t1['k']
731 nc = [k0 + k3, a*(k0 - k3)]
732 dc = [S.One, S.Zero, -a**2]
733 n = Add(
734 *[x*s**y for x, y in zip(_simpc(nc), range(len(nc))[::-1])])
735 d = Add(
736 *[x*s**y for x, y in zip(dc, range(len(dc))[::-1])])
737 debugf(' sypole: (%s) / (%s)', (n, d))
738 return n/d
740 def _simplepole(t1, s):
741 a, k0 = t1['a'], t1['k']
742 n = k0
743 d = s - a
744 debugf(' simplepole: (%s) / (%s)', (n, d))
745 return n/d
747 while len(xm) > 0:
748 t1 = xm.pop()
749 i_imagsym = None
750 i_realsym = None
751 i_pointsym = None
752 # The following code checks all remaining poles. If t1 is a pole at
753 # a+b*I, then we check for a-b*I, -a+b*I, and -a-b*I, and
754 # assign the respective indices to i_imagsym, i_realsym, i_pointsym.
755 # -a-b*I / i_pointsym only applies if both a and b are != 0.
756 for i in range(len(xm)):
757 real_eq = t1[re] == xm[i][re]
758 realsym = t1[re] == -xm[i][re]
759 imag_eq = t1[im] == xm[i][im]
760 imagsym = t1[im] == -xm[i][im]
761 if realsym and imagsym and t1[re] != 0 and t1[im] != 0:
762 i_pointsym = i
763 elif realsym and imag_eq and t1[re] != 0:
764 i_realsym = i
765 elif real_eq and imagsym and t1[im] != 0:
766 i_imagsym = i
768 # The next part looks for four possible pole constellations:
769 # quad: a+b*I, a-b*I, -a+b*I, -a-b*I
770 # cc: a+b*I, a-b*I (a may be zero)
771 # quad: a+b*I, -a+b*I (b may be zero)
772 # point: a+b*I, -a-b*I (a!=0 and b!=0 is needed, but that has been
773 # asserted when finding i_pointsym above.)
774 # If none apply, then t1 is a simple pole.
775 if (
776 i_imagsym is not None and i_realsym is not None
777 and i_pointsym is not None):
778 results.append(
779 _quadpole(t1,
780 xm[i_imagsym]['k'], xm[i_realsym]['k'],
781 xm[i_pointsym]['k'], s))
782 planes.append(Abs(re(t1['a'])))
783 # The three additional poles have now been used; to pop them
784 # easily we have to do it from the back.
785 indices_to_pop = [i_imagsym, i_realsym, i_pointsym]
786 indices_to_pop.sort(reverse=True)
787 for i in indices_to_pop:
788 xm.pop(i)
789 elif i_imagsym is not None:
790 results.append(_ccpole(t1, xm[i_imagsym]['k'], s))
791 planes.append(t1[re])
792 xm.pop(i_imagsym)
793 elif i_realsym is not None:
794 results.append(_rspole(t1, xm[i_realsym]['k'], s))
795 planes.append(Abs(t1[re]))
796 xm.pop(i_realsym)
797 elif i_pointsym is not None:
798 results.append(_sypole(t1, xm[i_pointsym]['k'], s))
799 planes.append(Abs(t1[re]))
800 xm.pop(i_pointsym)
801 else:
802 results.append(_simplepole(t1, s))
803 planes.append(t1[re])
805 return Add(*results), Max(*planes)
808def _laplace_rule_trig(fn, t_, s, doit=True, **hints):
809 """
810 This rule covers trigonometric factors by splitting everything into a
811 sum of exponential functions and collecting complex conjugate poles and
812 real symmetric poles.
813 """
814 t = Dummy('t', real=True)
816 if not fn.has(sin, cos, sinh, cosh):
817 return None
819 debugf('_laplace_rule_trig: (%s, %s, %s)', (fn, t_, s))
821 f, g = _laplace_trig_split(fn.subs(t_, t))
822 debugf(' f = %s\n g = %s', (f, g))
824 xm, xn = _laplace_trig_expsum(f, t)
825 debugf(' xm = %s\n xn = %s', (xm, xn))
827 if len(xn) > 0:
828 # not implemented yet
829 debug(' --> xn is not empty; giving up.')
830 return None
832 if not g.has(t):
833 r, p = _laplace_trig_ltex(xm, t, s)
834 return g*r, p, S.true
835 else:
836 # Just transform `g` and make frequency-shifted copies
837 planes = []
838 results = []
839 G, G_plane, G_cond = _laplace_transform(g, t, s)
840 for x1 in xm:
841 results.append(x1['k']*G.subs(s, s-x1['a']))
842 planes.append(G_plane+re(x1['a']))
843 return Add(*results).subs(t, t_), Max(*planes), G_cond
846def _laplace_rule_diff(f, t, s, doit=True, **hints):
847 """
848 This function looks for derivatives in the time domain and replaces it
849 by factors of `s` and initial conditions in the frequency domain. For
850 example, if it gets ``(diff(f(t), t), t, s)``, it will compute
851 ``s*LaplaceTransform(f(t), t, s) - f(0)``.
852 """
854 a = Wild('a', exclude=[t])
855 n = Wild('n', exclude=[t])
856 g = WildFunction('g')
857 ma1 = f.match(a*Derivative(g, (t, n)))
858 if ma1 and ma1[n].is_integer:
859 m = [z.has(t) for z in ma1[g].args]
860 if sum(m) == 1:
861 debug('_laplace_apply_rules match:')
862 debugf(' f, n: %s, %s', (f, ma1[n]))
863 debug(' rule: time derivative (4.1.8)')
864 d = []
865 for k in range(ma1[n]):
866 if k == 0:
867 y = ma1[g].subs(t, 0)
868 else:
869 y = Derivative(ma1[g], (t, k)).subs(t, 0)
870 d.append(s**(ma1[n]-k-1)*y)
871 r, pr, cr = _laplace_transform(ma1[g], t, s, simplify=False)
872 return (ma1[a]*(s**ma1[n]*r - Add(*d)), pr, cr)
873 return None
876def _laplace_rule_sdiff(f, t, s, doit=True, **hints):
877 """
878 This function looks for multiplications with polynoimials in `t` as they
879 correspond to differentiation in the frequency domain. For example, if it
880 gets ``(t*f(t), t, s)``, it will compute
881 ``-Derivative(LaplaceTransform(f(t), t, s), s)``.
882 """
884 if f.is_Mul:
885 pfac = [1]
886 ofac = [1]
887 for fac in Mul.make_args(f):
888 if fac.is_polynomial(t):
889 pfac.append(fac)
890 else:
891 ofac.append(fac)
892 if len(pfac) > 1:
893 pex = prod(pfac)
894 pc = Poly(pex, t).all_coeffs()
895 N = len(pc)
896 if N > 1:
897 debug('_laplace_apply_rules match:')
898 debugf(' f, n: %s, %s', (f, pfac))
899 debug(' rule: frequency derivative (4.1.6)')
900 oex = prod(ofac)
901 r_, p_, c_ = _laplace_transform(oex, t, s, simplify=False)
902 deri = [r_]
903 d1 = False
904 try:
905 d1 = -diff(deri[-1], s)
906 except ValueError:
907 d1 = False
908 if r_.has(LaplaceTransform):
909 for k in range(N-1):
910 deri.append((-1)**(k+1)*Derivative(r_, s, k+1))
911 else:
912 if d1:
913 deri.append(d1)
914 for k in range(N-2):
915 deri.append(-diff(deri[-1], s))
916 if d1:
917 r = Add(*[pc[N-n-1]*deri[n] for n in range(N)])
918 return (r, p_, c_)
919 return None
922def _laplace_expand(f, t, s, doit=True, **hints):
923 """
924 This function tries to expand its argument with successively stronger
925 methods: first it will expand on the top level, then it will expand any
926 multiplications in depth, then it will try all avilable expansion methods,
927 and finally it will try to expand trigonometric functions.
929 If it can expand, it will then compute the Laplace transform of the
930 expanded term.
931 """
933 if f.is_Add:
934 return None
935 r = expand(f, deep=False)
936 if r.is_Add:
937 return _laplace_transform(r, t, s, simplify=False)
938 r = expand_mul(f)
939 if r.is_Add:
940 return _laplace_transform(r, t, s, simplify=False)
941 r = expand(f)
942 if r.is_Add:
943 return _laplace_transform(r, t, s, simplify=False)
944 if r != f:
945 return _laplace_transform(r, t, s, simplify=False)
946 r = expand(expand_trig(f))
947 if r.is_Add:
948 return _laplace_transform(r, t, s, simplify=False)
949 return None
952def _laplace_apply_prog_rules(f, t, s):
953 """
954 This function applies all program rules and returns the result if one
955 of them gives a result.
956 """
958 prog_rules = [_laplace_rule_heaviside, _laplace_rule_delta,
959 _laplace_rule_timescale, _laplace_rule_exp,
960 _laplace_rule_trig,
961 _laplace_rule_diff, _laplace_rule_sdiff]
963 for p_rule in prog_rules:
964 if (L := p_rule(f, t, s)) is not None:
965 return L
966 return None
969def _laplace_apply_simple_rules(f, t, s):
970 """
971 This function applies all simple rules and returns the result if one
972 of them gives a result.
973 """
974 simple_rules, t_, s_ = _laplace_build_rules()
975 prep_old = ''
976 prep_f = ''
977 for t_dom, s_dom, check, plane, prep in simple_rules:
978 if prep_old != prep:
979 prep_f = prep(f.subs({t: t_}))
980 prep_old = prep
981 ma = prep_f.match(t_dom)
982 if ma:
983 try:
984 c = check.xreplace(ma)
985 except TypeError:
986 # This may happen if the time function has imaginary
987 # numbers in it. Then we give up.
988 continue
989 if c == S.true:
990 debug('_laplace_apply_simple_rules match:')
991 debugf(' f: %s', (f,))
992 debugf(' rule: %s o---o %s', (t_dom, s_dom))
993 debugf(' match: %s', (ma, ))
994 return (s_dom.xreplace(ma).subs({s_: s}),
995 plane.xreplace(ma), S.true)
996 return None
999def _laplace_transform(fn, t_, s_, simplify=True):
1000 """
1001 Front-end function of the Laplace transform. It tries to apply all known
1002 rules recursively, and if everything else fails, it tries to integrate.
1003 """
1004 debugf('[LT _l_t] (%s, %s, %s)', (fn, t_, s_))
1006 terms = Add.make_args(fn)
1007 terms_s = []
1008 planes = []
1009 conditions = []
1010 for ff in terms:
1011 k, ft = ff.as_independent(t_, as_Add=False)
1012 if (r := _laplace_apply_simple_rules(ft, t_, s_)) is not None:
1013 pass
1014 elif (r := _laplace_apply_prog_rules(ft, t_, s_)) is not None:
1015 pass
1016 elif (r := _laplace_expand(ft, t_, s_)) is not None:
1017 pass
1018 elif any(undef.has(t_) for undef in ft.atoms(AppliedUndef)):
1019 # If there are undefined functions f(t) then integration is
1020 # unlikely to do anything useful so we skip it and given an
1021 # unevaluated LaplaceTransform.
1022 r = (LaplaceTransform(ft, t_, s_), S.NegativeInfinity, True)
1023 elif (r := _laplace_transform_integration(
1024 ft, t_, s_, simplify=simplify)) is not None:
1025 pass
1026 else:
1027 r = (LaplaceTransform(ft, t_, s_), S.NegativeInfinity, True)
1028 (ri_, pi_, ci_) = r
1029 terms_s.append(k*ri_)
1030 planes.append(pi_)
1031 conditions.append(ci_)
1033 result = Add(*terms_s)
1034 if simplify:
1035 result = result.simplify(doit=False)
1036 plane = Max(*planes)
1037 condition = And(*conditions)
1039 return result, plane, condition
1042class LaplaceTransform(IntegralTransform):
1043 """
1044 Class representing unevaluated Laplace transforms.
1046 For usage of this class, see the :class:`IntegralTransform` docstring.
1048 For how to compute Laplace transforms, see the :func:`laplace_transform`
1049 docstring.
1051 If this is called with ``.doit()``, it returns the Laplace transform as an
1052 expression. If it is called with ``.doit(noconds=False)``, it returns a
1053 tuple containing the same expression, a convergence plane, and conditions.
1054 """
1056 _name = 'Laplace'
1058 def _compute_transform(self, f, t, s, **hints):
1059 _simplify = hints.get('simplify', False)
1060 LT = _laplace_transform_integration(f, t, s, simplify=_simplify)
1061 return LT
1063 def _as_integral(self, f, t, s):
1064 return Integral(f*exp(-s*t), (t, S.Zero, S.Infinity))
1066 def _collapse_extra(self, extra):
1067 conds = []
1068 planes = []
1069 for plane, cond in extra:
1070 conds.append(cond)
1071 planes.append(plane)
1072 cond = And(*conds)
1073 plane = Max(*planes)
1074 if cond == S.false:
1075 raise IntegralTransformError(
1076 'Laplace', None, 'No combined convergence.')
1077 return plane, cond
1079 def doit(self, **hints):
1080 """
1081 Try to evaluate the transform in closed form.
1083 Explanation
1084 ===========
1086 Standard hints are the following:
1087 - ``noconds``: if True, do not return convergence conditions. The
1088 default setting is `True`.
1089 - ``simplify``: if True, it simplifies the final result. The
1090 default setting is `False`.
1091 """
1092 _noconds = hints.get('noconds', True)
1093 _simplify = hints.get('simplify', False)
1095 debugf('[LT doit] (%s, %s, %s)', (self.function,
1096 self.function_variable,
1097 self.transform_variable))
1099 t_ = self.function_variable
1100 s_ = self.transform_variable
1101 fn = self.function
1103 r = _laplace_transform(fn, t_, s_, simplify=_simplify)
1105 if _noconds:
1106 return r[0]
1107 else:
1108 return r
1111def laplace_transform(f, t, s, legacy_matrix=True, **hints):
1112 r"""
1113 Compute the Laplace Transform `F(s)` of `f(t)`,
1115 .. math :: F(s) = \int_{0^{-}}^\infty e^{-st} f(t) \mathrm{d}t.
1117 Explanation
1118 ===========
1120 For all sensible functions, this converges absolutely in a
1121 half-plane
1123 .. math :: a < \operatorname{Re}(s)
1125 This function returns ``(F, a, cond)`` where ``F`` is the Laplace
1126 transform of ``f``, `a` is the half-plane of convergence, and `cond` are
1127 auxiliary convergence conditions.
1129 The implementation is rule-based, and if you are interested in which
1130 rules are applied, and whether integration is attempted, you can switch
1131 debug information on by setting ``sympy.SYMPY_DEBUG=True``. The numbers
1132 of the rules in the debug information (and the code) refer to Bateman's
1133 Tables of Integral Transforms [1].
1135 The lower bound is `0-`, meaning that this bound should be approached
1136 from the lower side. This is only necessary if distributions are involved.
1137 At present, it is only done if `f(t)` contains ``DiracDelta``, in which
1138 case the Laplace transform is computed implicitly as
1140 .. math ::
1141 F(s) = \lim_{\tau\to 0^{-}} \int_{\tau}^\infty e^{-st}
1142 f(t) \mathrm{d}t
1144 by applying rules.
1146 If the Laplace transform cannot be fully computed in closed form, this
1147 function returns expressions containing unevaluated
1148 :class:`LaplaceTransform` objects.
1150 For a description of possible hints, refer to the docstring of
1151 :func:`sympy.integrals.transforms.IntegralTransform.doit`. If
1152 ``noconds=True``, only `F` will be returned (i.e. not ``cond``, and also
1153 not the plane ``a``).
1155 .. deprecated:: 1.9
1156 Legacy behavior for matrices where ``laplace_transform`` with
1157 ``noconds=False`` (the default) returns a Matrix whose elements are
1158 tuples. The behavior of ``laplace_transform`` for matrices will change
1159 in a future release of SymPy to return a tuple of the transformed
1160 Matrix and the convergence conditions for the matrix as a whole. Use
1161 ``legacy_matrix=False`` to enable the new behavior.
1163 Examples
1164 ========
1166 >>> from sympy import DiracDelta, exp, laplace_transform
1167 >>> from sympy.abc import t, s, a
1168 >>> laplace_transform(t**4, t, s)
1169 (24/s**5, 0, True)
1170 >>> laplace_transform(t**a, t, s)
1171 (gamma(a + 1)/(s*s**a), 0, re(a) > -1)
1172 >>> laplace_transform(DiracDelta(t)-a*exp(-a*t), t, s, simplify=True)
1173 (s/(a + s), -re(a), True)
1175 References
1176 ==========
1178 .. [1] Erdelyi, A. (ed.), Tables of Integral Transforms, Volume 1,
1179 Bateman Manuscript Prooject, McGraw-Hill (1954), available:
1180 https://resolver.caltech.edu/CaltechAUTHORS:20140123-101456353
1182 See Also
1183 ========
1185 inverse_laplace_transform, mellin_transform, fourier_transform
1186 hankel_transform, inverse_hankel_transform
1188 """
1190 _noconds = hints.get('noconds', False)
1191 _simplify = hints.get('simplify', False)
1193 if isinstance(f, MatrixBase) and hasattr(f, 'applyfunc'):
1195 conds = not hints.get('noconds', False)
1197 if conds and legacy_matrix:
1198 adt = 'deprecated-laplace-transform-matrix'
1199 sympy_deprecation_warning(
1200 """
1201Calling laplace_transform() on a Matrix with noconds=False (the default) is
1202deprecated. Either noconds=True or use legacy_matrix=False to get the new
1203behavior.
1204 """,
1205 deprecated_since_version='1.9',
1206 active_deprecations_target=adt,
1207 )
1208 # Temporarily disable the deprecation warning for non-Expr objects
1209 # in Matrix
1210 with ignore_warnings(SymPyDeprecationWarning):
1211 return f.applyfunc(
1212 lambda fij: laplace_transform(fij, t, s, **hints))
1213 else:
1214 elements_trans = [laplace_transform(
1215 fij, t, s, **hints) for fij in f]
1216 if conds:
1217 elements, avals, conditions = zip(*elements_trans)
1218 f_laplace = type(f)(*f.shape, elements)
1219 return f_laplace, Max(*avals), And(*conditions)
1220 else:
1221 return type(f)(*f.shape, elements_trans)
1223 LT = LaplaceTransform(f, t, s).doit(noconds=False, simplify=_simplify)
1225 if not _noconds:
1226 return LT
1227 else:
1228 return LT[0]
1231def _inverse_laplace_transform_integration(F, s, t_, plane, simplify=True):
1232 """ The backend function for inverse Laplace transforms. """
1233 from sympy.integrals.meijerint import meijerint_inversion, _get_coeff_exp
1234 from sympy.integrals.transforms import inverse_mellin_transform
1236 # There are two strategies we can try:
1237 # 1) Use inverse mellin transform, related by a simple change of variables.
1238 # 2) Use the inversion integral.
1240 t = Dummy('t', real=True)
1242 def pw_simp(*args):
1243 """ Simplify a piecewise expression from hyperexpand. """
1244 # XXX we break modularity here!
1245 if len(args) != 3:
1246 return Piecewise(*args)
1247 arg = args[2].args[0].argument
1248 coeff, exponent = _get_coeff_exp(arg, t)
1249 e1 = args[0].args[0]
1250 e2 = args[1].args[0]
1251 return (
1252 Heaviside(1/Abs(coeff) - t**exponent)*e1 +
1253 Heaviside(t**exponent - 1/Abs(coeff))*e2)
1255 if F.is_rational_function(s):
1256 F = F.apart(s)
1258 if F.is_Add:
1259 f = Add(
1260 *[_inverse_laplace_transform_integration(X, s, t, plane, simplify)
1261 for X in F.args])
1262 return _simplify(f.subs(t, t_), simplify), True
1264 try:
1265 f, cond = inverse_mellin_transform(F, s, exp(-t), (None, S.Infinity),
1266 needeval=True, noconds=False)
1267 except IntegralTransformError:
1268 f = None
1269 if f is None:
1270 f = meijerint_inversion(F, s, t)
1271 if f is None:
1272 return None
1273 if f.is_Piecewise:
1274 f, cond = f.args[0]
1275 if f.has(Integral):
1276 return None
1277 else:
1278 cond = S.true
1279 f = f.replace(Piecewise, pw_simp)
1281 if f.is_Piecewise:
1282 # many of the functions called below can't work with piecewise
1283 # (b/c it has a bool in args)
1284 return f.subs(t, t_), cond
1286 u = Dummy('u')
1288 def simp_heaviside(arg, H0=S.Half):
1289 a = arg.subs(exp(-t), u)
1290 if a.has(t):
1291 return Heaviside(arg, H0)
1292 from sympy.solvers.inequalities import _solve_inequality
1293 rel = _solve_inequality(a > 0, u)
1294 if rel.lts == u:
1295 k = log(rel.gts)
1296 return Heaviside(t + k, H0)
1297 else:
1298 k = log(rel.lts)
1299 return Heaviside(-(t + k), H0)
1301 f = f.replace(Heaviside, simp_heaviside)
1303 def simp_exp(arg):
1304 return expand_complex(exp(arg))
1306 f = f.replace(exp, simp_exp)
1308 # TODO it would be nice to fix cosh and sinh ... simplify messes these
1309 # exponentials up
1311 return _simplify(f.subs(t, t_), simplify), cond
1314def _complete_the_square_in_denom(f, s):
1315 from sympy.simplify.radsimp import fraction
1316 [n, d] = fraction(f)
1317 if d.is_polynomial(s):
1318 cf = d.as_poly(s).all_coeffs()
1319 if len(cf) == 3:
1320 a, b, c = cf
1321 d = a*((s+b/(2*a))**2+c/a-(b/(2*a))**2)
1322 return n/d
1325@cacheit
1326def _inverse_laplace_build_rules():
1327 """
1328 This is an internal helper function that returns the table of inverse
1329 Laplace transform rules in terms of the time variable `t` and the
1330 frequency variable `s`. It is used by `_inverse_laplace_apply_rules`.
1331 """
1332 s = Dummy('s')
1333 t = Dummy('t')
1334 a = Wild('a', exclude=[s])
1335 b = Wild('b', exclude=[s])
1336 c = Wild('c', exclude=[s])
1338 debug('_inverse_laplace_build_rules is building rules')
1340 def _frac(f, s):
1341 try:
1342 return f.factor(s)
1343 except PolynomialError:
1344 return f
1346 def same(f): return f
1347 # This list is sorted according to the prep function needed.
1348 _ILT_rules = [
1349 (a/s, a, S.true, same, 1),
1350 (b*(s+a)**(-c), t**(c-1)*exp(-a*t)/gamma(c), c > 0, same, 1),
1351 (1/(s**2+a**2)**2, (sin(a*t) - a*t*cos(a*t))/(2*a**3),
1352 S.true, same, 1),
1353 # The next two rules must be there in that order. For the second
1354 # one, the condition would be a != 0 or, respectively, to take the
1355 # limit a -> 0 after the transform if a == 0. It is much simpler if
1356 # the case a == 0 has its own rule.
1357 (1/(s**b), t**(b - 1)/gamma(b), S.true, same, 1),
1358 (1/(s*(s+a)**b), lowergamma(b, a*t)/(a**b*gamma(b)),
1359 S.true, same, 1)
1360 ]
1361 return _ILT_rules, s, t
1364def _inverse_laplace_apply_simple_rules(f, s, t):
1365 """
1366 Helper function for the class InverseLaplaceTransform.
1367 """
1368 if f == 1:
1369 debug('_inverse_laplace_apply_simple_rules match:')
1370 debugf(' f: %s', (1,))
1371 debugf(' rule: 1 o---o DiracDelta(%s)', (t,))
1372 return DiracDelta(t), S.true
1374 _ILT_rules, s_, t_ = _inverse_laplace_build_rules()
1375 _prep = ''
1376 fsubs = f.subs({s: s_})
1378 for s_dom, t_dom, check, prep, fac in _ILT_rules:
1379 if _prep != (prep, fac):
1380 _F = prep(fsubs*fac)
1381 _prep = (prep, fac)
1382 ma = _F.match(s_dom)
1383 if ma:
1384 try:
1385 c = check.xreplace(ma)
1386 except TypeError:
1387 continue
1388 if c == S.true:
1389 debug('_inverse_laplace_apply_simple_rules match:')
1390 debugf(' f: %s', (f,))
1391 debugf(' rule: %s o---o %s', (s_dom, t_dom))
1392 debugf(' ma: %s', (ma,))
1393 return Heaviside(t)*t_dom.xreplace(ma).subs({t_: t}), S.true
1395 return None
1398def _inverse_laplace_time_shift(F, s, t, plane):
1399 """
1400 Helper function for the class InverseLaplaceTransform.
1401 """
1402 a = Wild('a', exclude=[s])
1403 g = Wild('g')
1405 if not F.has(s):
1406 return F*DiracDelta(t), S.true
1407 ma1 = F.match(exp(a*s))
1408 if ma1:
1409 if ma1[a].is_negative:
1410 debug('_inverse_laplace_time_shift match:')
1411 debugf(' f: %s', (F,))
1412 debug(' rule: exp(-a*s) o---o DiracDelta(t-a)')
1413 debugf(' ma: %s', (ma1,))
1414 return DiracDelta(t+ma1[a]), S.true
1415 else:
1416 debug('_inverse_laplace_time_shift match: negative time shift')
1417 return InverseLaplaceTransform(F, s, t, plane), S.true
1419 ma1 = F.match(exp(a*s)*g)
1420 if ma1:
1421 if ma1[a].is_negative:
1422 debug('_inverse_laplace_time_shift match:')
1423 debugf(' f: %s', (F,))
1424 debug(' rule: exp(-a*s)*F(s) o---o Heaviside(t-a)*f(t-a)')
1425 debugf(' ma: %s', (ma1,))
1426 return _inverse_laplace_transform(ma1[g], s, t+ma1[a], plane)
1427 else:
1428 debug('_inverse_laplace_time_shift match: negative time shift')
1429 return InverseLaplaceTransform(F, s, t, plane), S.true
1430 return None
1433def _inverse_laplace_time_diff(F, s, t, plane):
1434 """
1435 Helper function for the class InverseLaplaceTransform.
1436 """
1437 n = Wild('n', exclude=[s])
1438 g = Wild('g')
1440 ma1 = F.match(s**n*g)
1441 if ma1 and ma1[n].is_integer and ma1[n].is_positive:
1442 debug('_inverse_laplace_time_diff match:')
1443 debugf(' f: %s', (F,))
1444 debug(' rule: s**n*F(s) o---o diff(f(t), t, n)')
1445 debugf(' ma: %s', (ma1,))
1446 r, c = _inverse_laplace_transform(ma1[g], s, t, plane)
1447 r = r.replace(Heaviside(t), 1)
1448 if r.has(InverseLaplaceTransform):
1449 return diff(r, t, ma1[n]), c
1450 else:
1451 return Heaviside(t)*diff(r, t, ma1[n]), c
1452 return None
1455def _inverse_laplace_apply_prog_rules(F, s, t, plane):
1456 """
1457 Helper function for the class InverseLaplaceTransform.
1458 """
1459 prog_rules = [_inverse_laplace_time_shift,
1460 _inverse_laplace_time_diff]
1462 for p_rule in prog_rules:
1463 if (r := p_rule(F, s, t, plane)) is not None:
1464 return r
1465 return None
1468def _inverse_laplace_expand(fn, s, t, plane):
1469 """
1470 Helper function for the class InverseLaplaceTransform.
1471 """
1472 if fn.is_Add:
1473 return None
1474 r = expand(fn, deep=False)
1475 if r.is_Add:
1476 return _inverse_laplace_transform(r, s, t, plane)
1477 r = expand_mul(fn)
1478 if r.is_Add:
1479 return _inverse_laplace_transform(r, s, t, plane)
1480 r = expand(fn)
1481 if r.is_Add:
1482 return _inverse_laplace_transform(r, s, t, plane)
1483 if fn.is_rational_function(s):
1484 r = fn.apart(s).doit()
1485 if r.is_Add:
1486 return _inverse_laplace_transform(r, s, t, plane)
1487 return None
1490def _inverse_laplace_rational(fn, s, t, plane, simplify):
1491 """
1492 Helper function for the class InverseLaplaceTransform.
1493 """
1494 debugf('[ILT _i_l_r] (%s, %s, %s)', (fn, s, t))
1495 x_ = symbols('x_')
1496 f = fn.apart(s)
1497 terms = Add.make_args(f)
1498 terms_t = []
1499 conditions = [S.true]
1500 for term in terms:
1501 [n, d] = term.as_numer_denom()
1502 dc = d.as_poly(s).all_coeffs()
1503 dc_lead = dc[0]
1504 dc = [x/dc_lead for x in dc]
1505 nc = [x/dc_lead for x in n.as_poly(s).all_coeffs()]
1506 if len(dc) == 1:
1507 r = nc[0]*DiracDelta(t)
1508 terms_t.append(r)
1509 elif len(dc) == 2:
1510 r = nc[0]*exp(-dc[1]*t)
1511 terms_t.append(Heaviside(t)*r)
1512 elif len(dc) == 3:
1513 a = dc[1]/2
1514 b = (dc[2]-a**2).factor()
1515 if len(nc) == 1:
1516 nc = [S.Zero] + nc
1517 l, m = tuple(nc)
1518 if b == 0:
1519 r = (m*t+l*(1-a*t))*exp(-a*t)
1520 else:
1521 hyp = False
1522 if b.is_negative:
1523 b = -b
1524 hyp = True
1525 b2 = list(roots(x_**2-b, x_).keys())[0]
1526 bs = sqrt(b).simplify()
1527 if hyp:
1528 r = (
1529 l*exp(-a*t)*cosh(b2*t) + (m-a*l) /
1530 bs*exp(-a*t)*sinh(bs*t))
1531 else:
1532 r = l*exp(-a*t)*cos(b2*t) + (m-a*l)/bs*exp(-a*t)*sin(bs*t)
1533 terms_t.append(Heaviside(t)*r)
1534 else:
1535 ft, cond = _inverse_laplace_transform(
1536 fn, s, t, plane, simplify=True, dorational=False)
1537 terms_t.append(ft)
1538 conditions.append(cond)
1540 result = Add(*terms_t)
1541 if simplify:
1542 result = result.simplify(doit=False)
1543 debugf('[ILT _i_l_r] returns %s', (result,))
1544 return result, And(*conditions)
1547def _inverse_laplace_transform(
1548 fn, s_, t_, plane, simplify=True, dorational=True):
1549 """
1550 Front-end function of the inverse Laplace transform. It tries to apply all
1551 known rules recursively. If everything else fails, it tries to integrate.
1552 """
1553 terms = Add.make_args(fn)
1554 terms_t = []
1555 conditions = []
1557 debugf('[ILT _i_l_t] (%s, %s, %s)', (fn, s_, t_))
1559 for term in terms:
1560 k, f = term.as_independent(s_, as_Add=False)
1561 if (
1562 dorational and term.is_rational_function(s_) and
1563 (
1564 r := _inverse_laplace_rational(
1565 f, s_, t_, plane, simplify)) is not None):
1566 pass
1567 elif (r := _inverse_laplace_apply_simple_rules(f, s_, t_)) is not None:
1568 pass
1569 elif (r := _inverse_laplace_expand(f, s_, t_, plane)) is not None:
1570 pass
1571 elif (
1572 (r := _inverse_laplace_apply_prog_rules(f, s_, t_, plane))
1573 is not None):
1574 pass
1575 elif any(undef.has(s_) for undef in f.atoms(AppliedUndef)):
1576 # If there are undefined functions f(t) then integration is
1577 # unlikely to do anything useful so we skip it and given an
1578 # unevaluated LaplaceTransform.
1579 r = (InverseLaplaceTransform(f, s_, t_, plane), S.true)
1580 elif (
1581 r := _inverse_laplace_transform_integration(
1582 f, s_, t_, plane, simplify=simplify)) is not None:
1583 pass
1584 else:
1585 r = (InverseLaplaceTransform(f, s_, t_, plane), S.true)
1586 (ri_, ci_) = r
1587 terms_t.append(k*ri_)
1588 conditions.append(ci_)
1590 result = Add(*terms_t)
1591 if simplify:
1592 result = result.simplify(doit=False)
1593 condition = And(*conditions)
1595 return result, condition
1598class InverseLaplaceTransform(IntegralTransform):
1599 """
1600 Class representing unevaluated inverse Laplace transforms.
1602 For usage of this class, see the :class:`IntegralTransform` docstring.
1604 For how to compute inverse Laplace transforms, see the
1605 :func:`inverse_laplace_transform` docstring.
1606 """
1608 _name = 'Inverse Laplace'
1609 _none_sentinel = Dummy('None')
1610 _c = Dummy('c')
1612 def __new__(cls, F, s, x, plane, **opts):
1613 if plane is None:
1614 plane = InverseLaplaceTransform._none_sentinel
1615 return IntegralTransform.__new__(cls, F, s, x, plane, **opts)
1617 @property
1618 def fundamental_plane(self):
1619 plane = self.args[3]
1620 if plane is InverseLaplaceTransform._none_sentinel:
1621 plane = None
1622 return plane
1624 def _compute_transform(self, F, s, t, **hints):
1625 return _inverse_laplace_transform_integration(
1626 F, s, t, self.fundamental_plane, **hints)
1628 def _as_integral(self, F, s, t):
1629 c = self.__class__._c
1630 return (
1631 Integral(exp(s*t)*F, (s, c - S.ImaginaryUnit*S.Infinity,
1632 c + S.ImaginaryUnit*S.Infinity)) /
1633 (2*S.Pi*S.ImaginaryUnit))
1635 def doit(self, **hints):
1636 """
1637 Try to evaluate the transform in closed form.
1639 Explanation
1640 ===========
1642 Standard hints are the following:
1643 - ``noconds``: if True, do not return convergence conditions. The
1644 default setting is `True`.
1645 - ``simplify``: if True, it simplifies the final result. The
1646 default setting is `False`.
1647 """
1648 _noconds = hints.get('noconds', True)
1649 _simplify = hints.get('simplify', False)
1651 debugf('[ILT doit] (%s, %s, %s)', (self.function,
1652 self.function_variable,
1653 self.transform_variable))
1655 s_ = self.function_variable
1656 t_ = self.transform_variable
1657 fn = self.function
1658 plane = self.fundamental_plane
1660 r = _inverse_laplace_transform(fn, s_, t_, plane, simplify=_simplify)
1662 if _noconds:
1663 return r[0]
1664 else:
1665 return r
1668def inverse_laplace_transform(F, s, t, plane=None, **hints):
1669 r"""
1670 Compute the inverse Laplace transform of `F(s)`, defined as
1672 .. math ::
1673 f(t) = \frac{1}{2\pi i} \int_{c-i\infty}^{c+i\infty} e^{st}
1674 F(s) \mathrm{d}s,
1676 for `c` so large that `F(s)` has no singularites in the
1677 half-plane `\operatorname{Re}(s) > c-\epsilon`.
1679 Explanation
1680 ===========
1682 The plane can be specified by
1683 argument ``plane``, but will be inferred if passed as None.
1685 Under certain regularity conditions, this recovers `f(t)` from its
1686 Laplace Transform `F(s)`, for non-negative `t`, and vice
1687 versa.
1689 If the integral cannot be computed in closed form, this function returns
1690 an unevaluated :class:`InverseLaplaceTransform` object.
1692 Note that this function will always assume `t` to be real,
1693 regardless of the SymPy assumption on `t`.
1695 For a description of possible hints, refer to the docstring of
1696 :func:`sympy.integrals.transforms.IntegralTransform.doit`.
1698 Examples
1699 ========
1701 >>> from sympy import inverse_laplace_transform, exp, Symbol
1702 >>> from sympy.abc import s, t
1703 >>> a = Symbol('a', positive=True)
1704 >>> inverse_laplace_transform(exp(-a*s)/s, s, t)
1705 Heaviside(-a + t)
1707 See Also
1708 ========
1710 laplace_transform
1711 hankel_transform, inverse_hankel_transform
1712 """
1713 if isinstance(F, MatrixBase) and hasattr(F, 'applyfunc'):
1714 return F.applyfunc(
1715 lambda Fij: inverse_laplace_transform(Fij, s, t, plane, **hints))
1716 return InverseLaplaceTransform(F, s, t, plane).doit(**hints)
1719def _fast_inverse_laplace(e, s, t):
1720 """Fast inverse Laplace transform of rational function including RootSum"""
1721 a, b, n = symbols('a, b, n', cls=Wild, exclude=[s])
1723 def _ilt(e):
1724 if not e.has(s):
1725 return e
1726 elif e.is_Add:
1727 return _ilt_add(e)
1728 elif e.is_Mul:
1729 return _ilt_mul(e)
1730 elif e.is_Pow:
1731 return _ilt_pow(e)
1732 elif isinstance(e, RootSum):
1733 return _ilt_rootsum(e)
1734 else:
1735 raise NotImplementedError
1737 def _ilt_add(e):
1738 return e.func(*map(_ilt, e.args))
1740 def _ilt_mul(e):
1741 coeff, expr = e.as_independent(s)
1742 if expr.is_Mul:
1743 raise NotImplementedError
1744 return coeff * _ilt(expr)
1746 def _ilt_pow(e):
1747 match = e.match((a*s + b)**n)
1748 if match is not None:
1749 nm, am, bm = match[n], match[a], match[b]
1750 if nm.is_Integer and nm < 0:
1751 return t**(-nm-1)*exp(-(bm/am)*t)/(am**-nm*gamma(-nm))
1752 if nm == 1:
1753 return exp(-(bm/am)*t) / am
1754 raise NotImplementedError
1756 def _ilt_rootsum(e):
1757 expr = e.fun.expr
1758 [variable] = e.fun.variables
1759 return RootSum(e.poly, Lambda(variable, together(_ilt(expr))))
1761 return _ilt(e)