1"""
2Calculus of Variations - Euler-Lagrange Equation and Applications
3
4This script demonstrates:
5- Euler-Lagrange equation
6- Brachistochrone problem
7- Catenary curve
8- Geodesics on surfaces
9- Lagrangian mechanics (pendulum, spring)
10- Minimal surface of revolution
11- Isoperimetric problems
12"""
13
14import numpy as np
15
16try:
17 from scipy.optimize import minimize, fsolve
18 from scipy.integrate import odeint, solve_ivp
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 from mpl_toolkits.mplot3d import Axes3D
27 HAS_MATPLOTLIB = True
28except ImportError:
29 HAS_MATPLOTLIB = False
30 print("Warning: matplotlib not available, skipping visualizations")
31
32
33def euler_lagrange_geodesic_plane(x0, y0, x1, y1, n_points=100):
34 """
35 Solve for geodesic (straight line) on a plane
36 Minimize: ∫ √(1 + y'²) dx
37 Euler-Lagrange: d/dx(y'/√(1+y'²)) = 0 → y'' = 0
38 """
39 x = np.linspace(x0, x1, n_points)
40 # Straight line
41 y = y0 + (y1 - y0) * (x - x0) / (x1 - x0)
42 return x, y
43
44
45def brachistochrone_curve(x0, y0, x1, y1, n_points=100):
46 """
47 Brachistochrone problem: fastest descent curve under gravity
48 Solution is a cycloid: x = r(θ - sin(θ)), y = r(1 - cos(θ))
49 """
50 # Find r and theta_max that satisfy boundary conditions
51 # This is a transcendental equation, solve numerically
52
53 def equations(params):
54 r, theta_max = params
55 x_end = r * (theta_max - np.sin(theta_max))
56 y_end = r * (1 - np.cos(theta_max))
57 return [x_end - (x1 - x0), y_end - (y1 - y0)]
58
59 if HAS_SCIPY:
60 r_init = (x1 - x0) / np.pi
61 theta_init = np.pi
62 r, theta_max = fsolve(equations, [r_init, theta_init])
63
64 theta = np.linspace(0, theta_max, n_points)
65 x = x0 + r * (theta - np.sin(theta))
66 y = y0 + r * (1 - np.cos(theta))
67 else:
68 # Approximate with straight line
69 x = np.linspace(x0, x1, n_points)
70 y = y0 + (y1 - y0) * (x - x0) / (x1 - x0)
71
72 return x, y
73
74
75def catenary_curve(x0, x1, length, n_points=100):
76 """
77 Catenary: shape of hanging chain
78 Minimize: ∫ y√(1 + y'²) dx subject to fixed arc length
79 Solution: y = a*cosh(x/a) + b
80 """
81 # Find parameter 'a' such that arc length matches
82 def arc_length(a):
83 x = np.linspace(x0, x1, 1000)
84 y = a * np.cosh(x / a)
85 dy_dx = np.sinh(x / a)
86 ds = np.sqrt(1 + dy_dx**2)
87 return np.trapz(ds, x)
88
89 if HAS_SCIPY:
90 # Find 'a' that gives correct length
91 result = minimize(lambda a: (arc_length(a[0]) - length)**2, [1.0], bounds=[(0.1, 10)])
92 a = result.x[0]
93 else:
94 a = 1.0
95
96 x = np.linspace(x0, x1, n_points)
97 y = a * np.cosh(x / a)
98
99 return x, y
100
101
102def lagrangian_simple_pendulum(t, state, g=9.8, L=1.0):
103 """
104 Simple pendulum using Lagrangian mechanics
105 L = T - V = (1/2)mL²θ'² - mgL(1-cos(θ))
106 Euler-Lagrange: d/dt(∂L/∂θ') - ∂L/∂θ = 0
107 → θ'' + (g/L)sin(θ) = 0
108 """
109 theta, theta_dot = state
110 theta_ddot = -(g / L) * np.sin(theta)
111 return [theta_dot, theta_ddot]
112
113
114def lagrangian_spring_mass(t, state, k=1.0, m=1.0):
115 """
116 Spring-mass system using Lagrangian mechanics
117 L = T - V = (1/2)mx'² - (1/2)kx²
118 Euler-Lagrange: mx'' + kx = 0
119 """
120 x, x_dot = state
121 x_ddot = -(k / m) * x
122 return [x_dot, x_ddot]
123
124
125def minimal_surface_revolution(y0, y1, x_range, n_points=100):
126 """
127 Minimal surface of revolution
128 Minimize: ∫ y√(1 + y'²) dx
129 Euler-Lagrange: y'' = (1 + y'²)/y
130 Solution: catenoid y = c*cosh(x/c)
131 """
132 x0, x1 = x_range
133 x = np.linspace(x0, x1, n_points)
134
135 # For simplicity, use catenary with c=1
136 c = 1.0
137 y = c * np.cosh((x - (x0 + x1)/2) / c)
138
139 # Scale to match boundary conditions
140 y_scaled = y0 + (y1 - y0) * (y - y[0]) / (y[-1] - y[0])
141
142 return x, y_scaled
143
144
145# ============================================================================
146# MAIN DEMONSTRATIONS
147# ============================================================================
148
149print("=" * 70)
150print("CALCULUS OF VARIATIONS - EULER-LAGRANGE EQUATION")
151print("=" * 70)
152
153# Test 1: Geodesic on plane (straight line)
154print("\n1. GEODESIC ON PLANE - SHORTEST PATH")
155print("-" * 70)
156
157x0, y0 = 0, 0
158x1, y1 = 3, 2
159
160x_geo, y_geo = euler_lagrange_geodesic_plane(x0, y0, x1, y1)
161
162length_geo = np.sqrt((x1 - x0)**2 + (y1 - y0)**2)
163
164print(f"From ({x0}, {y0}) to ({x1}, {y1})")
165print(f"Geodesic length: {length_geo:.6f}")
166print(f"Expected (Euclidean): {length_geo:.6f}")
167
168# Compute actual path length
169dx = x_geo[1] - x_geo[0]
170dy_dx = np.gradient(y_geo, dx)
171ds = np.sqrt(1 + dy_dx**2)
172computed_length = np.trapz(ds, x_geo)
173print(f"Numerical verification: {computed_length:.6f}")
174
175# Test 2: Brachistochrone problem
176print("\n2. BRACHISTOCHRONE - FASTEST DESCENT")
177print("-" * 70)
178
179x0, y0 = 0, 0
180x1, y1 = 2, -1
181
182x_brach, y_brach = brachistochrone_curve(x0, y0, x1, y1)
183
184print(f"From ({x0}, {y0}) to ({x1}, {y1})")
185print("Solution is a cycloid curve")
186
187# Compute descent time (T = ∫ ds/v where v = √(2g|y|))
188g = 9.8
189dx = x_brach[1] - x_brach[0]
190dy_dx = np.gradient(y_brach, dx)
191ds = np.sqrt(dx**2 + (dy_dx * dx)**2)
192
193# Velocity v = √(2g|y - y0|)
194v = np.sqrt(2 * g * np.abs(y_brach - y0) + 1e-10)
195dt = ds / v
196time_brach = np.sum(dt)
197
198print(f"Descent time (brachistochrone): {time_brach:.6f} s")
199
200# Compare with straight line
201x_straight = np.linspace(x0, x1, len(x_brach))
202y_straight = y0 + (y1 - y0) * (x_straight - x0) / (x1 - x0)
203dy_dx_straight = np.gradient(y_straight, dx)
204ds_straight = np.sqrt(dx**2 + (dy_dx_straight * dx)**2)
205v_straight = np.sqrt(2 * g * np.abs(y_straight - y0) + 1e-10)
206dt_straight = ds_straight / v_straight
207time_straight = np.sum(dt_straight)
208
209print(f"Descent time (straight line): {time_straight:.6f} s")
210print(f"Time saved: {(time_straight - time_brach):.6f} s ({(time_straight/time_brach - 1)*100:.1f}%)")
211
212# Test 3: Catenary
213print("\n3. CATENARY - HANGING CHAIN")
214print("-" * 70)
215
216x0, x1 = -1, 1
217chain_length = 3.0
218
219x_cat, y_cat = catenary_curve(x0, x1, chain_length)
220
221print(f"Chain from x={x0} to x={x1}")
222print(f"Total length: {chain_length}")
223
224# Verify length
225dx = x_cat[1] - x_cat[0]
226dy_dx = np.gradient(y_cat, dx)
227ds = np.sqrt(1 + dy_dx**2)
228computed_length = np.trapz(ds, x_cat)
229
230print(f"Computed length: {computed_length:.6f}")
231print(f"Shape: y = a·cosh(x/a)")
232
233# Potential energy (proportional to ∫ y ds)
234potential_energy = np.trapz(y_cat * ds, x_cat)
235print(f"Potential energy (normalized): {potential_energy:.6f}")
236
237# Test 4: Lagrangian mechanics - Simple pendulum
238print("\n4. LAGRANGIAN MECHANICS - SIMPLE PENDULUM")
239print("-" * 70)
240
241g, L = 9.8, 1.0
242theta0 = np.pi / 4 # 45 degrees
243theta_dot0 = 0
244
245print(f"Pendulum length L = {L} m")
246print(f"Initial angle θ₀ = {np.degrees(theta0):.1f}°")
247print(f"Initial velocity θ'₀ = {theta_dot0} rad/s")
248
249if HAS_SCIPY:
250 t_span = (0, 5)
251 t_eval = np.linspace(0, 5, 500)
252
253 sol = solve_ivp(lagrangian_simple_pendulum, t_span, [theta0, theta_dot0],
254 t_eval=t_eval, args=(g, L))
255 t_pend = sol.t
256 theta = sol.y[0]
257 theta_dot = sol.y[1]
258
259 # Energy conservation check
260 T = 0.5 * L**2 * theta_dot**2 # Kinetic (with m=1)
261 V = g * L * (1 - np.cos(theta)) # Potential
262 E = T + V
263
264 print(f"\nEnergy conservation:")
265 print(f" E(t=0) = {E[0]:.6f} J")
266 print(f" E(t={t_eval[-1]}) = {E[-1]:.6f} J")
267 print(f" Relative change: {abs(E[-1] - E[0])/E[0] * 100:.2f}%")
268
269 # Small angle approximation period
270 period_small_angle = 2 * np.pi * np.sqrt(L / g)
271 print(f"\nSmall angle period: {period_small_angle:.4f} s")
272
273# Test 5: Lagrangian mechanics - Spring-mass
274print("\n5. LAGRANGIAN MECHANICS - SPRING-MASS SYSTEM")
275print("-" * 70)
276
277k, m = 4.0, 1.0
278x0 = 1.0
279v0 = 0.0
280
281omega = np.sqrt(k / m)
282period = 2 * np.pi / omega
283
284print(f"Spring constant k = {k} N/m")
285print(f"Mass m = {m} kg")
286print(f"Natural frequency ω = {omega:.4f} rad/s")
287print(f"Period T = {period:.4f} s")
288
289if HAS_SCIPY:
290 t_span = (0, 10)
291 t_eval = np.linspace(0, 10, 500)
292
293 sol = solve_ivp(lagrangian_spring_mass, t_span, [x0, v0],
294 t_eval=t_eval, args=(k, m))
295 t_spring = sol.t
296 x = sol.y[0]
297 x_dot = sol.y[1]
298
299 # Energy conservation
300 T = 0.5 * m * x_dot**2
301 V = 0.5 * k * x**2
302 E = T + V
303
304 print(f"\nEnergy conservation:")
305 print(f" E(t=0) = {E[0]:.6f} J")
306 print(f" E(t={t_eval[-1]}) = {E[-1]:.6f} J")
307 print(f" Relative change: {abs(E[-1] - E[0])/E[0] * 100:.2f}%")
308
309# Test 6: Minimal surface of revolution
310print("\n6. MINIMAL SURFACE OF REVOLUTION - CATENOID")
311print("-" * 70)
312
313y0, y1 = 1.0, 1.0
314x_range = (-1, 1)
315
316x_surf, y_surf = minimal_surface_revolution(y0, y1, x_range)
317
318print(f"Boundary conditions: y({x_range[0]}) = {y0}, y({x_range[1]}) = {y1}")
319print("Rotating curve around x-axis")
320print("Solution: catenoid (minimal surface)")
321
322# Surface area = 2π ∫ y√(1 + y'²) dx
323dx = x_surf[1] - x_surf[0]
324dy_dx = np.gradient(y_surf, dx)
325integrand = y_surf * np.sqrt(1 + dy_dx**2)
326surface_area = 2 * np.pi * np.trapz(integrand, x_surf)
327
328print(f"Surface area: {surface_area:.6f}")
329
330# Visualization
331if HAS_MATPLOTLIB:
332 print("\n" + "=" * 70)
333 print("VISUALIZATIONS")
334 print("=" * 70)
335
336 fig = plt.figure(figsize=(15, 10))
337
338 # Plot 1: Geodesic
339 ax1 = plt.subplot(2, 3, 1)
340 ax1.plot(x_geo, y_geo, 'b-', linewidth=2, label='Geodesic')
341 ax1.plot([x0, x1], [y0, y1], 'ro', markersize=8)
342 ax1.set_xlabel('x')
343 ax1.set_ylabel('y')
344 ax1.set_title('Geodesic on Plane')
345 ax1.legend()
346 ax1.grid(True, alpha=0.3)
347 ax1.set_aspect('equal')
348
349 # Plot 2: Brachistochrone
350 ax2 = plt.subplot(2, 3, 2)
351 ax2.plot(x_brach, y_brach, 'b-', linewidth=2, label='Brachistochrone')
352 ax2.plot(x_straight, y_straight, 'r--', linewidth=2, label='Straight line')
353 ax2.plot([x0, x1], [y0, y1], 'ko', markersize=8)
354 ax2.set_xlabel('x')
355 ax2.set_ylabel('y')
356 ax2.set_title('Brachistochrone: Fastest Descent')
357 ax2.legend()
358 ax2.grid(True, alpha=0.3)
359
360 # Plot 3: Catenary
361 ax3 = plt.subplot(2, 3, 3)
362 ax3.plot(x_cat, y_cat, 'b-', linewidth=2, label='Catenary')
363 ax3.plot([x0, x1], [y_cat[0], y_cat[-1]], 'ro', markersize=8)
364 ax3.set_xlabel('x')
365 ax3.set_ylabel('y')
366 ax3.set_title('Catenary: Hanging Chain')
367 ax3.legend()
368 ax3.grid(True, alpha=0.3)
369 ax3.invert_yaxis()
370
371 # Plot 4: Pendulum motion
372 if HAS_SCIPY:
373 ax4 = plt.subplot(2, 3, 4)
374 ax4.plot(t_pend, np.degrees(theta), 'b-', linewidth=2)
375 ax4.set_xlabel('Time (s)')
376 ax4.set_ylabel('Angle (degrees)')
377 ax4.set_title('Simple Pendulum: θ(t)')
378 ax4.grid(True, alpha=0.3)
379
380 # Plot 5: Pendulum phase portrait
381 ax5 = plt.subplot(2, 3, 5)
382 ax5.plot(np.degrees(theta), np.degrees(theta_dot), 'b-', linewidth=1.5)
383 ax5.plot(np.degrees(theta[0]), np.degrees(theta_dot[0]), 'go', markersize=8)
384 ax5.set_xlabel('θ (degrees)')
385 ax5.set_ylabel("θ' (degrees/s)")
386 ax5.set_title('Pendulum Phase Portrait')
387 ax5.grid(True, alpha=0.3)
388
389 # Plot 6: Spring-mass energy
390 ax6 = plt.subplot(2, 3, 6)
391 ax6.plot(t_spring, T, 'b-', linewidth=2, label='Kinetic')
392 ax6.plot(t_spring, V, 'r-', linewidth=2, label='Potential')
393 ax6.plot(t_spring, E, 'k--', linewidth=2, label='Total')
394 ax6.set_xlabel('Time (s)')
395 ax6.set_ylabel('Energy (J)')
396 ax6.set_title('Spring-Mass Energy Conservation')
397 ax6.legend()
398 ax6.grid(True, alpha=0.3)
399
400 plt.tight_layout()
401 plt.savefig('/opt/projects/01_Personal/03_Study/examples/Mathematical_Methods/14_calculus_of_variations.png', dpi=150)
402 print("Saved visualization: 14_calculus_of_variations.png")
403 plt.close()
404
405 # Additional 3D plot: Catenoid surface
406 if HAS_SCIPY:
407 fig = plt.figure(figsize=(10, 8))
408 ax = fig.add_subplot(111, projection='3d')
409
410 # Create surface of revolution
411 theta_surf = np.linspace(0, 2*np.pi, 50)
412 X_surf, Theta_surf = np.meshgrid(x_surf, theta_surf)
413
414 Y_surf_2d = np.tile(y_surf, (len(theta_surf), 1))
415 Y_surf_3d = Y_surf_2d * np.cos(Theta_surf)
416 Z_surf_3d = Y_surf_2d * np.sin(Theta_surf)
417
418 ax.plot_surface(X_surf, Y_surf_3d, Z_surf_3d, cmap='viridis', alpha=0.8)
419 ax.set_xlabel('x')
420 ax.set_ylabel('y')
421 ax.set_zlabel('z')
422 ax.set_title('Catenoid: Minimal Surface of Revolution')
423
424 plt.savefig('/opt/projects/01_Personal/03_Study/examples/Mathematical_Methods/14_catenoid_3d.png', dpi=150)
425 print("Saved visualization: 14_catenoid_3d.png")
426 plt.close()
427
428print("\n" + "=" * 70)
429print("DEMONSTRATION COMPLETE")
430print("=" * 70)