Coverage for /usr/lib/python3/dist-packages/mpmath/calculus/inverselaplace.py: 15%
170 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1# contributed to mpmath by Kristopher L. Kuhlman, February 2017
3class InverseLaplaceTransform(object):
4 r"""
5 Inverse Laplace transform methods are implemented using this
6 class, in order to simplify the code and provide a common
7 infrastructure.
9 Implement a custom inverse Laplace transform algorithm by
10 subclassing :class:`InverseLaplaceTransform` and implementing the
11 appropriate methods. The subclass can then be used by
12 :func:`~mpmath.invertlaplace` by passing it as the *method*
13 argument.
14 """
16 def __init__(self,ctx):
17 self.ctx = ctx
19 def calc_laplace_parameter(self,t,**kwargs):
20 r"""
21 Determine the vector of Laplace parameter values needed for an
22 algorithm, this will depend on the choice of algorithm (de
23 Hoog is default), the algorithm-specific parameters passed (or
24 default ones), and desired time.
25 """
26 raise NotImplementedError
28 def calc_time_domain_solution(self,fp):
29 r"""
30 Compute the time domain solution, after computing the
31 Laplace-space function evaluations at the abscissa required
32 for the algorithm. Abscissa computed for one algorithm are
33 typically not useful for another algorithm.
34 """
35 raise NotImplementedError
37class FixedTalbot(InverseLaplaceTransform):
39 def calc_laplace_parameter(self,t,**kwargs):
40 r"""The "fixed" Talbot method deforms the Bromwich contour towards
41 `-\infty` in the shape of a parabola. Traditionally the Talbot
42 algorithm has adjustable parameters, but the "fixed" version
43 does not. The `r` parameter could be passed in as a parameter,
44 if you want to override the default given by (Abate & Valko,
45 2004).
47 The Laplace parameter is sampled along a parabola opening
48 along the negative imaginary axis, with the base of the
49 parabola along the real axis at
50 `p=\frac{r}{t_\mathrm{max}}`. As the number of terms used in
51 the approximation (degree) grows, the abscissa required for
52 function evaluation tend towards `-\infty`, requiring high
53 precision to prevent overflow. If any poles, branch cuts or
54 other singularities exist such that the deformed Bromwich
55 contour lies to the left of the singularity, the method will
56 fail.
58 **Optional arguments**
60 :class:`~mpmath.calculus.inverselaplace.FixedTalbot.calc_laplace_parameter`
61 recognizes the following keywords
63 *tmax*
64 maximum time associated with vector of times
65 (typically just the time requested)
66 *degree*
67 integer order of approximation (M = number of terms)
68 *r*
69 abscissa for `p_0` (otherwise computed using rule
70 of thumb `2M/5`)
72 The working precision will be increased according to a rule of
73 thumb. If 'degree' is not specified, the working precision and
74 degree are chosen to hopefully achieve the dps of the calling
75 context. If 'degree' is specified, the working precision is
76 chosen to achieve maximum resulting precision for the
77 specified degree.
79 .. math ::
81 p_0=\frac{r}{t}
83 .. math ::
85 p_i=\frac{i r \pi}{Mt_\mathrm{max}}\left[\cot\left(
86 \frac{i\pi}{M}\right) + j \right] \qquad 1\le i <M
88 where `j=\sqrt{-1}`, `r=2M/5`, and `t_\mathrm{max}` is the
89 maximum specified time.
91 """
93 # required
94 # ------------------------------
95 # time of desired approximation
96 self.t = self.ctx.convert(t)
98 # optional
99 # ------------------------------
100 # maximum time desired (used for scaling) default is requested
101 # time.
102 self.tmax = self.ctx.convert(kwargs.get('tmax',self.t))
104 # empirical relationships used here based on a linear fit of
105 # requested and delivered dps for exponentially decaying time
106 # functions for requested dps up to 512.
108 if 'degree' in kwargs:
109 self.degree = kwargs['degree']
110 self.dps_goal = self.degree
111 else:
112 self.dps_goal = int(1.72*self.ctx.dps)
113 self.degree = max(12,int(1.38*self.dps_goal))
115 M = self.degree
117 # this is adjusting the dps of the calling context hopefully
118 # the caller doesn't monkey around with it between calling
119 # this routine and calc_time_domain_solution()
120 self.dps_orig = self.ctx.dps
121 self.ctx.dps = self.dps_goal
123 # Abate & Valko rule of thumb for r parameter
124 self.r = kwargs.get('r',self.ctx.fraction(2,5)*M)
126 self.theta = self.ctx.linspace(0.0, self.ctx.pi, M+1)
128 self.cot_theta = self.ctx.matrix(M,1)
129 self.cot_theta[0] = 0 # not used
131 # all but time-dependent part of p
132 self.delta = self.ctx.matrix(M,1)
133 self.delta[0] = self.r
135 for i in range(1,M):
136 self.cot_theta[i] = self.ctx.cot(self.theta[i])
137 self.delta[i] = self.r*self.theta[i]*(self.cot_theta[i] + 1j)
139 self.p = self.ctx.matrix(M,1)
140 self.p = self.delta/self.tmax
142 # NB: p is complex (mpc)
144 def calc_time_domain_solution(self,fp,t,manual_prec=False):
145 r"""The fixed Talbot time-domain solution is computed from the
146 Laplace-space function evaluations using
148 .. math ::
150 f(t,M)=\frac{2}{5t}\sum_{k=0}^{M-1}\Re \left[
151 \gamma_k \bar{f}(p_k)\right]
153 where
155 .. math ::
157 \gamma_0 = \frac{1}{2}e^{r}\bar{f}(p_0)
159 .. math ::
161 \gamma_k = e^{tp_k}\left\lbrace 1 + \frac{jk\pi}{M}\left[1 +
162 \cot \left( \frac{k \pi}{M} \right)^2 \right] - j\cot\left(
163 \frac{k \pi}{M}\right)\right \rbrace \qquad 1\le k<M.
165 Again, `j=\sqrt{-1}`.
167 Before calling this function, call
168 :class:`~mpmath.calculus.inverselaplace.FixedTalbot.calc_laplace_parameter`
169 to set the parameters and compute the required coefficients.
171 **References**
173 1. Abate, J., P. Valko (2004). Multi-precision Laplace
174 transform inversion. *International Journal for Numerical
175 Methods in Engineering* 60:979-993,
176 http://dx.doi.org/10.1002/nme.995
177 2. Talbot, A. (1979). The accurate numerical inversion of
178 Laplace transforms. *IMA Journal of Applied Mathematics*
179 23(1):97, http://dx.doi.org/10.1093/imamat/23.1.97
180 """
182 # required
183 # ------------------------------
184 self.t = self.ctx.convert(t)
186 # assume fp was computed from p matrix returned from
187 # calc_laplace_parameter(), so is already a list or matrix of
188 # mpmath 'mpc' types
190 # these were computed in previous call to
191 # calc_laplace_parameter()
192 theta = self.theta
193 delta = self.delta
194 M = self.degree
195 p = self.p
196 r = self.r
198 ans = self.ctx.matrix(M,1)
199 ans[0] = self.ctx.exp(delta[0])*fp[0]/2
201 for i in range(1,M):
202 ans[i] = self.ctx.exp(delta[i])*fp[i]*(
203 1 + 1j*theta[i]*(1 + self.cot_theta[i]**2) -
204 1j*self.cot_theta[i])
206 result = self.ctx.fraction(2,5)*self.ctx.fsum(ans)/self.t
208 # setting dps back to value when calc_laplace_parameter was
209 # called, unless flag is set.
210 if not manual_prec:
211 self.ctx.dps = self.dps_orig
213 return result.real
215# ****************************************
217class Stehfest(InverseLaplaceTransform):
219 def calc_laplace_parameter(self,t,**kwargs):
220 r"""
221 The Gaver-Stehfest method is a discrete approximation of the
222 Widder-Post inversion algorithm, rather than a direct
223 approximation of the Bromwich contour integral.
225 The method abscissa along the real axis, and therefore has
226 issues inverting oscillatory functions (which have poles in
227 pairs away from the real axis).
229 The working precision will be increased according to a rule of
230 thumb. If 'degree' is not specified, the working precision and
231 degree are chosen to hopefully achieve the dps of the calling
232 context. If 'degree' is specified, the working precision is
233 chosen to achieve maximum resulting precision for the
234 specified degree.
236 .. math ::
238 p_k = \frac{k \log 2}{t} \qquad 1 \le k \le M
239 """
241 # required
242 # ------------------------------
243 # time of desired approximation
244 self.t = self.ctx.convert(t)
246 # optional
247 # ------------------------------
249 # empirical relationships used here based on a linear fit of
250 # requested and delivered dps for exponentially decaying time
251 # functions for requested dps up to 512.
253 if 'degree' in kwargs:
254 self.degree = kwargs['degree']
255 self.dps_goal = int(1.38*self.degree)
256 else:
257 self.dps_goal = int(2.93*self.ctx.dps)
258 self.degree = max(16,self.dps_goal)
260 # _coeff routine requires even degree
261 if self.degree%2 > 0:
262 self.degree += 1
264 M = self.degree
266 # this is adjusting the dps of the calling context
267 # hopefully the caller doesn't monkey around with it
268 # between calling this routine and calc_time_domain_solution()
269 self.dps_orig = self.ctx.dps
270 self.ctx.dps = self.dps_goal
272 self.V = self._coeff()
273 self.p = self.ctx.matrix(self.ctx.arange(1,M+1))*self.ctx.ln2/self.t
275 # NB: p is real (mpf)
277 def _coeff(self):
278 r"""Salzer summation weights (aka, "Stehfest coefficients")
279 only depend on the approximation order (M) and the precision"""
281 M = self.degree
282 M2 = int(M/2) # checked earlier that M is even
284 V = self.ctx.matrix(M,1)
286 # Salzer summation weights
287 # get very large in magnitude and oscillate in sign,
288 # if the precision is not high enough, there will be
289 # catastrophic cancellation
290 for k in range(1,M+1):
291 z = self.ctx.matrix(min(k,M2)+1,1)
292 for j in range(int((k+1)/2),min(k,M2)+1):
293 z[j] = (self.ctx.power(j,M2)*self.ctx.fac(2*j)/
294 (self.ctx.fac(M2-j)*self.ctx.fac(j)*
295 self.ctx.fac(j-1)*self.ctx.fac(k-j)*
296 self.ctx.fac(2*j-k)))
297 V[k-1] = self.ctx.power(-1,k+M2)*self.ctx.fsum(z)
299 return V
301 def calc_time_domain_solution(self,fp,t,manual_prec=False):
302 r"""Compute time-domain Stehfest algorithm solution.
304 .. math ::
306 f(t,M) = \frac{\log 2}{t} \sum_{k=1}^{M} V_k \bar{f}\left(
307 p_k \right)
309 where
311 .. math ::
313 V_k = (-1)^{k + N/2} \sum^{\min(k,N/2)}_{i=\lfloor(k+1)/2 \rfloor}
314 \frac{i^{\frac{N}{2}}(2i)!}{\left(\frac{N}{2}-i \right)! \, i! \,
315 \left(i-1 \right)! \, \left(k-i\right)! \, \left(2i-k \right)!}
317 As the degree increases, the abscissa (`p_k`) only increase
318 linearly towards `\infty`, but the Stehfest coefficients
319 (`V_k`) alternate in sign and increase rapidly in sign,
320 requiring high precision to prevent overflow or loss of
321 significance when evaluating the sum.
323 **References**
325 1. Widder, D. (1941). *The Laplace Transform*. Princeton.
326 2. Stehfest, H. (1970). Algorithm 368: numerical inversion of
327 Laplace transforms. *Communications of the ACM* 13(1):47-49,
328 http://dx.doi.org/10.1145/361953.361969
330 """
332 # required
333 self.t = self.ctx.convert(t)
335 # assume fp was computed from p matrix returned from
336 # calc_laplace_parameter(), so is already
337 # a list or matrix of mpmath 'mpf' types
339 result = self.ctx.fdot(self.V,fp)*self.ctx.ln2/self.t
341 # setting dps back to value when calc_laplace_parameter was called
342 if not manual_prec:
343 self.ctx.dps = self.dps_orig
345 # ignore any small imaginary part
346 return result.real
348# ****************************************
350class deHoog(InverseLaplaceTransform):
352 def calc_laplace_parameter(self,t,**kwargs):
353 r"""the de Hoog, Knight & Stokes algorithm is an
354 accelerated form of the Fourier series numerical
355 inverse Laplace transform algorithms.
357 .. math ::
359 p_k = \gamma + \frac{jk}{T} \qquad 0 \le k < 2M+1
361 where
363 .. math ::
365 \gamma = \alpha - \frac{\log \mathrm{tol}}{2T},
367 `j=\sqrt{-1}`, `T = 2t_\mathrm{max}` is a scaled time,
368 `\alpha=10^{-\mathrm{dps\_goal}}` is the real part of the
369 rightmost pole or singularity, which is chosen based on the
370 desired accuracy (assuming the rightmost singularity is 0),
371 and `\mathrm{tol}=10\alpha` is the desired tolerance, which is
372 chosen in relation to `\alpha`.`
374 When increasing the degree, the abscissa increase towards
375 `j\infty`, but more slowly than the fixed Talbot
376 algorithm. The de Hoog et al. algorithm typically does better
377 with oscillatory functions of time, and less well-behaved
378 functions. The method tends to be slower than the Talbot and
379 Stehfest algorithsm, especially so at very high precision
380 (e.g., `>500` digits precision).
382 """
384 # required
385 # ------------------------------
386 self.t = self.ctx.convert(t)
388 # optional
389 # ------------------------------
390 self.tmax = kwargs.get('tmax',self.t)
392 # empirical relationships used here based on a linear fit of
393 # requested and delivered dps for exponentially decaying time
394 # functions for requested dps up to 512.
396 if 'degree' in kwargs:
397 self.degree = kwargs['degree']
398 self.dps_goal = int(1.38*self.degree)
399 else:
400 self.dps_goal = int(self.ctx.dps*1.36)
401 self.degree = max(10,self.dps_goal)
403 # 2*M+1 terms in approximation
404 M = self.degree
406 # adjust alpha component of abscissa of convergence for higher
407 # precision
408 tmp = self.ctx.power(10.0,-self.dps_goal)
409 self.alpha = self.ctx.convert(kwargs.get('alpha',tmp))
411 # desired tolerance (here simply related to alpha)
412 self.tol = self.ctx.convert(kwargs.get('tol',self.alpha*10.0))
413 self.np = 2*self.degree+1 # number of terms in approximation
415 # this is adjusting the dps of the calling context
416 # hopefully the caller doesn't monkey around with it
417 # between calling this routine and calc_time_domain_solution()
418 self.dps_orig = self.ctx.dps
419 self.ctx.dps = self.dps_goal
421 # scaling factor (likely tun-able, but 2 is typical)
422 self.scale = kwargs.get('scale',2)
423 self.T = self.ctx.convert(kwargs.get('T',self.scale*self.tmax))
425 self.p = self.ctx.matrix(2*M+1,1)
426 self.gamma = self.alpha - self.ctx.log(self.tol)/(self.scale*self.T)
427 self.p = (self.gamma + self.ctx.pi*
428 self.ctx.matrix(self.ctx.arange(self.np))/self.T*1j)
430 # NB: p is complex (mpc)
432 def calc_time_domain_solution(self,fp,t,manual_prec=False):
433 r"""Calculate time-domain solution for
434 de Hoog, Knight & Stokes algorithm.
436 The un-accelerated Fourier series approach is:
438 .. math ::
440 f(t,2M+1) = \frac{e^{\gamma t}}{T} \sum_{k=0}^{2M}{}^{'}
441 \Re\left[\bar{f}\left( p_k \right)
442 e^{i\pi t/T} \right],
444 where the prime on the summation indicates the first term is halved.
446 This simplistic approach requires so many function evaluations
447 that it is not practical. Non-linear acceleration is
448 accomplished via Pade-approximation and an analytic expression
449 for the remainder of the continued fraction. See the original
450 paper (reference 2 below) a detailed description of the
451 numerical approach.
453 **References**
455 1. Davies, B. (2005). *Integral Transforms and their
456 Applications*, Third Edition. Springer.
457 2. de Hoog, F., J. Knight, A. Stokes (1982). An improved
458 method for numerical inversion of Laplace transforms. *SIAM
459 Journal of Scientific and Statistical Computing* 3:357-366,
460 http://dx.doi.org/10.1137/0903022
462 """
464 M = self.degree
465 np = self.np
466 T = self.T
468 self.t = self.ctx.convert(t)
470 # would it be useful to try re-using
471 # space between e&q and A&B?
472 e = self.ctx.zeros(np,M+1)
473 q = self.ctx.matrix(2*M,M)
474 d = self.ctx.matrix(np,1)
475 A = self.ctx.zeros(np+1,1)
476 B = self.ctx.ones(np+1,1)
478 # initialize Q-D table
479 e[:,0] = 0.0 + 0j
480 q[0,0] = fp[1]/(fp[0]/2)
481 for i in range(1,2*M):
482 q[i,0] = fp[i+1]/fp[i]
484 # rhombus rule for filling triangular Q-D table (e & q)
485 for r in range(1,M+1):
486 # start with e, column 1, 0:2*M-2
487 mr = 2*(M-r) + 1
488 e[0:mr,r] = q[1:mr+1,r-1] - q[0:mr,r-1] + e[1:mr+1,r-1]
489 if not r == M:
490 rq = r+1
491 mr = 2*(M-rq)+1 + 2
492 for i in range(mr):
493 q[i,rq-1] = q[i+1,rq-2]*e[i+1,rq-1]/e[i,rq-1]
495 # build up continued fraction coefficients (d)
496 d[0] = fp[0]/2
497 for r in range(1,M+1):
498 d[2*r-1] = -q[0,r-1] # even terms
499 d[2*r] = -e[0,r] # odd terms
501 # seed A and B for recurrence
502 A[0] = 0.0 + 0.0j
503 A[1] = d[0]
504 B[0:2] = 1.0 + 0.0j
506 # base of the power series
507 z = self.ctx.expjpi(self.t/T) # i*pi is already in fcn
509 # coefficients of Pade approximation (A & B)
510 # using recurrence for all but last term
511 for i in range(1,2*M):
512 A[i+1] = A[i] + d[i]*A[i-1]*z
513 B[i+1] = B[i] + d[i]*B[i-1]*z
515 # "improved remainder" to continued fraction
516 brem = (1 + (d[2*M-1] - d[2*M])*z)/2
517 # powm1(x,y) computes x^y - 1 more accurately near zero
518 rem = brem*self.ctx.powm1(1 + d[2*M]*z/brem,
519 self.ctx.fraction(1,2))
521 # last term of recurrence using new remainder
522 A[np] = A[2*M] + rem*A[2*M-1]
523 B[np] = B[2*M] + rem*B[2*M-1]
525 # diagonal Pade approximation
526 # F=A/B represents accelerated trapezoid rule
527 result = self.ctx.exp(self.gamma*self.t)/T*(A[np]/B[np]).real
529 # setting dps back to value when calc_laplace_parameter was called
530 if not manual_prec:
531 self.ctx.dps = self.dps_orig
533 return result
535# ****************************************
537class LaplaceTransformInversionMethods(object):
538 def __init__(ctx, *args, **kwargs):
539 ctx._fixed_talbot = FixedTalbot(ctx)
540 ctx._stehfest = Stehfest(ctx)
541 ctx._de_hoog = deHoog(ctx)
543 def invertlaplace(ctx, f, t, **kwargs):
544 r"""Computes the numerical inverse Laplace transform for a
545 Laplace-space function at a given time. The function being
546 evaluated is assumed to be a real-valued function of time.
548 The user must supply a Laplace-space function `\bar{f}(p)`,
549 and a desired time at which to estimate the time-domain
550 solution `f(t)`.
552 A few basic examples of Laplace-space functions with known
553 inverses (see references [1,2]) :
555 .. math ::
557 \mathcal{L}\left\lbrace f(t) \right\rbrace=\bar{f}(p)
559 .. math ::
561 \mathcal{L}^{-1}\left\lbrace \bar{f}(p) \right\rbrace = f(t)
563 .. math ::
565 \bar{f}(p) = \frac{1}{(p+1)^2}
567 .. math ::
569 f(t) = t e^{-t}
571 >>> from mpmath import *
572 >>> mp.dps = 15; mp.pretty = True
573 >>> tt = [0.001, 0.01, 0.1, 1, 10]
574 >>> fp = lambda p: 1/(p+1)**2
575 >>> ft = lambda t: t*exp(-t)
576 >>> ft(tt[0]),ft(tt[0])-invertlaplace(fp,tt[0],method='talbot')
577 (0.000999000499833375, 8.57923043561212e-20)
578 >>> ft(tt[1]),ft(tt[1])-invertlaplace(fp,tt[1],method='talbot')
579 (0.00990049833749168, 3.27007646698047e-19)
580 >>> ft(tt[2]),ft(tt[2])-invertlaplace(fp,tt[2],method='talbot')
581 (0.090483741803596, -1.75215800052168e-18)
582 >>> ft(tt[3]),ft(tt[3])-invertlaplace(fp,tt[3],method='talbot')
583 (0.367879441171442, 1.2428864009344e-17)
584 >>> ft(tt[4]),ft(tt[4])-invertlaplace(fp,tt[4],method='talbot')
585 (0.000453999297624849, 4.04513489306658e-20)
587 The methods also work for higher precision:
589 >>> mp.dps = 100; mp.pretty = True
590 >>> nstr(ft(tt[0]),15),nstr(ft(tt[0])-invertlaplace(fp,tt[0],method='talbot'),15)
591 ('0.000999000499833375', '-4.96868310693356e-105')
592 >>> nstr(ft(tt[1]),15),nstr(ft(tt[1])-invertlaplace(fp,tt[1],method='talbot'),15)
593 ('0.00990049833749168', '1.23032291513122e-104')
595 .. math ::
597 \bar{f}(p) = \frac{1}{p^2+1}
599 .. math ::
601 f(t) = \mathrm{J}_0(t)
603 >>> mp.dps = 15; mp.pretty = True
604 >>> fp = lambda p: 1/sqrt(p*p + 1)
605 >>> ft = lambda t: besselj(0,t)
606 >>> ft(tt[0]),ft(tt[0])-invertlaplace(fp,tt[0])
607 (0.999999750000016, -6.09717765032273e-18)
608 >>> ft(tt[1]),ft(tt[1])-invertlaplace(fp,tt[1])
609 (0.99997500015625, -5.61756281076169e-17)
611 .. math ::
613 \bar{f}(p) = \frac{\log p}{p}
615 .. math ::
617 f(t) = -\gamma -\log t
619 >>> mp.dps = 15; mp.pretty = True
620 >>> fp = lambda p: log(p)/p
621 >>> ft = lambda t: -euler-log(t)
622 >>> ft(tt[0]),ft(tt[0])-invertlaplace(fp,tt[0],method='stehfest')
623 (6.3305396140806, -1.92126634837863e-16)
624 >>> ft(tt[1]),ft(tt[1])-invertlaplace(fp,tt[1],method='stehfest')
625 (4.02795452108656, -4.81486093200704e-16)
627 **Options**
629 :func:`~mpmath.invertlaplace` recognizes the following optional keywords
630 valid for all methods:
632 *method*
633 Chooses numerical inverse Laplace transform algorithm
634 (described below).
635 *degree*
636 Number of terms used in the approximation
638 **Algorithms**
640 Mpmath implements three numerical inverse Laplace transform
641 algorithms, attributed to: Talbot, Stehfest, and de Hoog,
642 Knight and Stokes. These can be selected by using
643 *method='talbot'*, *method='stehfest'*, or *method='dehoog'*
644 or by passing the classes *method=FixedTalbot*,
645 *method=Stehfest*, or *method=deHoog*. The functions
646 :func:`~mpmath.invlaptalbot`, :func:`~mpmath.invlapstehfest`,
647 and :func:`~mpmath.invlapdehoog` are also available as
648 shortcuts.
650 All three algorithms implement a heuristic balance between the
651 requested precision and the precision used internally for the
652 calculations. This has been tuned for a typical exponentially
653 decaying function and precision up to few hundred decimal
654 digits.
656 The Laplace transform converts the variable time (i.e., along
657 a line) into a parameter given by the right half of the
658 complex `p`-plane. Singularities, poles, and branch cuts in
659 the complex `p`-plane contain all the information regarding
660 the time behavior of the corresponding function. Any numerical
661 method must therefore sample `p`-plane "close enough" to the
662 singularities to accurately characterize them, while not
663 getting too close to have catastrophic cancellation, overflow,
664 or underflow issues. Most significantly, if one or more of the
665 singularities in the `p`-plane is not on the left side of the
666 Bromwich contour, its effects will be left out of the computed
667 solution, and the answer will be completely wrong.
669 *Talbot*
671 The fixed Talbot method is high accuracy and fast, but the
672 method can catastrophically fail for certain classes of time-domain
673 behavior, including a Heaviside step function for positive
674 time (e.g., `H(t-2)`), or some oscillatory behaviors. The
675 Talbot method usually has adjustable parameters, but the
676 "fixed" variety implemented here does not. This method
677 deforms the Bromwich integral contour in the shape of a
678 parabola towards `-\infty`, which leads to problems
679 when the solution has a decaying exponential in it (e.g., a
680 Heaviside step function is equivalent to multiplying by a
681 decaying exponential in Laplace space).
683 *Stehfest*
685 The Stehfest algorithm only uses abscissa along the real axis
686 of the complex `p`-plane to estimate the time-domain
687 function. Oscillatory time-domain functions have poles away
688 from the real axis, so this method does not work well with
689 oscillatory functions, especially high-frequency ones. This
690 method also depends on summation of terms in a series that
691 grows very large, and will have catastrophic cancellation
692 during summation if the working precision is too low.
694 *de Hoog et al.*
696 The de Hoog, Knight, and Stokes method is essentially a
697 Fourier-series quadrature-type approximation to the Bromwich
698 contour integral, with non-linear series acceleration and an
699 analytical expression for the remainder term. This method is
700 typically the most robust and is therefore the default
701 method. This method also involves the greatest amount of
702 overhead, so it is typically the slowest of the three methods
703 at high precision.
705 **Singularities**
707 All numerical inverse Laplace transform methods have problems
708 at large time when the Laplace-space function has poles,
709 singularities, or branch cuts to the right of the origin in
710 the complex plane. For simple poles in `\bar{f}(p)` at the
711 `p`-plane origin, the time function is constant in time (e.g.,
712 `\mathcal{L}\left\lbrace 1 \right\rbrace=1/p` has a pole at
713 `p=0`). A pole in `\bar{f}(p)` to the left of the origin is a
714 decreasing function of time (e.g., `\mathcal{L}\left\lbrace
715 e^{-t/2} \right\rbrace=1/(p+1/2)` has a pole at `p=-1/2`), and
716 a pole to the right of the origin leads to an increasing
717 function in time (e.g., `\mathcal{L}\left\lbrace t e^{t/4}
718 \right\rbrace = 1/(p-1/4)^2` has a pole at `p=1/4`). When
719 singularities occur off the real `p` axis, the time-domain
720 function is oscillatory. For example `\mathcal{L}\left\lbrace
721 \mathrm{J}_0(t) \right\rbrace=1/\sqrt{p^2+1}` has a branch cut
722 starting at `p=j=\sqrt{-1}` and is a decaying oscillatory
723 function, This range of behaviors is illustrated in Duffy [3]
724 Figure 4.10.4, p. 228.
726 In general as `p \rightarrow \infty` `t \rightarrow 0` and
727 vice-versa. All numerical inverse Laplace transform methods
728 require their abscissa to shift closer to the origin for
729 larger times. If the abscissa shift left of the rightmost
730 singularity in the Laplace domain, the answer will be
731 completely wrong (the effect of singularities to the right of
732 the Bromwich contour are not included in the results).
734 For example, the following exponentially growing function has
735 a pole at `p=3`:
737 .. math ::
739 \bar{f}(p)=\frac{1}{p^2-9}
741 .. math ::
743 f(t)=\frac{1}{3}\sinh 3t
745 >>> mp.dps = 15; mp.pretty = True
746 >>> fp = lambda p: 1/(p*p-9)
747 >>> ft = lambda t: sinh(3*t)/3
748 >>> tt = [0.01,0.1,1.0,10.0]
749 >>> ft(tt[0]),invertlaplace(fp,tt[0],method='talbot')
750 (0.0100015000675014, 0.0100015000675014)
751 >>> ft(tt[1]),invertlaplace(fp,tt[1],method='talbot')
752 (0.101506764482381, 0.101506764482381)
753 >>> ft(tt[2]),invertlaplace(fp,tt[2],method='talbot')
754 (3.33929164246997, 3.33929164246997)
755 >>> ft(tt[3]),invertlaplace(fp,tt[3],method='talbot')
756 (1781079096920.74, -1.61331069624091e-14)
758 **References**
760 1. [DLMF]_ section 1.14 (http://dlmf.nist.gov/1.14T4)
761 2. Cohen, A.M. (2007). Numerical Methods for Laplace Transform
762 Inversion, Springer.
763 3. Duffy, D.G. (1998). Advanced Engineering Mathematics, CRC Press.
765 **Numerical Inverse Laplace Transform Reviews**
767 1. Bellman, R., R.E. Kalaba, J.A. Lockett (1966). *Numerical
768 inversion of the Laplace transform: Applications to Biology,
769 Economics, Engineering, and Physics*. Elsevier.
770 2. Davies, B., B. Martin (1979). Numerical inversion of the
771 Laplace transform: a survey and comparison of methods. *Journal
772 of Computational Physics* 33:1-32,
773 http://dx.doi.org/10.1016/0021-9991(79)90025-1
774 3. Duffy, D.G. (1993). On the numerical inversion of Laplace
775 transforms: Comparison of three new methods on characteristic
776 problems from applications. *ACM Transactions on Mathematical
777 Software* 19(3):333-359, http://dx.doi.org/10.1145/155743.155788
778 4. Kuhlman, K.L., (2013). Review of Inverse Laplace Transform
779 Algorithms for Laplace-Space Numerical Approaches, *Numerical
780 Algorithms*, 63(2):339-355.
781 http://dx.doi.org/10.1007/s11075-012-9625-3
783 """
785 rule = kwargs.get('method','dehoog')
786 if type(rule) is str:
787 lrule = rule.lower()
788 if lrule == 'talbot':
789 rule = ctx._fixed_talbot
790 elif lrule == 'stehfest':
791 rule = ctx._stehfest
792 elif lrule == 'dehoog':
793 rule = ctx._de_hoog
794 else:
795 raise ValueError("unknown invlap algorithm: %s" % rule)
796 else:
797 rule = rule(ctx)
799 # determine the vector of Laplace-space parameter
800 # needed for the requested method and desired time
801 rule.calc_laplace_parameter(t,**kwargs)
803 # compute the Laplace-space function evalutations
804 # at the required abscissa.
805 fp = [f(p) for p in rule.p]
807 # compute the time-domain solution from the
808 # Laplace-space function evaluations
809 return rule.calc_time_domain_solution(fp,t)
811 # shortcuts for the above function for specific methods
812 def invlaptalbot(ctx, *args, **kwargs):
813 kwargs['method'] = 'talbot'
814 return ctx.invertlaplace(*args, **kwargs)
816 def invlapstehfest(ctx, *args, **kwargs):
817 kwargs['method'] = 'stehfest'
818 return ctx.invertlaplace(*args, **kwargs)
820 def invlapdehoog(ctx, *args, **kwargs):
821 kwargs['method'] = 'dehoog'
822 return ctx.invertlaplace(*args, **kwargs)
824# ****************************************
826if __name__ == '__main__':
827 import doctest
828 doctest.testmod()