Package jline.api.lti

Class Laplace_invert

java.lang.Object

jline.api.lti.Laplace_invert

public final class Laplace_invert
extends Object

Numerical inversion of a Laplace transform.

Euler, Talbot, Gaver-Stehfest and CME are delegated to iltcme, whose Abate-Whitt evaluator and coefficient tables are already in the tree. WEEKS (the Laguerre series) is implemented here and had no counterpart in any codebase: jline.lib.lti.laguerre is Gauss-Laguerre QUADRATURE, which is a different thing, and the native Python lib_lti_laguerre silently delegated to Talbot.

Reference for Weeks: W. Weeks, J. ACM 13, 1966; J. Abate, G. Choudhury and W. Whitt, INFORMS J. Computing 8(4), 1996; P. G. Harrison and W. J. Knottenbelt, "Passage Time Distributions in Large Markov Chains", 2002, Sec. 4.1-4.3, whose Fig. 1 is the automatic scaling search below.

Field Summary

Fields

Modifier and Type	Field	Description
static final int	WEEKS_DEFAULT_P0	Default number of trapezoid pairs in the Laguerre quadrature.

Method Summary

Modifier and Type	Method	Description
static double	laplace_invert(UnaryOperator<org.apache.commons.math3.complex.Complex> F, double t, String method, int n)	Invert F at t by the named method: "euler", "talbot", "gaver-stehfest", "cme" or "weeks".
static double[]	laplace_invert_cdf(UnaryOperator<org.apache.commons.math3.complex.Complex> F, double[] t, String method, int n)	The DISTRIBUTION on a grid, from the transform of the DENSITY.
static double	laplace_invert_cme(UnaryOperator<org.apache.commons.math3.complex.Complex> F, double t, int n)	Invert F at t by the Concentrated Matrix Exponential method.
static double	laplace_invert_euler(UnaryOperator<org.apache.commons.math3.complex.Complex> F, double t, int n)	Invert F at t by the Euler (Abate-Whitt) method.
static double	laplace_invert_gaver_stehfest(UnaryOperator<org.apache.commons.math3.complex.Complex> F, double t, int n)	Invert F at t by Gaver-Stehfest.
static double[]	laplace_invert_pdf(UnaryOperator<org.apache.commons.math3.complex.Complex> F, double[] t, String method, int n)	The DENSITY on a grid: the inversion clamped at zero.
static double	laplace_invert_talbot(UnaryOperator<org.apache.commons.math3.complex.Complex> F, double t, int n)	Invert F at t by Talbot's deformed contour.
static double	laplace_invert_weeks(UnaryOperator<org.apache.commons.math3.complex.Complex> F, double t)	Convenience: build the parameters, then invert at one point.
static double	laplace_invert_weeks(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).
static double[]	laplace_weeks_coeffs(UnaryOperator<org.apache.commons.math3.complex.Complex> F, double sigma, double b, int 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 generating function is
static WeeksParams	laplace_weeks_scaling(UnaryOperator<org.apache.commons.math3.complex.Complex> F)	The search at the default p0 = 200 and tolerance 1e-10.
static WeeksParams	laplace_weeks_scaling(UnaryOperator<org.apache.commons.math3.complex.Complex> F, int p0, double tol)	The automatic (sigma, b) search of Fig.
static org.apache.commons.math3.complex.Complex[]	talbot_get_alpha(int n)	Talbot contour nodes alpha_i, i = 1..n.
static org.apache.commons.math3.complex.Complex[]	talbot_get_omega(int n, org.apache.commons.math3.complex.Complex[] alpha)	Talbot contour weights omega_i, i = 1..n, for the nodes of ALPHA.

Methods inherited from class java.lang.Object

clone, equals, finalize, getClass, hashCode, notify, notifyAll, toString, wait, wait, wait

Field Details

WEEKS_DEFAULT_P0

public static final int WEEKS_DEFAULT_P0

Default number of trapezoid pairs in the Laguerre quadrature.

See Also:

• Constant Field Values

Method Details

laplace_weeks_coeffs

public static double[] laplace_weeks_coeffs(UnaryOperator<org.apache.commons.math3.complex.Complex> F,
 double sigma,
 double b,
 int 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 generating function is

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

NOTE ON THE PAPER. Eq. 10 as printed carries the factor (1-z) rather than 1/(1-z). The scaled form 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) is wrong at every t (163 per cent at t = 0.1 on Exp(2)).

Sec. 4.3 fixes the trapezoid count at 2*p0 and the radius at r = 0.1^(4/p0) for every n, so the quadrature is one discrete Fourier transform of Q sampled on the circle and the transform is evaluated 2*p0 times IN TOTAL rather than per coefficient. The DFT is evaluated directly: at 2*p0 = 400 points that is 160k complex multiplies, which is not worth a dependency.

laplace_weeks_scaling

public static WeeksParams laplace_weeks_scaling(UnaryOperator<org.apache.commons.math3.complex.Complex> F,
 int p0,
 double tol)

The automatic (sigma, b) search of Fig. 1: accept the first pair at which the coefficients have decayed by term p0, doubling sigma from 0.001 and stepping b by 4 whenever sigma passes 0.2.

REFUSES BY NAME when the box is exhausted. Raising b further is counterproductive and excessive damping is unstable in finite precision, and a density with a discontinuity in itself or its derivatives has no usable Laguerre representation at all (Sec. 4.2). Returning the last iterate would report noise as an answer; Euler handles those cases.

laplace_weeks_scaling

public static WeeksParams laplace_weeks_scaling(UnaryOperator<org.apache.commons.math3.complex.Complex> F)

The search at the default p0 = 200 and tolerance 1e-10.

laplace_invert_weeks

public static double laplace_invert_weeks(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).

Unlike Euler and Talbot the coefficients do not depend on t, so ONE parameter set serves an arbitrary number of time points: the transform is evaluated 2*p0 times in total, not 2*p0 times per t. That is the property this method is here for, so build the WeeksParams once and reuse it on a grid.

laplace_invert_weeks

public static double laplace_invert_weeks(UnaryOperator<org.apache.commons.math3.complex.Complex> F,
 double t)

Convenience: build the parameters, then invert at one point.

talbot_get_alpha

public static org.apache.commons.math3.complex.Complex[] talbot_get_alpha(int n)

Talbot contour nodes alpha_i, i = 1..n.

The nodes and weights themselves come from talbot, which already carried them; only the inverter below was missing, so the class was unreachable. The Euler and Gaver-Stehfest nodes have no counterpart on this surface on purpose: they are built inline by iltcme.abateWhittWeights, the one implementation the JAR carries for those two. Talbot deforms the contour rather than shifting the real axis, so it is not an Abate-Whitt scheme and that evaluator cannot host it.

talbot_get_omega

public static org.apache.commons.math3.complex.Complex[] talbot_get_omega(int n,
 org.apache.commons.math3.complex.Complex[] alpha)

Talbot contour weights omega_i, i = 1..n, for the nodes of ALPHA.

laplace_invert_talbot

public static double laplace_invert_talbot(UnaryOperator<org.apache.commons.math3.complex.Complex> F,
 double t,
 int n)

Invert F at t by Talbot's deformed contour. N defaults to 32.

Talbot samples F off the real axis, so F must accept a genuinely complex argument; unlike Gaver-Stehfest it cannot be fed a real-only transform.

laplace_invert_euler

public static double laplace_invert_euler(UnaryOperator<org.apache.commons.math3.complex.Complex> F,
 double t,
 int n)

Invert F at t by the Euler (Abate-Whitt) method. N defaults to 41.

laplace_invert_gaver_stehfest

public static double laplace_invert_gaver_stehfest(UnaryOperator<org.apache.commons.math3.complex.Complex> F,
 double t,
 int n)

Invert F at t by Gaver-Stehfest. N defaults to 12.

laplace_invert_cme

public static double laplace_invert_cme(UnaryOperator<org.apache.commons.math3.complex.Complex> F,
 double t,
 int n)

Invert F at t by the Concentrated Matrix Exponential method. N defaults to 25.

laplace_invert

public static double laplace_invert(UnaryOperator<org.apache.commons.math3.complex.Complex> F,
 double t,
 String method,
 int n)

Invert F at t by the named method: "euler", "talbot", "gaver-stehfest", "cme" or "weeks".

Parameters:

n - number of terms; 0 takes the method's own default (41 Euler, 32 Talbot, 12 Gaver-Stehfest, 25 CME, 200 the Weeks p0)

laplace_invert_pdf

public static double[] laplace_invert_pdf(UnaryOperator<org.apache.commons.math3.complex.Complex> F,
 double[] t,
 String method,
 int n)

The DENSITY on a grid: the inversion clamped at zero.

A density cannot be negative, and a numerical inversion can undershoot near the origin or in a tail.

laplace_invert_cdf

public static double[] laplace_invert_cdf(UnaryOperator<org.apache.commons.math3.complex.Complex> F,
 double[] t,
 String method,
 int n)

The DISTRIBUTION on a grid, from the transform of the DENSITY.

F(s)/s is the transform of the CDF, so that is what is inverted -- passing the CDF's own transform here would invert it twice. The result is clamped into [0,1] and made monotone by a running maximum, because a numerical inversion is pointwise and nothing in it enforces either property; a non-monotone "CDF" then yields negative probabilities downstream.