LINE Solver (C++)
Templated C++ port of the LINE queueing solver
Loading...
Searching...
No Matches
laplace_invert.h File Reference
Algorithm layer » Transform inversion (lti)

Numerical inversion of a Laplace transform: Euler, Talbot, Gaver-Stehfest. More...

#include <algorithm>
#include <cmath>
#include <complex>
#include <cstddef>
#include <functional>
#include <string>
#include <vector>
#include "line/api/mam/matlab_ilt.h"
#include "line/util/error.h"
Include dependency graph for laplace_invert.h:

Go to the source code of this file.

Classes

struct  line::lti::WeeksParams
 A Laguerre expansion: the damping, the scaling and the coefficients. More...

Namespaces

namespace  line
namespace  line::lti

Typedefs

using line::lti::Cplx = std::complex<double>
using line::lti::LaplaceFn = std::function<Cplx(Cplx)>
 The transform, evaluated at complex argument.
using line::lti::RealLaplaceFn = std::function<double(double)>
 A transform that can only be evaluated on the real axis.

Enumerations

enum class  line::lti::LaplaceMethod {
  line::lti::Euler = 0 , line::lti::Talbot , line::lti::GaverStehfest , line::lti::Cme ,
  line::lti::Weeks
}
 The methods laplace_invert accepts. More...

Functions

std::vector< Cplx > line::lti::euler_get_alpha (std::size_t n)
 Euler nodes: a vertical line at Re = (n-1) log(10) / 6.
std::vector< double > line::lti::euler_get_eta (std::size_t n)
 Euler weights before the alternating sign and the scale.
std::vector< Cplx > line::lti::euler_get_omega (std::size_t n)
 Euler weights: eta, alternating in sign, scaled by 10^((n-1)/6).
std::vector< Cplx > line::lti::talbot_get_alpha (std::size_t n)
 Talbot nodes: the cotangent contour, bending into the left half plane.
std::vector< Cplx > line::lti::talbot_get_omega (std::size_t n, const std::vector< Cplx > &alpha)
 Talbot weights, which carry the contour's own derivative.
std::vector< double > line::lti::gaver_stehfest_get_alpha (std::size_t n)
 Gaver-Stehfest nodes: k log 2, on the REAL axis.
std::vector< double > line::lti::gaver_stehfest_get_omega (std::size_t n)
 Gaver-Stehfest weights.
double line::lti::laplace_invert_euler (const LaplaceFn &F, double t, std::size_t n=41)
 Euler inversion of F at t.
double line::lti::laplace_invert_talbot (const LaplaceFn &F, double t, std::size_t n=32)
 Talbot inversion of F at t.
double line::lti::laplace_invert_gaver_stehfest (const RealLaplaceFn &F, double t, std::size_t n=12)
 Gaver-Stehfest inversion of F at t.
std::vector< double > line::lti::laplace_weeks_coeffs (const LaplaceFn &F, double sigma=0.0, double b=1.0, std::size_t p0=200)
 Laguerre coefficients q_n, n = 0..2*p0-1, of f_{sigma,b}(t) = exp(-sigma t) f(t/b), whose generating function is.
WeeksParams line::lti::laplace_weeks_scaling (const LaplaceFn &F, std::size_t p0=200, double tol=1e-10)
 The automatic (sigma, b) search of Fig.
double line::lti::laplace_invert_weeks (const WeeksParams &w, double t)
 Invert by the Laguerre series f(t) = sum_n q_n l_n(t), recovered as exp(sigma*b*t) f_{sigma,b}(b*t).
double line::lti::laplace_invert_weeks (const LaplaceFn &F, double t, std::size_t p0=200)
 Convenience overload: build the parameters, then invert at one point.
LaplaceMethod line::lti::laplace_method (const std::string &s)
 Parse the reference's method names, including its two Gaver spellings.
double line::lti::laplace_invert (const LaplaceFn &F, double t, LaplaceMethod method=LaplaceMethod::Euler, std::size_t n=0)
 Invert F at t by the named method.
std::vector< double > line::lti::laplace_invert_pdf (const LaplaceFn &F, const std::vector< double > &t, LaplaceMethod method=LaplaceMethod::Euler, std::size_t n=0)
 The DENSITY on a grid: the inversion clamped at zero.
std::vector< double > line::lti::laplace_invert_cdf (const LaplaceFn &F, const std::vector< double > &t, LaplaceMethod method=LaplaceMethod::Euler, std::size_t n=0)
 The DISTRIBUTION on a grid, from the transform of the DENSITY.

Detailed Description

Numerical inversion of a Laplace transform: Euler, Talbot, Gaver-Stehfest.

Port of python/line_solver/api/lti/__init__.py. THIS IS PYTHON-ONLY: MATLAB carries only the CME method (matlab/lib/thirdparty/iltcme/matlab_ilt.m, already ported as api/mam/matlab_ilt.h), so native Python is the reference for the other three.

ALL FOUR ARE THE SAME FRAMEWORK. Abate-Whitt writes 

f(t) ~ (1/t) sum_k Re[ omega_k F(alpha_k / t) ],

and a method IS its (alpha, omega) pair – nothing else differs, which is why they share one evaluator here. What differs is where the nodes sit:

  • EULER puts them on a vertical line and accelerates an alternating series with binomial (Euler) weights. Odd n only; the reference silently rounds an even n up, and so does this. ITS DEFAULT IS NOT THE REFERENCE'S – see the note on laplace_invert_euler.
  • TALBOT deforms the contour into the left half plane, where the transform decays, so it needs far fewer nodes – 32 against Euler's 99. It requires F to be analytic there, which a rational transform is and a transform with a branch cut is not.
  • GAVER-STEHFEST samples F on the REAL axis only, which is what makes it the one usable method when the transform cannot be evaluated at complex argument. It pays for that in conditioning: the weights alternate in sign and grow, so it needs high precision and n even (again rounded).
  • CME is api/mam/matlab_ilt.h, whose coefficient table is vendored; it is not duplicated here, and laplace_invert dispatches to it.

The defaults are the reference's own for Talbot (32) and Gaver-Stehfest (12). EULER'S IS NOT: the old default was 99, which is unusable in double precision; both this port and native Python now default to 41. The measurement is at laplace_invert_euler.

ARITHMETIC: double. Every one of these is a floating-point quadrature.

Definition in file laplace_invert.h.