api.lti

talbot_get_omega(n, alpha)

ARR = TALBOT_GET_OMEGA(N, ALPHA)

Talbot contour weights. Twin of the native Python api.lti.talbot_get_omega.

talbot_get_alpha(n)

ARR = TALBOT_GET_ALPHA(N)

Talbot contour nodes. Twin of the native Python api.lti.talbot_get_alpha.

laplace_weeks_scaling(F, p0, tol)

[SIGMA, B, Q] = LAPLACE_WEEKS_SCALING(F, P0, TOL)

Automatic selection of the exponential damping SIGMA and the scaling B for the Laguerre inversion, following the algorithm of Fig. 1 of P. G. Harrison and W. J. Knottenbelt, “Passage Time Distributions in Large Markov Chains”, 2002. The search accepts the first (sigma, b) at which the coefficients have decayed by term P0:

sigma = 0; b = 1 while |q_p0| > tol or |q_{p0+1}| > tol

sigma = 0.001 if sigma == 0 else 2*sigma if sigma > 0.2

b = b + 4; if b > 10, no suitable parameters exist sigma = 0

Raising b too far is numerically counterproductive and excessive damping is unstable in finite precision, which is why the search is bounded rather than unbounded. When it exhausts the box the failure is REFUSED BY NAME: a density with a discontinuity in itself or its derivatives has no usable Laguerre representation (Sec. 4.2), and returning the last iterate would report noise as an answer. Use ‘euler’ for those, at roughly 50 transform evaluations per time point.

laplace_weeks_coeffs(F, sigma, b, p0)

Q = LAPLACE_WEEKS_COEFFS(F, SIGMA, B, P0)

Laguerre coefficients q_n, n = 0..2*p0-1, of the damped and scaled function

f_{sigma,b}(t) = exp(-sigma t) f(t/b)

whose Laguerre generating function is (Harrison and Knottenbelt 2002, Sec. 4.2)

Q_{sigma,b}(z) = b/(1-z) * L( b(1+z)/(2(1-z)) + b*sigma )

and q_n = (1/(2 pi i)) contour-integral Q(z)/z^{n+1} dz on |z| = r.

NOTE ON THE PAPER. Eq. 10 as printed carries the factor (1-z) rather than 1/(1-z). The scaled form quoted above, printed later in the same section, carries 1/(1-z) and is the correct one: with l_n(t) = exp(-t/2) L_n(t) the transform of l_n is (s-1/2)^n/(s+1/2)^{n+1}, so L(s) = Q(z)/(s+1/2) with z = (s-1/2)/(s+1/2), and s+1/2 = 1/(1-z). Implementing the printed (1-z) gives a wrong answer at every t (163 per cent at t = 0.1 on Exp(2)), so the discrepancy is a typo and not a convention.

Sec. 4.3 fixes the number of trapezoids at 2*p0 and the radius at r = 0.1^(4/p0) for every n, rather than letting them grow with n. The resulting quadrature is a discrete Fourier transform of Q sampled on the circle, so all 2*p0 coefficients come out of one FFT and the transform is evaluated 2*p0 times in total.

laplace_invert_weeks(F, t, q, sigma, b)

F = LAPLACE_INVERT_WEEKS(F, T, Q, SIGMA, B)

Invert a Laplace transform by the Laguerre (Weeks) series

f(t) = sum_{n>=0} q_n l_n(t)

of W. Weeks, “Numerical inversion of Laplace transforms using Laguerre functions”, J. ACM 13, 1966, in the form used by P. G. Harrison and W. J. Knottenbelt, “Passage Time Distributions in Large Markov Chains”, 2002, Sec. 4.1-4.3, after J. Abate, G. Choudhury and W. Whitt, “On the Laguerre method for numerically inverting Laplace transforms”, INFORMS J. Computing 8(4), 1996.

F transform handle, called with a COMPLEX argument T time points (vector); t <= 0 returns 0 Q Laguerre coefficients from LAPLACE_WEEKS_COEFFS; omit or leave empty

to have LAPLACE_WEEKS_SCALING pick sigma and b and compute them

SIGMA, B the exponential damping and scaling parameters that Q was built: with; required when Q is supplied

Unlike Euler and Talbot the coefficients do not depend on t, so one coefficient set serves an arbitrary number of time points. That is the property this method is here for: the transform is evaluated 2*p0 times in total, not 2*p0 times per t.

laplace_invert_talbot(F, t, n)

RES = LAPLACE_INVERT_TALBOT(F, T, N)

Invert a Laplace transform at T by Talbot’s deformed contour. F is a handle called with a COMPLEX argument. N defaults to 32. Twin of the native Python api.lti.laplace_invert_talbot.

laplace_invert_pdf(F, tset, method, n)

RES = LAPLACE_INVERT_PDF(F, TSET, METHOD, N)

Invert the transform of a DENSITY on a grid. Values at t <= 0 are 0 and the result is clamped at zero. Twin of the native Python api.lti.laplace_invert_pdf.

laplace_invert_gaver_stehfest(F, t, n)

RES = LAPLACE_INVERT_GAVER_STEHFEST(F, T, N)

Invert a Laplace transform at T by Gaver-Stehfest. This method samples the REAL axis only, so F may be a real-argument handle. N defaults to 12 and is rounded DOWN to even. Twin of the native Python api.lti.laplace_invert_gaver_stehfest.

laplace_invert_euler(F, t, n)

RES = LAPLACE_INVERT_EULER(F, T, N)

Invert a Laplace transform at T by the Euler (Abate-Whitt) method. F is a handle called with a COMPLEX argument.

N defaults to 41 and is rounded UP to odd. THE LARGER N IS NOT THE BETTER ONE: the weights carry a factor 10^((n-1)/6) against an ALTERNATING sum, so accuracy is a race between the series converging and the cancellation eating the mantissa. Worst relative error on F(s) = 2/(s+2) over t in {0.1, 0.5, 1, 2}:

n = 11 21 31 41 51 71 99 err = 4.4e-3 2.1e-6 1.6e-9 1.6e-10 4.7e-8 1.7e-4 1.4e+0

At 99 the result is 140 per cent wrong and negative at some t. In double precision 41 is the optimum; the native Python twin carries the same table.

laplace_invert_cme(F, t, maxFnEvals)

RES = LAPLACE_INVERT_CME(F, T, MAXFNEVALS)

Invert a Laplace transform at T by the Concentrated Matrix Exponential method of Horvath, Horvath and Telek, using the pre-computed coefficient tables shipped with LINE. Delegates to the vendored MATLAB_ILT, which owns the tables; this wrapper exists so that api/lti presents the same surface in all four codebases.

laplace_invert_cdf(F, tset, method, n)

RES = LAPLACE_INVERT_CDF(F, TSET, METHOD, N)

Invert the transform of a CDF on a grid. F IS THE TRANSFORM OF THE DENSITY, not of the CDF: this function forms F(s)/s itself. Values at t <= 0 are 0; the result is clipped to [0,1] and made monotone, since neither property is guaranteed by a truncated inversion. Twin of the native Python api.lti.laplace_invert_cdf.

laplace_invert(F, t, method, n)

RES = LAPLACE_INVERT(F, T, METHOD, N)

Invert a Laplace transform at a single time point T. METHOD is one of ‘euler’ (default), ‘talbot’, ‘gaver-stehfest’, ‘cme’ or ‘weeks’. N is the number of terms and defaults per method: 41 euler, 32 talbot, 12 gaver-stehfest, 25 cme, 200 weeks (the p0 of the coefficient search).

‘weeks’ recomputes its coefficients on every call, which wastes the one property that method has. Call LAPLACE_WEEKS_SCALING once and pass the coefficients to LAPLACE_INVERT_WEEKS when inverting on a grid.

Twin of the native Python api.lti.laplace_invert.

gaver_stehfest_get_omega(n)

RES = GAVER_STEHFEST_GET_OMEGA(N)

Gaver-Stehfest weights. N is rounded DOWN to even. Twin of the native Python api.lti.gaver_stehfest_get_omega.

gaver_stehfest_get_alpha(n)

RES = GAVER_STEHFEST_GET_ALPHA(N)

Gaver-Stehfest nodes alpha_k = k log 2, k = 1..N. N is rounded DOWN to even (the method is defined for even N only); Euler rounds UP instead, and the two conventions must not be swapped. Twin of the native Python api.lti.gaver_stehfest_get_alpha.

euler_get_omega(n)

RES = EULER_GET_OMEGA(N)

Euler (Abate-Whitt) weights omega_i = 10^((n-1)/6) (-1)^(i-1) eta_i. Twin of the native Python api.lti.euler_get_omega.

euler_get_eta(n)

RES = EULER_GET_ETA(N)

Euler (Abate-Whitt) weights eta_i, i = 1..N, N odd. Twin of the native Python api.lti.euler_get_eta.

euler_get_alpha(n)

RES = EULER_GET_ALPHA(N)

Euler (Abate-Whitt) nodes alpha_i = (n-1)log(10)/6 + i*pi*1i, i = 0..N-1. Twin of the native Python api.lti.euler_get_alpha.