deutsch 
Numeric Laplace transform inversion algorithms
Interface definition
/*
Copyright (C) 2017 Torsten Brischalle
email: torsten@brischalle.de
web: http://www.aaabbb.de
Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation
files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy,
modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software
is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR
IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
using System;
using System.Collections.Generic;
using System.Linq;
using System.Numerics;
using System.Text;
using System.Threading.Tasks;
namespace LaplaceInversion
{
public delegate Complex LaplaceFunction(Complex s);
public interface ILaplaceInversion
{
string Name { get; }
double[] Calculate(LaplaceFunction F, double[] t);
}
}
Koizumi algorithm (FFT based)
/*
Copyright (C) 2017 Torsten Brischalle
email: torsten@brischalle.de
web: http://www.aaabbb.de
ported from:
numerical laplace transform inversion algorithm by
Prof. Dr. Helmut Weber
web: https://www.cs.hs-rm.de/~weber/
Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation
files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy,
modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software
is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR
IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
using System;
using System.Collections.Generic;
using System.Linq;
using System.Numerics;
using System.Text;
using System.Threading.Tasks;
namespace LaplaceInversion
{
/* --------------------------------------------------------------------
Function: calculates the numercal inverse laplace transform
/00
F(s) = | exp(-s*t)f(t)dt, Re(s) > 0
0/
with the Koizumi algorithm.
Input: n - number of fourier terms
Literature: 1. W. Squire, The Numerical Treatment of Laplace Transforms:
I. The Koizumi Inversion Method, International Journal for
Numerical Methods in Engineering, 20, 1697-1702 (1984)
2. S. Koizumi, A new method of evaluation of the Heaviside
operational expression by Fourier sines, Phil. Mag. 19 (1935)
---------------------------------------------------------------------
*/
public class LaplaceInversionKoizumi : ILaplaceInversion
{
private int _n;
private double[] _coeff;
public string Name { get { return "Koizumi"; } }
public LaplaceInversionKoizumi(int n = 2048)
{
if (n < 1)
throw new Exception("n < 1");
_n = n;
_coeff = new double[n];
}
public double[] Calculate(LaplaceFunction F, double[] tArr)
{
double[] yArr = new double[tArr.Length];
double T = tArr.Max() * 4;
double v1 = Math.PI / (2 * T);
double v2 = 2.0 / T;
int i;
// calculate coefficients
for (i = 0; i < _n; i++)
{
Complex z;
z = F(Complex.ImaginaryOne * (1.0 - 2.0 * (i + 1)) * v1);
_coeff[i] = z.Imaginary * v2;
}
// evaluate fourier series for each t
for (i = 0; i < tArr.Length; i++)
{
double t = tArr[i];
if (t <= 0)
yArr[i] = 0;
else
{
double arg = t * v1;
double ct = 2.0 * Math.Cos(2.0 * arg);
double c2 = 0.0;
double c1 = _coeff[_n - 1];
int k;
for (k = _n - 2; k >= 0; k--)
{
double ck = ct * c1 - c2 + _coeff[k];
c2 = c1;
c1 = ck;
}
yArr[i] = (c1 + c2) * Math.Sin(arg);
}
}
return yArr;
}
}
}
Huddleston algorithm (FFT based)
/*
Copyright (C) 2017 Torsten Brischalle
email: torsten@brischalle.de
web: http://www.aaabbb.de
ported from:
numerical inverse laplace transform algorithm by
Fausto Arinos de Almeida Barbuto (Calgary, Canada)
E-mail: fausto_barbuto@yahoo.ca
web: http://code.activestate.com/recipes/128243-numerical-inversion-of-laplace-transforms-using-th/
Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation
files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy,
modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software
is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR
IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
using System;
using System.Collections.Generic;
using System.Linq;
using System.Numerics;
using System.Text;
using System.Threading.Tasks;
namespace LaplaceInversion
{
/*
#################################################################
# Function InvLap(t,omega,sigma,nint), numerically inverts a #
# Laplace transform F(s) into f(t) using the Fast Fourier #
# Transform (FFT) algorithm for a specific time "t", an #
# upper frequency limit "omega", a real parameter "sigma" #
# and the number of integration intervals "nint" . #
# #
# Function F(s) is defined in separate as Fs(s) (see code #
# below). Fs(s) has to be changed accordingly everytime the #
# user wants to invert a different function. #
# #
# I suggest to use omega>100 and nint=50*omega. The higher #
# the values for omega, the more accurate the results will be #
# in general, but at the expense of longer processing times. #
# #
# Sigma is a real number which must be a little bigger than #
# the real part of rightmost pole of the function F(s). For #
# example, F(s) = 1/s + 1/(s-2) + 1/(s+1) has poles for s=0, #
# s=2 and s=-1. Hence, sigma must be made equal to, say, #
# 2.05 so as to keep all poles at the left of this value. #
# The analytical inverse for this simple function is #
# f(t) = 1 + exp(-t) + exp(2t). For t=1.25, omega=200, #
# nint=10000 and sigma=2.05, the numerical inversion yields #
# f(1.25) ~= 13.456844516, or -0.09% away from the actual #
# analytical result, 13.468998757 (results truncated to 9 #
# decimal places). This is the example used in this code. #
# #
# Creator: Fausto Arinos de Almeida Barbuto (Calgary, Canada) #
# Date: May 18, 2002 #
# E-mail: fausto_barbuto@yahoo.ca #
# #
# Reference: #
# Huddleston, T. and Byrne, P: "Numerical Inversion of #
# Laplace Transforms", University of South Alabama, April #
# 1999 (found at http://www.eng.usouthal.edu/huddleston/ #
# SoftwareSupport/Download/Inversion99.doc) #
# #
# Usage: invoke InvLap(t,omega,sigma,nint), for t>0. #
# #
#################################################################
*/
public class LaplaceInversionHuddlestonByrne : ILaplaceInversion
{
private double _sigma;
private double _omega;
private int _n;
public string Name { get { return "Huddleston/Byrne"; } }
//-------------------------------------------------------------------------------------------------
public LaplaceInversionHuddlestonByrne(double sigma, double omega = 200, int n = 10000)
{
this._sigma = sigma;
this._omega = omega;
this._n = n;
}
//-------------------------------------------------------------------------------------------------
public double[] Calculate(LaplaceFunction F, double[] tArr)
{
double[] yArr = new double[tArr.Length];
int idx;
for (idx = 0; idx < tArr.Length; idx++)
{
double t = tArr[idx];
if (t == 0)
{
yArr[idx] = 0;
continue;
}
double sum = 0.0;
double delta = this._omega / this._n;
double wi = 0.0;
// The for-loop below computes the FFT Inversion Algorithm.
// It is in fact the trapezoidal rule for numerical integration.
int i;
for (i = 0; i < this._n; i++)
{
Complex witi = new Complex(0, wi * t);
double wf = wi + delta;
Complex wfti = new Complex(0, wf * t);
double fi = (Complex.Exp(witi) * F(new Complex(this._sigma, wi))).Real;
double ff = (Complex.Exp(wfti) * F(new Complex(this._sigma, wf))).Real;
sum = sum + 0.5 * (wf - wi) * (fi + ff);
wi = wf;
}
yArr[idx] = sum * Math.Exp(this._sigma * t) / Math.PI;
}
return yArr;
}
//-------------------------------------------------------------------------------------------------
}
}
Talbot algorithm
/*
Copyright (C) 2017 Torsten Brischalle
email: torsten@brischalle.de
web: http://www.aaabbb.de
ported from:
numerical inverse laplace transform algorithm by
Fernando Damian Nieuwveldt
email:fdnieuwveldt@gmail.com
web: http://code.activestate.com/recipes/576934-numerical-inversion-of-the-laplace-transform-using/
Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation
files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy,
modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software
is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR
IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
using System;
using System.Collections.Generic;
using System.Linq;
using System.Numerics;
using System.Text;
using System.Threading.Tasks;
namespace LaplaceInversion
{
/*
# Talbot suggested that the Bromwich line be deformed into a contour that begins
# and ends in the left half plane, i.e., z → −∞ at both ends.
# Due to the exponential factor the integrand decays rapidly
# on such a contour. In such situations the trapezoidal rule converge
# extraordinarily rapidly.
# For example here we compute the inverse transform of F(s) = 1/(s+1) at t = 1
#
# >>> error = Talbot(1,24)-exp(-1)
# >>> error
# (3.3306690738754696e-015+0j)
#
# Talbot method is very powerful here we see an error of 3.3e-015
# with only 24 function evaluations
#
# Created by Fernando Damian Nieuwveldt
# email:fdnieuwveldt@gmail.com
# Date : 25 October 2009
#
# Reference
# L.N.Trefethen, J.A.C.Weideman, and T.Schmelzer. Talbot quadratures
# and rational approximations. BIT. Numerical Mathematics,
# 46(3):653 670, 2006.
*/
public class LaplaceInversionTalbot : ILaplaceInversion
{
private int _n;
private double _shift = 0.0;
public string Name { get { return "Talbot"; } }
//-------------------------------------------------------------------------------------------------
public LaplaceInversionTalbot(int n = 128)
{
_n = n;
}
//-------------------------------------------------------------------------------------------------
public double[] Calculate(LaplaceFunction F, double[] tArr)
{
double[] yArr = new double[tArr.Length];
int idx;
for (idx = 0; idx < tArr.Length; idx++)
{
double t = tArr[idx];
if (t <= 0)
{
yArr[idx] = 0;
continue;
}
// Initiate the stepsize
double h = 2 * Math.PI / _n;
// Shift contour to the right in case there is a pole on the positive real axis : Note the contour will
// not be optimal since it was originally devoloped for function with
// singularities on the negative real axis
// For example take F(s) = 1/(s-1), it has a pole at s = 1, the contour needs to be shifted with one
// unit, i.e shift = 1. But in the test example no shifting is necessary
Complex ans = new Complex(0, 0);
int k;
double c1 = 0.5017;
double c2 = 0.6407;
double c3 = 0.6122;
Complex c4 = new Complex(0, 0.2645);
// The for loop is evaluating the Laplace inversion at each point theta which is based on the trapezoidal rule
for (k = 0; k <= _n; k++)
{
double theta = -Math.PI + (k + 0.5) * h;
Complex z = _shift + _n / t * (c1 * theta / Math.Tan(c2 * theta) - c3 + c4 * theta);
Complex dz = _n / t * (-c1 * c2 * theta / Sqr(Math.Sin(c2 * theta)) + c1 / Math.Tan(c2 * theta) + c4);
ans += Complex.Exp(z * t) * F(z) * dz;
}
yArr[idx] = ((h / (2 * Complex.ImaginaryOne * Math.PI)) * ans).Real;
}
return yArr;
}
//-------------------------------------------------------------------------------------------------
private double Sqr(double d)
{
return d * d;
}
//-------------------------------------------------------------------------------------------------
}
}
Zakian algorithm
/*
Copyright (C) 2017 Torsten Brischalle
email: torsten@brischalle.de
web: http://www.aaabbb.de
ported from:
numerical inverse laplace transform algorithm by
Fausto Arinos de Almeida Barbuto (Calgary, Canada)
E-mail: fausto_barbuto@yahoo.ca
web: http://code.activestate.com/recipes/127469-numerical-inversion-of-laplace-transforms-through-/
Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation
files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy,
modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software
is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR
IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
using System;
using System.Collections.Generic;
using System.Linq;
using System.Numerics;
using System.Text;
using System.Threading.Tasks;
namespace LaplaceInversion
{
/*
#################################################################
# Function Zakian(t), numerically inverts a Laplace transform #
# F(s) into f(t) for a specific time "t" using the Zakian #
# algorithm. #
# #
# Function F(s) is defined in separate as Fs(s) (see code #
# below). In the present case, F(s) = 1/(s-1), whose #
# analytical inversion is f(t) = exp(t), and t=1.0 thus #
# yielding f(1) = 2.71828180837 what is "e" with an error #
# of -7.4e-007%. Oscillatory systems may become innacurate #
# after the second cycle, though. #
# #
# Fs(s) has to be changed accordingly everytime the user #
# wants to invert a different function. #
# #
# Creator: Fausto Arinos de Almeida Barbuto (Calgary, Canada) #
# Date: May 16, 2002 #
# E-mail: fausto_barbuto@yahoo.ca #
# #
# Reference: #
# Huddleston, T. and Byrne, P: "Numerical Inversion of #
# Laplace Transforms", University of South Alabama, April #
# 1999 (found at http://www.eng.usouthal.edu/huddleston/ #
# SoftwareSupport/Download/Inversion99.doc) #
# #
#################################################################
*/
public class LaplaceInversionZakian : ILaplaceInversion
{
private Complex[] _a;
private Complex[] _k;
public string Name { get { return "Zakian"; } }
//-------------------------------------------------------------------------------------------------
public LaplaceInversionZakian()
{
this._a = new Complex[] {
new Complex(12.83767675, 1.666063445),
new Complex(12.22613209, 5.012718792),
new Complex(10.93430308, 8.409673116),
new Complex(8.776434715, 11.92185389),
new Complex(5.225453361, 15.72952905)
};
this._k = new Complex[] {
new Complex(-36902.08210, 196990.4257),
new Complex(61277.02524, -95408.62551),
new Complex(-28916.56288, 18169.18531),
new Complex(4655.361138, -1.901528642),
new Complex(-118.7414011, -141.3036911)
};
}
//-------------------------------------------------------------------------------------------------
public double[] Calculate(LaplaceFunction F, double[] tArr)
{
double[] yArr = new double[tArr.Length];
int idx;
for (idx = 0; idx < tArr.Length; idx++)
{
double t = tArr[idx];
if (t <= 0)
{
yArr[idx] = 0;
continue;
}
double sum = 0;
int i;
Complex c = new Complex();
for (i = 0; i < 5; i++)
{
c = this._a[i] / t;
c = F(c) * this._k[i];
sum += c.Real;
}
yArr[idx] = 2 * sum / t;
}
return yArr;
}
//-------------------------------------------------------------------------------------------------
}
}
Stehfest algorithm
/*
Copyright (C) 2017 Torsten Brischalle
email: torsten@brischalle.de
web: http://www.aaabbb.de
Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation
files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy,
modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software
is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR
IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
using System;
using System.Collections.Generic;
using System.Linq;
using System.Numerics;
using System.Text;
using System.Threading.Tasks;
namespace LaplaceInversion
{
public class LaplaceInversionStehfest : ILaplaceInversion
{
private readonly int _N = 16;
private double _ln2;
private double[] _valArr;
public string Name { get { return "Stehfest"; } }
//-------------------------------------------------------------------------------------------------
public LaplaceInversionStehfest()
{
_ln2 = Math.Log(2.0);
_valArr = new double[_N];
int n2 = _N / 2;
double sign = 1;
int i;
if ((n2 % 2) != 0)
sign = -1;
for (i = 0; i < this._N; i++)
{
int kMin = (i + 2) / 2;
int kMax = i + 1;
int k;
if (kMax > n2)
kMax = n2;
_valArr[i] = 0;
sign = -sign;
for (k = kMin; k <= kMax; k++)
{
_valArr[i] += (Math.Pow(k, n2) / Factorial(k)) * (Factorial(2 * k)
/ Factorial(2 * k - i - 1)) / Factorial(n2 - k)
/ Factorial(k - 1) / Factorial(i + 1 - k);
}
_valArr[i] = sign * _valArr[i];
}
}
//-------------------------------------------------------------------------------------------------
public double[] Calculate(LaplaceFunction F, double[] tArr)
{
double[] yArr = new double[tArr.Length];
int idx;
for (idx = 0; idx < tArr.Length; idx++)
{
double t = tArr[idx];
if (t <= 0)
{
yArr[idx] = 0;
continue;
}
double ln2t = this._ln2 / t;
double x = 0;
double y = 0;
int i;
for (i = 0; i < this._valArr.Length; i++)
{
x += ln2t;
y += this._valArr[i] * F(x).Real;
}
yArr[idx] = ln2t * y;
}
return yArr;
}
//-------------------------------------------------------------------------------------------------
private double Factorial(int n)
{
double x = 1;
if (n > 1)
{
int i;
for (i = 2; i <= n; i++)
x *= i;
}
return x;
}
//-------------------------------------------------------------------------------------------------
}
}