1"""
2Ordinary Differential Equations (ODEs) - Numerical Methods and Analysis
3
4This script demonstrates:
5- Euler's method
6- Runge-Kutta 4th order (RK4)
7- scipy.integrate.solve_ivp
8- Harmonic oscillator (second-order ODE)
9- Damped oscillator
10- Lorenz system (chaotic system)
11- Phase portraits
12"""
13
14import numpy as np
15
16try:
17 from scipy.integrate import solve_ivp
18 HAS_SCIPY = True
19except ImportError:
20 HAS_SCIPY = False
21 print("Warning: scipy not available")
22
23try:
24 import matplotlib.pyplot as plt
25 from mpl_toolkits.mplot3d import Axes3D
26 HAS_MATPLOTLIB = True
27except ImportError:
28 HAS_MATPLOTLIB = False
29 print("Warning: matplotlib not available, skipping visualizations")
30
31
32def euler_method(f, y0, t_span, dt):
33 """
34 Solve ODE dy/dt = f(t, y) using Euler's method
35 y_{n+1} = y_n + dt * f(t_n, y_n)
36 """
37 t0, tf = t_span
38 t = np.arange(t0, tf + dt, dt)
39 y = np.zeros((len(t), len(np.atleast_1d(y0))))
40 y[0] = y0
41
42 for i in range(len(t) - 1):
43 y[i + 1] = y[i] + dt * f(t[i], y[i])
44
45 return t, y
46
47
48def rk4_method(f, y0, t_span, dt):
49 """
50 Solve ODE dy/dt = f(t, y) using Runge-Kutta 4th order
51 """
52 t0, tf = t_span
53 t = np.arange(t0, tf + dt, dt)
54 y = np.zeros((len(t), len(np.atleast_1d(y0))))
55 y[0] = y0
56
57 for i in range(len(t) - 1):
58 k1 = dt * f(t[i], y[i])
59 k2 = dt * f(t[i] + dt/2, y[i] + k1/2)
60 k3 = dt * f(t[i] + dt/2, y[i] + k2/2)
61 k4 = dt * f(t[i] + dt, y[i] + k3)
62
63 y[i + 1] = y[i] + (k1 + 2*k2 + 2*k3 + k4) / 6
64
65 return t, y
66
67
68def harmonic_oscillator(t, y, omega=1.0):
69 """
70 Simple harmonic oscillator: y'' + ω²y = 0
71 State: y = [y, y']
72 """
73 return np.array([y[1], -omega**2 * y[0]])
74
75
76def damped_oscillator(t, y, gamma=0.5, omega=1.0):
77 """
78 Damped harmonic oscillator: y'' + 2γy' + ω²y = 0
79 State: y = [y, y']
80 """
81 return np.array([y[1], -2*gamma*y[1] - omega**2 * y[0]])
82
83
84def lorenz_system(t, y, sigma=10, rho=28, beta=8/3):
85 """
86 Lorenz system (chaotic):
87 dx/dt = σ(y - x)
88 dy/dt = x(ρ - z) - y
89 dz/dt = xy - βz
90 """
91 x, y_val, z = y
92 return np.array([
93 sigma * (y_val - x),
94 x * (rho - z) - y_val,
95 x * y_val - beta * z
96 ])
97
98
99def van_der_pol(t, y, mu=1.0):
100 """
101 Van der Pol oscillator (limit cycle):
102 y'' - μ(1 - y²)y' + y = 0
103 """
104 return np.array([y[1], mu * (1 - y[0]**2) * y[1] - y[0]])
105
106
107def pendulum(t, y, g=9.8, L=1.0, damping=0.0):
108 """
109 Nonlinear pendulum: θ'' + (g/L)sin(θ) + damping*θ' = 0
110 State: y = [θ, θ']
111 """
112 return np.array([y[1], -(g/L) * np.sin(y[0]) - damping * y[1]])
113
114
115# ============================================================================
116# MAIN DEMONSTRATIONS
117# ============================================================================
118
119print("=" * 70)
120print("ORDINARY DIFFERENTIAL EQUATIONS - NUMERICAL METHODS")
121print("=" * 70)
122
123# Test 1: Exponential decay - Euler vs RK4
124print("\n1. EXPONENTIAL DECAY: dy/dt = -y, y(0) = 1")
125print("-" * 70)
126
127# Exact solution: y(t) = e^(-t)
128f_decay = lambda t, y: -y
129y0 = np.array([1.0])
130t_span = (0, 5)
131dt = 0.5
132
133t_euler, y_euler = euler_method(f_decay, y0, t_span, dt)
134t_rk4, y_rk4 = rk4_method(f_decay, y0, t_span, dt)
135y_exact = np.exp(-t_euler)
136
137print(f"Time step dt = {dt}")
138print(f"\n{'t':>6s} {'Exact':>10s} {'Euler':>10s} {'RK4':>10s} {'Euler Err':>12s} {'RK4 Err':>12s}")
139print("-" * 70)
140for i in range(0, len(t_euler), 2):
141 t = t_euler[i]
142 exact = y_exact[i]
143 euler_val = y_euler[i, 0]
144 rk4_val = y_rk4[i, 0]
145 euler_err = abs(exact - euler_val)
146 rk4_err = abs(exact - rk4_val)
147 print(f"{t:6.1f} {exact:10.6f} {euler_val:10.6f} {rk4_val:10.6f} "
148 f"{euler_err:12.2e} {rk4_err:12.2e}")
149
150# Test 2: Harmonic oscillator
151print("\n2. HARMONIC OSCILLATOR: y'' + ω²y = 0")
152print("-" * 70)
153omega = 2.0
154y0 = np.array([1.0, 0.0]) # y(0)=1, y'(0)=0
155t_span = (0, 10)
156dt = 0.05
157
158t_rk4, y_rk4 = rk4_method(lambda t, y: harmonic_oscillator(t, y, omega),
159 y0, t_span, dt)
160
161# Exact solution: y(t) = cos(ωt)
162y_exact = np.cos(omega * t_rk4)
163
164print(f"ω = {omega}, y(0) = {y0[0]}, y'(0) = {y0[1]}")
165print(f"\nEnergy conservation check (E = 0.5(y'² + ω²y²)):")
166for i in [0, len(t_rk4)//4, len(t_rk4)//2, -1]:
167 t = t_rk4[i]
168 y_val = y_rk4[i, 0]
169 y_dot = y_rk4[i, 1]
170 energy = 0.5 * (y_dot**2 + omega**2 * y_val**2)
171 print(f" t={t:5.2f}: E = {energy:.6f}")
172
173# Test 3: Damped oscillator
174print("\n3. DAMPED HARMONIC OSCILLATOR")
175print("-" * 70)
176gamma_vals = [0.1, 0.5, 1.0, 2.0]
177omega = 1.0
178y0 = np.array([1.0, 0.0])
179t_span = (0, 20)
180dt = 0.05
181
182print("Damping regimes:")
183for gamma in gamma_vals:
184 discriminant = gamma**2 - omega**2
185 if discriminant < 0:
186 regime = "Underdamped"
187 elif discriminant == 0:
188 regime = "Critically damped"
189 else:
190 regime = "Overdamped"
191
192 t_rk4, y_rk4 = rk4_method(lambda t, y: damped_oscillator(t, y, gamma, omega),
193 y0, t_span, dt)
194
195 # Find amplitude decay
196 max_val_initial = np.max(np.abs(y_rk4[:100, 0]))
197 max_val_final = np.max(np.abs(y_rk4[-100:, 0]))
198
199 print(f" γ={gamma:.1f}: {regime:18s}, "
200 f"initial amp={max_val_initial:.4f}, "
201 f"final amp={max_val_final:.4f}")
202
203# Test 4: Lorenz system
204print("\n4. LORENZ SYSTEM (CHAOTIC ATTRACTOR)")
205print("-" * 70)
206sigma, rho, beta = 10, 28, 8/3
207y0 = np.array([1.0, 1.0, 1.0])
208t_span = (0, 50)
209dt = 0.01
210
211if HAS_SCIPY:
212 sol = solve_ivp(lambda t, y: lorenz_system(t, y, sigma, rho, beta),
213 t_span, y0, dense_output=True, max_step=dt)
214 t_lorenz = sol.t
215 y_lorenz = sol.y.T
216else:
217 t_lorenz, y_lorenz = rk4_method(lambda t, y: lorenz_system(t, y, sigma, rho, beta),
218 y0, t_span, dt)
219
220print(f"σ={sigma}, ρ={rho}, β={beta:.3f}")
221print(f"Initial conditions: x={y0[0]}, y={y0[1]}, z={y0[2]}")
222print(f"\nTrajectory statistics:")
223print(f" x: min={np.min(y_lorenz[:, 0]):7.3f}, max={np.max(y_lorenz[:, 0]):7.3f}")
224print(f" y: min={np.min(y_lorenz[:, 1]):7.3f}, max={np.max(y_lorenz[:, 1]):7.3f}")
225print(f" z: min={np.min(y_lorenz[:, 2]):7.3f}, max={np.max(y_lorenz[:, 2]):7.3f}")
226
227# Test sensitivity to initial conditions
228y0_perturb = y0 + np.array([0.001, 0, 0])
229if HAS_SCIPY:
230 sol_perturb = solve_ivp(lambda t, y: lorenz_system(t, y, sigma, rho, beta),
231 t_span, y0_perturb, dense_output=True, max_step=dt)
232 y_perturb = sol_perturb.y.T
233else:
234 _, y_perturb = rk4_method(lambda t, y: lorenz_system(t, y, sigma, rho, beta),
235 y0_perturb, t_span, dt)
236
237divergence = np.linalg.norm(y_lorenz - y_perturb, axis=1)
238print(f"\nSensitivity to initial conditions (Δx₀ = 0.001):")
239print(f" t={t_lorenz[len(t_lorenz)//4]:.1f}: distance = {divergence[len(divergence)//4]:.4f}")
240print(f" t={t_lorenz[len(t_lorenz)//2]:.1f}: distance = {divergence[len(divergence)//2]:.4f}")
241
242# Test 5: Van der Pol oscillator
243print("\n5. VAN DER POL OSCILLATOR (LIMIT CYCLE)")
244print("-" * 70)
245mu_vals = [0.1, 1.0, 5.0]
246t_span = (0, 50)
247dt = 0.05
248
249print("Nonlinearity parameter μ:")
250for mu in mu_vals:
251 y0 = np.array([2.0, 0.0])
252 t_vdp, y_vdp = rk4_method(lambda t, y: van_der_pol(t, y, mu),
253 y0, t_span, dt)
254
255 # Check limit cycle amplitude (look at later times)
256 y_steady = y_vdp[len(y_vdp)//2:, 0]
257 amplitude = (np.max(y_steady) - np.min(y_steady)) / 2
258
259 print(f" μ={mu:.1f}: limit cycle amplitude ≈ {amplitude:.3f}")
260
261# Test 6: Nonlinear pendulum
262print("\n6. NONLINEAR PENDULUM")
263print("-" * 70)
264g, L = 9.8, 1.0
265t_span = (0, 10)
266dt = 0.01
267
268print("Initial angle θ₀:")
269for theta0 in [0.1, 1.0, 3.0]: # Small, medium, large angle
270 y0 = np.array([theta0, 0.0])
271 t_pend, y_pend = rk4_method(lambda t, y: pendulum(t, y, g, L),
272 y0, t_span, dt)
273
274 # Compute period (time between zero crossings from same direction)
275 crossings = []
276 for i in range(1, len(y_pend)):
277 if y_pend[i-1, 0] < 0 and y_pend[i, 0] >= 0 and y_pend[i, 1] > 0:
278 crossings.append(t_pend[i])
279
280 if len(crossings) >= 2:
281 period = np.mean(np.diff(crossings))
282 else:
283 period = np.nan
284
285 # Small angle approximation: T = 2π√(L/g)
286 period_approx = 2 * np.pi * np.sqrt(L / g)
287
288 print(f" θ₀={theta0:.1f} rad: T={period:.4f}s "
289 f"(small angle: {period_approx:.4f}s)")
290
291# Visualization
292if HAS_MATPLOTLIB:
293 print("\n" + "=" * 70)
294 print("VISUALIZATIONS")
295 print("=" * 70)
296
297 fig = plt.figure(figsize=(15, 10))
298
299 # Plot 1: Euler vs RK4 comparison
300 ax1 = plt.subplot(2, 3, 1)
301 f_decay = lambda t, y: -y
302 y0 = np.array([1.0])
303 t_span = (0, 5)
304 dt = 0.5
305
306 t_euler, y_euler = euler_method(f_decay, y0, t_span, dt)
307 t_rk4, y_rk4 = rk4_method(f_decay, y0, t_span, dt)
308 t_exact = np.linspace(0, 5, 100)
309 y_exact = np.exp(-t_exact)
310
311 ax1.plot(t_exact, y_exact, 'k-', linewidth=2, label='Exact')
312 ax1.plot(t_euler, y_euler[:, 0], 'ro--', markersize=6, label='Euler')
313 ax1.plot(t_rk4, y_rk4[:, 0], 'bs--', markersize=6, label='RK4')
314 ax1.set_xlabel('t')
315 ax1.set_ylabel('y')
316 ax1.set_title('Exponential Decay: dy/dt = -y')
317 ax1.legend()
318 ax1.grid(True, alpha=0.3)
319
320 # Plot 2: Harmonic oscillator
321 ax2 = plt.subplot(2, 3, 2)
322 omega = 2.0
323 y0 = np.array([1.0, 0.0])
324 t_span = (0, 10)
325 dt = 0.05
326
327 t_rk4, y_rk4 = rk4_method(lambda t, y: harmonic_oscillator(t, y, omega),
328 y0, t_span, dt)
329 y_exact = np.cos(omega * t_rk4)
330
331 ax2.plot(t_rk4, y_exact, 'k--', linewidth=2, alpha=0.5, label='Exact')
332 ax2.plot(t_rk4, y_rk4[:, 0], 'b-', linewidth=1.5, label='RK4')
333 ax2.set_xlabel('t')
334 ax2.set_ylabel('y')
335 ax2.set_title(f'Harmonic Oscillator (ω={omega})')
336 ax2.legend()
337 ax2.grid(True, alpha=0.3)
338
339 # Plot 3: Phase portrait - harmonic oscillator
340 ax3 = plt.subplot(2, 3, 3)
341 ax3.plot(y_rk4[:, 0], y_rk4[:, 1], 'b-', linewidth=1.5)
342 ax3.plot(y_rk4[0, 0], y_rk4[0, 1], 'go', markersize=8, label='Start')
343 ax3.set_xlabel('y')
344 ax3.set_ylabel("y'")
345 ax3.set_title('Phase Portrait: Harmonic Oscillator')
346 ax3.legend()
347 ax3.grid(True, alpha=0.3)
348 ax3.set_aspect('equal')
349
350 # Plot 4: Damped oscillators
351 ax4 = plt.subplot(2, 3, 4)
352 for gamma in [0.1, 0.5, 1.0, 2.0]:
353 t_rk4, y_rk4 = rk4_method(lambda t, y: damped_oscillator(t, y, gamma, 1.0),
354 np.array([1.0, 0.0]), (0, 20), 0.05)
355 ax4.plot(t_rk4, y_rk4[:, 0], linewidth=1.5, label=f'γ={gamma}')
356
357 ax4.set_xlabel('t')
358 ax4.set_ylabel('y')
359 ax4.set_title('Damped Oscillators')
360 ax4.legend()
361 ax4.grid(True, alpha=0.3)
362
363 # Plot 5: Lorenz attractor (3D)
364 ax5 = fig.add_subplot(2, 3, 5, projection='3d')
365 ax5.plot(y_lorenz[:, 0], y_lorenz[:, 1], y_lorenz[:, 2],
366 'b-', linewidth=0.5, alpha=0.7)
367 ax5.plot([y_lorenz[0, 0]], [y_lorenz[0, 1]], [y_lorenz[0, 2]],
368 'go', markersize=6)
369 ax5.set_xlabel('x')
370 ax5.set_ylabel('y')
371 ax5.set_zlabel('z')
372 ax5.set_title('Lorenz Attractor')
373
374 # Plot 6: Van der Pol limit cycle
375 ax6 = plt.subplot(2, 3, 6)
376 for mu in [0.5, 2.0]:
377 t_vdp, y_vdp = rk4_method(lambda t, y: van_der_pol(t, y, mu),
378 np.array([2.0, 0.0]), (0, 50), 0.05)
379 # Plot only steady-state portion
380 ax6.plot(y_vdp[1000:, 0], y_vdp[1000:, 1],
381 linewidth=1.5, label=f'μ={mu}')
382
383 ax6.set_xlabel('y')
384 ax6.set_ylabel("y'")
385 ax6.set_title('Van der Pol Limit Cycles')
386 ax6.legend()
387 ax6.grid(True, alpha=0.3)
388 ax6.set_aspect('equal')
389
390 plt.tight_layout()
391 plt.savefig('/opt/projects/01_Personal/03_Study/examples/Mathematical_Methods/07_ode.png', dpi=150)
392 print("Saved visualization: 07_ode.png")
393 plt.close()
394
395print("\n" + "=" * 70)
396print("DEMONSTRATION COMPLETE")
397print("=" * 70)