1"""
2Laplace Transform - Transform Pairs, Inverse Transform, and ODE Solutions
3
4This script demonstrates:
5- Laplace transform pairs
6- Inverse Laplace transform (numerical)
7- Solving ODEs using Laplace transform
8- Transfer functions
9- Step response
10- Frequency response
11- Convolution theorem
12"""
13
14import numpy as np
15
16try:
17 from scipy.integrate import quad
18 from scipy.signal import lti, step, impulse, bode
19 HAS_SCIPY = True
20except ImportError:
21 HAS_SCIPY = False
22 print("Warning: scipy not available, using limited implementations")
23
24try:
25 import matplotlib.pyplot as plt
26 HAS_MATPLOTLIB = True
27except ImportError:
28 HAS_MATPLOTLIB = False
29 print("Warning: matplotlib not available, skipping visualizations")
30
31
32def laplace_transform(f, s, t_max=10, n_points=1000):
33 """
34 Compute Laplace transform: F(s) = ∫₀^∞ f(t)e^(-st) dt
35 Numerically integrate from 0 to t_max
36 """
37 if HAS_SCIPY:
38 integrand = lambda t: f(t) * np.exp(-s * t)
39 result, error = quad(integrand, 0, t_max)
40 return result
41 else:
42 t = np.linspace(0, t_max, n_points)
43 integrand = f(t) * np.exp(-s * t)
44 result = np.trapz(integrand, t)
45 return result
46
47
48def inverse_laplace_bromwich(F, t, gamma=0.5, n_points=100):
49 """
50 Inverse Laplace transform using Bromwich integral (simplified)
51 f(t) = (1/2πi) ∫_{γ-i∞}^{γ+i∞} F(s)e^(st) ds
52
53 Use finite limits for practical computation
54 """
55 omega = np.linspace(-50, 50, n_points)
56 s = gamma + 1j * omega
57 integrand = F(s) * np.exp(s * t)
58
59 # Trapezoidal integration
60 result = np.trapz(integrand, omega) / (2 * np.pi)
61 return result.real
62
63
64class LaplacePairs:
65 """Common Laplace transform pairs"""
66
67 @staticmethod
68 def unit_step():
69 """u(t) ↔ 1/s"""
70 f = lambda t: 1.0
71 F = lambda s: 1 / s
72 return f, F
73
74 @staticmethod
75 def exponential(a):
76 """e^(at) ↔ 1/(s-a)"""
77 f = lambda t: np.exp(a * t)
78 F = lambda s: 1 / (s - a)
79 return f, F
80
81 @staticmethod
82 def sine(omega):
83 """sin(ωt) ↔ ω/(s²+ω²)"""
84 f = lambda t: np.sin(omega * t)
85 F = lambda s: omega / (s**2 + omega**2)
86 return f, F
87
88 @staticmethod
89 def cosine(omega):
90 """cos(ωt) ↔ s/(s²+ω²)"""
91 f = lambda t: np.cos(omega * t)
92 F = lambda s: s / (s**2 + omega**2)
93 return f, F
94
95 @staticmethod
96 def damped_sine(a, omega):
97 """e^(-at)sin(ωt) ↔ ω/((s+a)²+ω²)"""
98 f = lambda t: np.exp(-a * t) * np.sin(omega * t)
99 F = lambda s: omega / ((s + a)**2 + omega**2)
100 return f, F
101
102 @staticmethod
103 def power(n):
104 """t^n ↔ n!/s^(n+1)"""
105 f = lambda t: t**n
106 F = lambda s: np.math.factorial(n) / s**(n+1)
107 return f, F
108
109
110def solve_ode_laplace_first_order(a, b, y0):
111 """
112 Solve first-order ODE: y' + ay = b using Laplace transform
113 Returns solution function y(t)
114 """
115 # Laplace transform: sY(s) - y(0) + aY(s) = b/s
116 # Y(s) = (y0 + b/s) / (s + a)
117
118 def y(t):
119 # Inverse transform: y(t) = y0*e^(-at) + (b/a)(1 - e^(-at))
120 return y0 * np.exp(-a * t) + (b / a) * (1 - np.exp(-a * t))
121
122 return y
123
124
125def solve_ode_laplace_second_order(a, b, c, y0, yp0):
126 """
127 Solve second-order ODE: y'' + ay' + by = c using Laplace transform
128 Returns solution function y(t)
129 """
130 # Laplace: s²Y(s) - sy(0) - y'(0) + a(sY(s) - y(0)) + bY(s) = c/s
131 # Y(s) = [sy(0) + y'(0) + ay(0) + c/s] / (s² + as + b)
132
133 # Characteristic equation: s² + as + b = 0
134 discriminant = a**2 - 4*b
135
136 if discriminant > 0:
137 # Overdamped
138 r1 = (-a + np.sqrt(discriminant)) / 2
139 r2 = (-a - np.sqrt(discriminant)) / 2
140
141 def y(t):
142 # Particular solution: y_p = c/b
143 y_p = c / b
144 # Homogeneous solution
145 A = (yp0 - r2*(y0 - y_p)) / (r1 - r2)
146 B = (r1*(y0 - y_p) - yp0) / (r1 - r2)
147 return A * np.exp(r1 * t) + B * np.exp(r2 * t) + y_p
148
149 elif discriminant == 0:
150 # Critically damped
151 r = -a / 2
152
153 def y(t):
154 y_p = c / b
155 A = y0 - y_p
156 B = yp0 - r * A
157 return (A + B * t) * np.exp(r * t) + y_p
158
159 else:
160 # Underdamped
161 sigma = -a / 2
162 omega = np.sqrt(-discriminant) / 2
163
164 def y(t):
165 y_p = c / b
166 A = y0 - y_p
167 B = (yp0 - sigma * A) / omega
168 return np.exp(sigma * t) * (A * np.cos(omega * t) + B * np.sin(omega * t)) + y_p
169
170 return y
171
172
173# ============================================================================
174# MAIN DEMONSTRATIONS
175# ============================================================================
176
177print("=" * 70)
178print("LAPLACE TRANSFORM - PAIRS, INVERSE, AND ODE SOLUTIONS")
179print("=" * 70)
180
181# Test 1: Verify transform pairs
182print("\n1. LAPLACE TRANSFORM PAIRS")
183print("-" * 70)
184
185pairs = [
186 ("Unit step", LaplacePairs.unit_step()),
187 ("e^(-2t)", LaplacePairs.exponential(-2)),
188 ("sin(3t)", LaplacePairs.sine(3)),
189 ("cos(2t)", LaplacePairs.cosine(2)),
190 ("t²", LaplacePairs.power(2)),
191]
192
193s_test = 2.0
194print(f"Testing at s = {s_test}")
195print(f"\n{'Function':20s} {'F(s) formula':20s} {'Numerical':15s}")
196print("-" * 70)
197
198for name, (f, F) in pairs:
199 F_formula = F(s_test)
200 F_numerical = laplace_transform(f, s_test, t_max=10)
201 print(f"{name:20s} {F_formula:20.6f} {F_numerical:15.6f}")
202
203# Test 2: Inverse Laplace transform
204print("\n2. INVERSE LAPLACE TRANSFORM")
205print("-" * 70)
206
207# F(s) = 1/(s+1) → f(t) = e^(-t)
208F_exp = lambda s: 1 / (s + 1)
209t_test = 1.0
210f_inverse = inverse_laplace_bromwich(F_exp, t_test, gamma=0.5)
211f_exact = np.exp(-t_test)
212
213print(f"F(s) = 1/(s+1), expected: f(t) = e^(-t)")
214print(f"At t = {t_test}:")
215print(f" Inverse transform: {f_inverse:.6f}")
216print(f" Exact: {f_exact:.6f}")
217print(f" Error: {abs(f_inverse - f_exact):.2e}")
218
219# F(s) = ω/(s²+ω²) → f(t) = sin(ωt)
220omega = 2.0
221F_sin = lambda s: omega / (s**2 + omega**2)
222t_test = np.pi / (2 * omega) # Peak of sine
223f_inverse = inverse_laplace_bromwich(F_sin, t_test, gamma=0.5)
224f_exact = np.sin(omega * t_test)
225
226print(f"\nF(s) = {omega}/(s²+{omega**2}), expected: f(t) = sin({omega}t)")
227print(f"At t = {t_test:.4f}:")
228print(f" Inverse transform: {f_inverse:.6f}")
229print(f" Exact: {f_exact:.6f}")
230
231# Test 3: Solving first-order ODE
232print("\n3. SOLVING FIRST-ORDER ODE: y' + 2y = 4, y(0) = 0")
233print("-" * 70)
234
235a, b, y0 = 2.0, 4.0, 0.0
236y_solution = solve_ode_laplace_first_order(a, b, y0)
237
238print(f"ODE: y' + {a}y = {b}")
239print(f"Initial condition: y(0) = {y0}")
240print(f"\nSolution: y(t) = 2(1 - e^(-2t))")
241
242t_vals = [0, 0.5, 1.0, 2.0, 5.0]
243print(f"\n{'t':>6s} {'y(t)':>12s}")
244print("-" * 20)
245for t in t_vals:
246 print(f"{t:6.1f} {y_solution(t):12.6f}")
247
248print(f"\nSteady-state (t→∞): y = {b/a:.6f}")
249
250# Test 4: Solving second-order ODE
251print("\n4. SOLVING SECOND-ORDER ODE: y'' + 3y' + 2y = 4, y(0)=0, y'(0)=0")
252print("-" * 70)
253
254a, b, c = 3.0, 2.0, 4.0
255y0, yp0 = 0.0, 0.0
256y_solution = solve_ode_laplace_second_order(a, b, c, y0, yp0)
257
258print(f"ODE: y'' + {a}y' + {b}y = {c}")
259print(f"Initial conditions: y(0) = {y0}, y'(0) = {yp0}")
260
261# Characteristic roots
262discriminant = a**2 - 4*b
263r1 = (-a + np.sqrt(discriminant)) / 2
264r2 = (-a - np.sqrt(discriminant)) / 2
265print(f"Characteristic roots: r₁ = {r1:.4f}, r₂ = {r2:.4f}")
266print(f"System is overdamped")
267
268t_vals = [0, 0.5, 1.0, 2.0, 5.0]
269print(f"\n{'t':>6s} {'y(t)':>12s}")
270print("-" * 20)
271for t in t_vals:
272 print(f"{t:6.1f} {y_solution(t):12.6f}")
273
274print(f"\nSteady-state (t→∞): y = {c/b:.6f}")
275
276# Test 5: Underdamped oscillator
277print("\n5. UNDERDAMPED OSCILLATOR: y'' + 0.4y' + 4y = 0, y(0)=1, y'(0)=0")
278print("-" * 70)
279
280a, b, c = 0.4, 4.0, 0.0
281y0, yp0 = 1.0, 0.0
282y_solution = solve_ode_laplace_second_order(a, b, c, y0, yp0)
283
284discriminant = a**2 - 4*b
285print(f"Discriminant: {discriminant:.4f} < 0 (underdamped)")
286
287sigma = -a / 2
288omega_d = np.sqrt(-discriminant) / 2
289print(f"Damping coefficient: σ = {sigma:.4f}")
290print(f"Damped frequency: ωd = {omega_d:.4f}")
291
292t_vals = [0, 0.5, 1.0, 2.0, 5.0]
293print(f"\n{'t':>6s} {'y(t)':>12s}")
294print("-" * 20)
295for t in t_vals:
296 print(f"{t:6.1f} {y_solution(t):12.6f}")
297
298# Test 6: Transfer function (if scipy available)
299if HAS_SCIPY:
300 print("\n6. TRANSFER FUNCTION AND FREQUENCY RESPONSE")
301 print("-" * 70)
302
303 # Second-order system: H(s) = ω_n² / (s² + 2ζω_n s + ω_n²)
304 omega_n = 2.0 # Natural frequency
305 zeta = 0.3 # Damping ratio
306
307 num = [omega_n**2]
308 den = [1, 2*zeta*omega_n, omega_n**2]
309
310 print(f"Transfer function: H(s) = {omega_n**2}/(s² + {2*zeta*omega_n}s + {omega_n**2})")
311 print(f"Natural frequency ωn = {omega_n} rad/s")
312 print(f"Damping ratio ζ = {zeta}")
313
314 # Create system
315 system = lti(num, den)
316
317 # Peak frequency
318 if zeta < 1/np.sqrt(2):
319 omega_peak = omega_n * np.sqrt(1 - 2*zeta**2)
320 print(f"Resonant peak at ω = {omega_peak:.4f} rad/s")
321
322# Test 7: Convolution theorem
323print("\n7. CONVOLUTION THEOREM")
324print("-" * 70)
325
326print("Theorem: L{f*g} = F(s)·G(s)")
327print("Example: f(t) = e^(-t), g(t) = e^(-2t)")
328
329# Functions
330f = lambda t: np.exp(-t)
331g = lambda t: np.exp(-2*t)
332
333# Laplace transforms
334F = lambda s: 1 / (s + 1)
335G = lambda s: 1 / (s + 2)
336
337# Product in s-domain
338FG = lambda s: F(s) * G(s)
339
340# Analytical convolution: (f*g)(t) = ∫₀^t e^(-τ)e^(-2(t-τ))dτ = e^(-t) - e^(-2t)
341convolution_exact = lambda t: np.exp(-t) - np.exp(-2*t)
342
343# Verify at s=3
344s_test = 3.0
345FG_value = FG(s_test)
346L_convolution = laplace_transform(convolution_exact, s_test, t_max=10)
347
348print(f"\nAt s = {s_test}:")
349print(f" F(s)·G(s) = {FG_value:.6f}")
350print(f" L{{f*g}} = {L_convolution:.6f}")
351print(f" Error: {abs(FG_value - L_convolution):.2e}")
352
353# Visualization
354if HAS_MATPLOTLIB:
355 print("\n" + "=" * 70)
356 print("VISUALIZATIONS")
357 print("=" * 70)
358
359 fig = plt.figure(figsize=(15, 10))
360
361 # Plot 1: Common transform pairs
362 ax1 = plt.subplot(2, 3, 1)
363 t = np.linspace(0, 5, 500)
364
365 functions = [
366 (lambda t: np.exp(-t), 'e^(-t)'),
367 (lambda t: np.exp(-2*t), 'e^(-2t)'),
368 (lambda t: t * np.exp(-t), 'te^(-t)'),
369 ]
370
371 for f, label in functions:
372 ax1.plot(t, f(t), linewidth=2, label=label)
373
374 ax1.set_xlabel('t')
375 ax1.set_ylabel('f(t)')
376 ax1.set_title('Time Domain Functions')
377 ax1.legend()
378 ax1.grid(True, alpha=0.3)
379
380 # Plot 2: Step response of first-order system
381 ax2 = plt.subplot(2, 3, 2)
382 t = np.linspace(0, 5, 500)
383
384 for a in [0.5, 1.0, 2.0]:
385 y_sol = solve_ode_laplace_first_order(a, 1.0, 0.0)
386 y = np.array([y_sol(ti) for ti in t])
387 ax2.plot(t, y, linewidth=2, label=f'a={a}')
388
389 ax2.set_xlabel('t')
390 ax2.set_ylabel('y(t)')
391 ax2.set_title("First-Order Step Response: y' + ay = 1")
392 ax2.legend()
393 ax2.grid(True, alpha=0.3)
394
395 # Plot 3: Second-order responses
396 ax3 = plt.subplot(2, 3, 3)
397 t = np.linspace(0, 10, 500)
398
399 # Different damping ratios
400 configs = [
401 (0.1, 'Underdamped ζ=0.05'),
402 (0.4, 'Underdamped ζ=0.1'),
403 (4.0, 'Critically damped'),
404 ]
405
406 for a, label in configs:
407 y_sol = solve_ode_laplace_second_order(a, 4.0, 4.0, 0.0, 0.0)
408 y = np.array([y_sol(ti) for ti in t])
409 ax3.plot(t, y, linewidth=2, label=label)
410
411 ax3.axhline(y=1.0, color='k', linestyle='--', alpha=0.3)
412 ax3.set_xlabel('t')
413 ax3.set_ylabel('y(t)')
414 ax3.set_title("Second-Order Step Response")
415 ax3.legend()
416 ax3.grid(True, alpha=0.3)
417
418 # Plot 4: Damped oscillation
419 ax4 = plt.subplot(2, 3, 4)
420 t = np.linspace(0, 15, 500)
421
422 y_sol = solve_ode_laplace_second_order(0.4, 4.0, 0.0, 1.0, 0.0)
423 y = np.array([y_sol(ti) for ti in t])
424
425 # Envelope
426 sigma = -0.4 / 2
427 envelope_pos = np.exp(sigma * t)
428 envelope_neg = -np.exp(sigma * t)
429
430 ax4.plot(t, y, 'b-', linewidth=2, label='y(t)')
431 ax4.plot(t, envelope_pos, 'r--', alpha=0.5, label='Envelope')
432 ax4.plot(t, envelope_neg, 'r--', alpha=0.5)
433 ax4.set_xlabel('t')
434 ax4.set_ylabel('y(t)')
435 ax4.set_title('Underdamped Oscillator')
436 ax4.legend()
437 ax4.grid(True, alpha=0.3)
438
439 # Plot 5: Frequency response (if scipy available)
440 if HAS_SCIPY:
441 ax5 = plt.subplot(2, 3, 5)
442
443 omega_n = 2.0
444 zeta = 0.3
445 num = [omega_n**2]
446 den = [1, 2*zeta*omega_n, omega_n**2]
447 system = lti(num, den)
448
449 w, mag, phase = bode(system, np.logspace(-1, 1, 100))
450
451 ax5.semilogx(w, mag, 'b-', linewidth=2)
452 ax5.set_xlabel('Frequency (rad/s)')
453 ax5.set_ylabel('Magnitude (dB)')
454 ax5.set_title(f'Bode Plot: ζ={zeta}, ωn={omega_n}')
455 ax5.grid(True, alpha=0.3, which='both')
456
457 # Plot 6: Phase response
458 ax6 = plt.subplot(2, 3, 6)
459 ax6.semilogx(w, phase, 'r-', linewidth=2)
460 ax6.set_xlabel('Frequency (rad/s)')
461 ax6.set_ylabel('Phase (deg)')
462 ax6.set_title('Phase Response')
463 ax6.grid(True, alpha=0.3, which='both')
464 else:
465 # Alternative plots
466 ax5 = plt.subplot(2, 3, 5)
467 t = np.linspace(0, 5, 500)
468 f1 = np.exp(-t)
469 f2 = np.exp(-2*t)
470
471 # Numerical convolution
472 dt = t[1] - t[0]
473 conv = np.convolve(f1, f2, mode='same') * dt
474
475 ax5.plot(t, conv, 'b-', linewidth=2, label='f*g (numerical)')
476 ax5.plot(t, np.exp(-t) - np.exp(-2*t), 'r--', linewidth=2, label='Analytical')
477 ax5.set_xlabel('t')
478 ax5.set_ylabel('(f*g)(t)')
479 ax5.set_title('Convolution: e^(-t) * e^(-2t)')
480 ax5.legend()
481 ax5.grid(True, alpha=0.3)
482
483 plt.tight_layout()
484 plt.savefig('/opt/projects/01_Personal/03_Study/examples/Mathematical_Methods/12_laplace_transform.png', dpi=150)
485 print("Saved visualization: 12_laplace_transform.png")
486 plt.close()
487
488print("\n" + "=" * 70)
489print("DEMONSTRATION COMPLETE")
490print("=" * 70)