diff --git a/maths/laplace_transformation.py b/maths/laplace_transformation.py new file mode 100644 index 000000000000..089d681cee4f --- /dev/null +++ b/maths/laplace_transformation.py @@ -0,0 +1,66 @@ +""" +This module provides a numerical implementation of the Laplace Transform. + +https://en.wikipedia.org/wiki/Laplace_transform + +""" + +import numpy as np + + +def laplace_transform( + function_values: np.ndarray, s_value: float, delta_t: float +) -> float: + """ + Calculate the numerical Laplace Transform of a function given its values over time. + + This implementation supports only real-valued, non-negative Laplace + parameters ``s``. + + Args: + function_values: A numpy array of the function values f(t). + s_value: The real-valued Laplace parameter ``s``. Must be non-negative. + delta_t: The time step between samples. + + Returns: + The approximate real-valued value of the Laplace transform at s_value. + + Example: For f(t) = 1, the Laplace transform L{1} = 1/s. + If s = 2, L{1} should be 0.5. + + >>> t = np.linspace(0, 50, 10000) + >>> f_t = np.ones_like(t) # f(t) = 1 + >>> res = laplace_transform(f_t, s_value=2.0, delta_t=50/10000) + >>> abs(res - 0.5) < 1e-3 + True + + Example: For f(t) = e^(-t), the Laplace transform L{e^-t} = 1/(s+1). + If s = 1, L{e^-t} should be 0.5. + + >>> t = np.linspace(0, 50, 10000) + >>> f_t = np.exp(-t) + >>> res = laplace_transform(f_t, s_value=1.0, delta_t=50/10000) + >>> abs(res - 0.5) < 1e-3 + True + """ + if delta_t <= 0: + raise ValueError("delta_t must be a positive value.") + if function_values.size == 0: + raise ValueError("function_values array cannot be empty.") + + # Time vector corresponding to the function values + time_vector = np.arange(len(function_values)) * delta_t + + # The integrand: f(t) * e^(-s*t) + integrand = function_values * np.exp(-s_value * time_vector) + + # Numerical integration using the trapezoidal rule + result = np.trapezoid(integrand, dx=delta_t) + + return float(result) + + +if __name__ == "__main__": + import doctest + + doctest.testmod()