Skip to content

Commit abcbec3

Browse files
feat: add numerical laplace transform
1 parent 791deb4 commit abcbec3

1 file changed

Lines changed: 63 additions & 0 deletions

File tree

maths/laplace_transformation.py

Lines changed: 63 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,63 @@
1+
"""
2+
The Laplace Transform is defined as: L{f(t)} = integral from 0 to infinity of e^(-st) * f(t) dt.
3+
4+
Wiki: https://en.wikipedia.org/wiki/Laplace_transform
5+
6+
"""
7+
8+
import numpy as np
9+
10+
11+
def laplace_transform(
12+
function_values: np.ndarray, s_value: float, delta_t: float
13+
) -> float:
14+
"""
15+
Calculate the numerical Laplace Transform of a function given its values over time.
16+
17+
Args:
18+
function_values: A numpy array of the function values f(t).
19+
s_value: The complex frequency parameter 's' (modeled here as a float).
20+
delta_t: The time step between samples.
21+
22+
Returns:
23+
The approximate value of the Laplace transform at s_value.
24+
25+
Example: For f(t) = 1, the Laplace transform L{1} = 1/s.
26+
If s = 2, L{1} should be 0.5.
27+
28+
>>> t = np.linspace(0, 50, 10000)
29+
>>> f_t = np.ones_like(t) # f(t) = 1
30+
>>> res = laplace_transform(f_t, s_value=2.0, delta_t=50/10000)
31+
>>> abs(res - 0.5) < 1e-3
32+
True
33+
34+
Example: For f(t) = e^(-t), the Laplace transform L{e^-t} = 1/(s+1).
35+
If s = 1, L{e^-t} should be 0.5.
36+
37+
>>> t = np.linspace(0, 50, 10000)
38+
>>> f_t = np.exp(-t)
39+
>>> res = laplace_transform(f_t, s_value=1.0, delta_t=50/10000)
40+
>>> abs(res - 0.5) < 1e-3
41+
True
42+
"""
43+
if s_value < 0:
44+
raise ValueError("s_value must be non-negative for convergence.")
45+
46+
# Time vector corresponding to the function values
47+
time_vector = np.arange(len(function_values)) * delta_t
48+
49+
# The integrand: f(t) * e^(-s*t)
50+
integrand = function_values * np.exp(-s_value * time_vector)
51+
52+
# Numerical integration using the trapezoidal rule
53+
result = np.trapezoid(integrand, dx=delta_t)
54+
55+
return float(result)
56+
57+
58+
if __name__ == "__main__":
59+
import doctest
60+
61+
62+
doctest.testmod()
63+

0 commit comments

Comments
 (0)