1"""
2Vector Analysis - Gradient, Divergence, Curl, and Integral Theorems
3
4This script demonstrates:
5- Gradient of scalar fields
6- Divergence of vector fields
7- Curl of vector fields
8- Line integrals
9- Surface integrals
10- Green's theorem
11- Stokes' theorem
12- Divergence theorem verification
13"""
14
15import numpy as np
16
17try:
18 import matplotlib.pyplot as plt
19 from mpl_toolkits.mplot3d import Axes3D
20 HAS_MATPLOTLIB = True
21except ImportError:
22 HAS_MATPLOTLIB = False
23 print("Warning: matplotlib not available, skipping visualizations")
24
25
26def gradient_2d(f, x, y, h=1e-5):
27 """Compute gradient of scalar field f using finite differences"""
28 df_dx = (f(x + h, y) - f(x - h, y)) / (2 * h)
29 df_dy = (f(x, y + h) - f(x, y - h)) / (2 * h)
30 return np.array([df_dx, df_dy])
31
32
33def gradient_3d(f, x, y, z, h=1e-5):
34 """Compute gradient of scalar field f(x,y,z) using finite differences"""
35 df_dx = (f(x + h, y, z) - f(x - h, y, z)) / (2 * h)
36 df_dy = (f(x, y + h, z) - f(x, y - h, z)) / (2 * h)
37 df_dz = (f(x, y, z + h) - f(x, y, z - h)) / (2 * h)
38 return np.array([df_dx, df_dy, df_dz])
39
40
41def divergence_2d(Fx, Fy, x, y, h=1e-5):
42 """Compute divergence of 2D vector field F = (Fx, Fy)"""
43 dFx_dx = (Fx(x + h, y) - Fx(x - h, y)) / (2 * h)
44 dFy_dy = (Fy(x, y + h) - Fy(x, y - h)) / (2 * h)
45 return dFx_dx + dFy_dy
46
47
48def divergence_3d(Fx, Fy, Fz, x, y, z, h=1e-5):
49 """Compute divergence of 3D vector field F = (Fx, Fy, Fz)"""
50 dFx_dx = (Fx(x + h, y, z) - Fx(x - h, y, z)) / (2 * h)
51 dFy_dy = (Fy(x, y + h, z) - Fy(x, y - h, z)) / (2 * h)
52 dFz_dz = (Fz(x, y, z + h) - Fz(x, y, z - h)) / (2 * h)
53 return dFx_dx + dFy_dy + dFz_dz
54
55
56def curl_3d(Fx, Fy, Fz, x, y, z, h=1e-5):
57 """Compute curl of 3D vector field F = (Fx, Fy, Fz)"""
58 dFz_dy = (Fz(x, y + h, z) - Fz(x, y - h, z)) / (2 * h)
59 dFy_dz = (Fy(x, y, z + h) - Fy(x, y, z - h)) / (2 * h)
60
61 dFx_dz = (Fx(x, y, z + h) - Fx(x, y, z - h)) / (2 * h)
62 dFz_dx = (Fz(x + h, y, z) - Fz(x - h, y, z)) / (2 * h)
63
64 dFy_dx = (Fy(x + h, y, z) - Fy(x - h, y, z)) / (2 * h)
65 dFx_dy = (Fx(x, y + h, z) - Fx(x, y - h, z)) / (2 * h)
66
67 curl_x = dFz_dy - dFy_dz
68 curl_y = dFx_dz - dFz_dx
69 curl_z = dFy_dx - dFx_dy
70
71 return np.array([curl_x, curl_y, curl_z])
72
73
74def line_integral_2d(Fx, Fy, curve_x, curve_y, t_vals):
75 """Compute line integral of F·dr along parametric curve"""
76 integral = 0
77 for i in range(len(t_vals) - 1):
78 t = t_vals[i]
79 dt = t_vals[i + 1] - t_vals[i]
80
81 # Midpoint
82 t_mid = t + dt / 2
83 x_mid = curve_x(t_mid)
84 y_mid = curve_y(t_mid)
85
86 # Tangent vector dr/dt
87 dx_dt = (curve_x(t_mid + 1e-5) - curve_x(t_mid - 1e-5)) / (2e-5)
88 dy_dt = (curve_y(t_mid + 1e-5) - curve_y(t_mid - 1e-5)) / (2e-5)
89
90 # F dot dr
91 F_dot_dr = Fx(x_mid, y_mid) * dx_dt + Fy(x_mid, y_mid) * dy_dt
92 integral += F_dot_dr * dt
93
94 return integral
95
96
97def surface_integral_divergence(Fx, Fy, Fz, x_range, y_range, z_val, n_points=20):
98 """
99 Compute surface integral of F·n over a horizontal surface z=z_val
100 where n = (0, 0, 1) is the upward normal
101 """
102 x_vals = np.linspace(x_range[0], x_range[1], n_points)
103 y_vals = np.linspace(y_range[0], y_range[1], n_points)
104
105 integral = 0
106 dx = (x_range[1] - x_range[0]) / (n_points - 1)
107 dy = (y_range[1] - y_range[0]) / (n_points - 1)
108
109 for x in x_vals:
110 for y in y_vals:
111 # F·n = Fz for upward normal (0,0,1)
112 integral += Fz(x, y, z_val) * dx * dy
113
114 return integral
115
116
117# ============================================================================
118# MAIN DEMONSTRATIONS
119# ============================================================================
120
121print("=" * 70)
122print("VECTOR ANALYSIS - GRADIENT, DIVERGENCE, CURL, AND THEOREMS")
123print("=" * 70)
124
125# Test 1: Gradient of scalar field
126print("\n1. GRADIENT OF SCALAR FIELD")
127print("-" * 70)
128# f(x,y) = x^2 + y^2
129f = lambda x, y: x**2 + y**2
130x, y = 1.0, 2.0
131grad_f = gradient_2d(f, x, y)
132print(f"f(x,y) = x² + y²")
133print(f"∇f at ({x}, {y}) = {grad_f}")
134print(f"Expected: [2x, 2y] = [2, 4]")
135
136# 3D example: f(x,y,z) = x^2 + y^2 + z^2
137f_3d = lambda x, y, z: x**2 + y**2 + z**2
138x, y, z = 1.0, 1.0, 1.0
139grad_f_3d = gradient_3d(f_3d, x, y, z)
140print(f"\nf(x,y,z) = x² + y² + z²")
141print(f"∇f at ({x}, {y}, {z}) = {grad_f_3d}")
142print(f"Expected: [2x, 2y, 2z] = [2, 2, 2]")
143
144# Test 2: Divergence of vector field
145print("\n2. DIVERGENCE OF VECTOR FIELD")
146print("-" * 70)
147# F = (x, y)
148Fx = lambda x, y: x
149Fy = lambda x, y: y
150x, y = 1.0, 2.0
151div_F = divergence_2d(Fx, Fy, x, y)
152print(f"F(x,y) = (x, y)")
153print(f"∇·F at ({x}, {y}) = {div_F:.6f}")
154print(f"Expected: ∂x/∂x + ∂y/∂y = 2")
155
156# 3D example: F = (x, y, z)
157Fx_3d = lambda x, y, z: x
158Fy_3d = lambda x, y, z: y
159Fz_3d = lambda x, y, z: z
160x, y, z = 1.0, 1.0, 1.0
161div_F_3d = divergence_3d(Fx_3d, Fy_3d, Fz_3d, x, y, z)
162print(f"\nF(x,y,z) = (x, y, z)")
163print(f"∇·F at ({x}, {y}, {z}) = {div_F_3d:.6f}")
164print(f"Expected: 3")
165
166# Test 3: Curl of vector field
167print("\n3. CURL OF VECTOR FIELD")
168print("-" * 70)
169# F = (-y, x, 0) - rotation field
170Fx_rot = lambda x, y, z: -y
171Fy_rot = lambda x, y, z: x
172Fz_rot = lambda x, y, z: 0
173x, y, z = 1.0, 0.0, 0.0
174curl_F = curl_3d(Fx_rot, Fy_rot, Fz_rot, x, y, z)
175print(f"F(x,y,z) = (-y, x, 0)")
176print(f"∇×F at ({x}, {y}, {z}) = {curl_F}")
177print(f"Expected: (0, 0, 2)")
178
179# Conservative field: F = ∇(xy)
180Fx_cons = lambda x, y, z: y
181Fy_cons = lambda x, y, z: x
182Fz_cons = lambda x, y, z: 0
183curl_cons = curl_3d(Fx_cons, Fy_cons, Fz_cons, 1.0, 1.0, 0.0)
184print(f"\nConservative field F = ∇(xy) = (y, x, 0)")
185print(f"∇×F = {curl_cons}")
186print(f"Expected: (0, 0, 0) for conservative field")
187
188# Test 4: Line integral
189print("\n4. LINE INTEGRAL")
190print("-" * 70)
191# Integrate F = (y, x) along circle x=cos(t), y=sin(t)
192Fx_field = lambda x, y: y
193Fy_field = lambda x, y: x
194curve_x = lambda t: np.cos(t)
195curve_y = lambda t: np.sin(t)
196t_vals = np.linspace(0, 2*np.pi, 100)
197
198line_int = line_integral_2d(Fx_field, Fy_field, curve_x, curve_y, t_vals)
199print(f"F(x,y) = (y, x)")
200print(f"Curve: circle x=cos(t), y=sin(t), t∈[0,2π]")
201print(f"∮ F·dr = {line_int:.6f}")
202print(f"Expected: 2π = {2*np.pi:.6f}")
203
204# Test 5: Green's theorem
205print("\n5. GREEN'S THEOREM")
206print("-" * 70)
207# ∮ P dx + Q dy = ∬ (∂Q/∂x - ∂P/∂y) dA
208# F = (-y, x), over unit circle
209P = lambda x, y: -y
210Q = lambda x, y: x
211
212# Line integral (already computed above)
213print(f"Field: F = (-y, x)")
214print(f"Region: unit circle")
215
216# Compute ∂Q/∂x - ∂P/∂y
217dQ_dx = lambda x, y: 0 # ∂x/∂x = 1, but we need derivative of Q w.r.t. x
218dP_dy = lambda x, y: 0 # ∂(-y)/∂y = -1
219
220# Actually: ∂Q/∂x - ∂P/∂y = ∂(x)/∂x - ∂(-y)/∂y = 1 - (-1) = 2
221integrand_value = 2
222area = np.pi * 1**2 # unit circle
223double_integral = integrand_value * area
224
225print(f"∮ F·dr (line integral) = {-line_int:.6f}")
226print(f"∬ (∂Q/∂x - ∂P/∂y) dA = 2 × π = {double_integral:.6f}")
227print(f"Green's theorem verified: {np.abs(-line_int - double_integral) < 0.1}")
228
229# Test 6: Divergence theorem
230print("\n6. DIVERGENCE THEOREM")
231print("-" * 70)
232# F = (x, y, z), over cube [0,1]^3
233print(f"Field: F = (x, y, z)")
234print(f"Region: unit cube [0,1]³")
235
236# Surface integral over 6 faces
237flux = 0
238
239# Face x=1 (outward normal = (1,0,0))
240y_vals = np.linspace(0, 1, 10)
241z_vals = np.linspace(0, 1, 10)
242for y in y_vals:
243 for z in z_vals:
244 flux += Fx_3d(1, y, z) * (1/10) * (1/10)
245
246# Face x=0 (outward normal = (-1,0,0))
247for y in y_vals:
248 for z in z_vals:
249 flux -= Fx_3d(0, y, z) * (1/10) * (1/10)
250
251# Face y=1
252for x in np.linspace(0, 1, 10):
253 for z in z_vals:
254 flux += Fy_3d(x, 1, z) * (1/10) * (1/10)
255
256# Face y=0
257for x in np.linspace(0, 1, 10):
258 for z in z_vals:
259 flux -= Fy_3d(x, 0, z) * (1/10) * (1/10)
260
261# Face z=1
262for x in np.linspace(0, 1, 10):
263 for y in y_vals:
264 flux += Fz_3d(x, y, 1) * (1/10) * (1/10)
265
266# Face z=0
267for x in np.linspace(0, 1, 10):
268 for y in y_vals:
269 flux -= Fz_3d(x, y, 0) * (1/10) * (1/10)
270
271print(f"Surface integral ∬ F·n dS = {flux:.6f}")
272
273# Volume integral of divergence
274# div(F) = 3, volume = 1
275volume_integral = 3 * 1
276print(f"Volume integral ∭ ∇·F dV = {volume_integral:.6f}")
277print(f"Divergence theorem verified: {np.abs(flux - volume_integral) < 0.2}")
278
279# Test 7: Stokes' theorem
280print("\n7. STOKES' THEOREM")
281print("-" * 70)
282# F = (-y, x, 0), over unit disk in xy-plane
283print(f"Field: F = (-y, x, 0)")
284print(f"Surface: unit disk in xy-plane")
285print(f"Boundary: unit circle")
286
287# Line integral around boundary (circle)
288Fx_stokes = lambda x, y: -y
289Fy_stokes = lambda x, y: x
290line_int_stokes = line_integral_2d(Fx_stokes, Fy_stokes, curve_x, curve_y, t_vals)
291print(f"∮ F·dr (line integral) = {line_int_stokes:.6f}")
292
293# Surface integral of curl·n
294# curl(F) = (0, 0, 2), n = (0, 0, 1)
295# curl·n = 2
296surface_int_stokes = 2 * np.pi * 1**2
297print(f"∬ (∇×F)·n dS = 2 × π = {surface_int_stokes:.6f}")
298print(f"Stokes' theorem verified: {np.abs(line_int_stokes - surface_int_stokes) < 0.1}")
299
300# Visualization
301if HAS_MATPLOTLIB:
302 print("\n" + "=" * 70)
303 print("VISUALIZATIONS")
304 print("=" * 70)
305
306 fig = plt.figure(figsize=(15, 10))
307
308 # Plot 1: Gradient field
309 ax1 = plt.subplot(2, 3, 1)
310 x = np.linspace(-2, 2, 20)
311 y = np.linspace(-2, 2, 20)
312 X, Y = np.meshgrid(x, y)
313
314 f_scalar = X**2 + Y**2
315 U = 2*X # ∂f/∂x
316 V = 2*Y # ∂f/∂y
317
318 contour = ax1.contour(X, Y, f_scalar, levels=10, cmap='viridis', alpha=0.3)
319 ax1.quiver(X, Y, U, V, color='red', alpha=0.6)
320 ax1.set_xlabel('x')
321 ax1.set_ylabel('y')
322 ax1.set_title('Gradient Field ∇f (f=x²+y²)')
323 ax1.set_aspect('equal')
324 ax1.grid(True, alpha=0.3)
325
326 # Plot 2: Divergence visualization
327 ax2 = plt.subplot(2, 3, 2)
328 x = np.linspace(-2, 2, 20)
329 y = np.linspace(-2, 2, 20)
330 X, Y = np.meshgrid(x, y)
331 U = X
332 V = Y
333
334 div = np.ones_like(X) * 2 # div(F) = 2
335 im = ax2.contourf(X, Y, div, levels=10, cmap='RdBu_r', alpha=0.6)
336 ax2.quiver(X, Y, U, V, color='black', alpha=0.6)
337 ax2.set_xlabel('x')
338 ax2.set_ylabel('y')
339 ax2.set_title('Divergence: F=(x,y), ∇·F=2')
340 ax2.set_aspect('equal')
341 plt.colorbar(im, ax=ax2)
342
343 # Plot 3: Curl visualization
344 ax3 = plt.subplot(2, 3, 3)
345 x = np.linspace(-2, 2, 20)
346 y = np.linspace(-2, 2, 20)
347 X, Y = np.meshgrid(x, y)
348 U = -Y
349 V = X
350
351 ax3.quiver(X, Y, U, V, color='blue', alpha=0.6)
352
353 # Add circular streamlines
354 theta = np.linspace(0, 2*np.pi, 100)
355 for r in [0.5, 1.0, 1.5]:
356 ax3.plot(r*np.cos(theta), r*np.sin(theta), 'r--', alpha=0.3, linewidth=0.5)
357
358 ax3.set_xlabel('x')
359 ax3.set_ylabel('y')
360 ax3.set_title('Curl: F=(-y,x,0), ∇×F=(0,0,2)')
361 ax3.set_aspect('equal')
362 ax3.grid(True, alpha=0.3)
363
364 # Plot 4: Line integral path
365 ax4 = plt.subplot(2, 3, 4)
366 t = np.linspace(0, 2*np.pi, 100)
367 x_circle = np.cos(t)
368 y_circle = np.sin(t)
369
370 # Draw vector field
371 x_grid = np.linspace(-1.5, 1.5, 15)
372 y_grid = np.linspace(-1.5, 1.5, 15)
373 X_grid, Y_grid = np.meshgrid(x_grid, y_grid)
374 U_grid = Y_grid
375 V_grid = X_grid
376
377 ax4.quiver(X_grid, Y_grid, U_grid, V_grid, alpha=0.4)
378 ax4.plot(x_circle, y_circle, 'r-', linewidth=3, label='Integration path')
379 ax4.arrow(1, 0, 0, 0.3, head_width=0.1, head_length=0.1, fc='red', ec='red')
380 ax4.set_xlabel('x')
381 ax4.set_ylabel('y')
382 ax4.set_title('Line Integral: F=(y,x)')
383 ax4.legend()
384 ax4.set_aspect('equal')
385 ax4.grid(True, alpha=0.3)
386
387 # Plot 5: Green's theorem
388 ax5 = plt.subplot(2, 3, 5)
389 x = np.linspace(-1.5, 1.5, 20)
390 y = np.linspace(-1.5, 1.5, 20)
391 X, Y = np.meshgrid(x, y)
392
393 # ∂Q/∂x - ∂P/∂y = 2 for F=(-y,x)
394 integrand = np.ones_like(X) * 2
395
396 # Mask to show only inside circle
397 mask = X**2 + Y**2 <= 1
398 integrand_masked = np.where(mask, integrand, np.nan)
399
400 im = ax5.contourf(X, Y, integrand_masked, levels=10, cmap='YlOrRd', alpha=0.7)
401 ax5.plot(x_circle, y_circle, 'b-', linewidth=3, label='∂D')
402 ax5.set_xlabel('x')
403 ax5.set_ylabel('y')
404 ax5.set_title("Green's Theorem: ∂Q/∂x - ∂P/∂y = 2")
405 ax5.legend()
406 ax5.set_aspect('equal')
407 ax5.grid(True, alpha=0.3)
408
409 # Plot 6: Conservative vs non-conservative field
410 ax6 = plt.subplot(2, 3, 6)
411 x = np.linspace(-2, 2, 15)
412 y = np.linspace(-2, 2, 15)
413 X, Y = np.meshgrid(x, y)
414
415 # Conservative: F = ∇(xy) = (y, x)
416 U_cons = Y
417 V_cons = X
418
419 ax6.quiver(X, Y, U_cons, V_cons, color='green', alpha=0.6, label='Conservative')
420
421 # Draw potential lines
422 for c in [-2, -1, 0, 1, 2]:
423 x_pot = np.linspace(-2, 2, 100)
424 y_pot = c / x_pot
425 y_pot = np.clip(y_pot, -2, 2)
426 mask = np.abs(x_pot) > 0.1
427 ax6.plot(x_pot[mask], y_pot[mask], 'g--', alpha=0.3, linewidth=0.5)
428
429 ax6.set_xlabel('x')
430 ax6.set_ylabel('y')
431 ax6.set_title('Conservative Field: F=∇(xy)=(y,x)')
432 ax6.set_aspect('equal')
433 ax6.grid(True, alpha=0.3)
434
435 plt.tight_layout()
436 plt.savefig('/opt/projects/01_Personal/03_Study/examples/Mathematical_Methods/05_vector_analysis.png', dpi=150)
437 print("Saved visualization: 05_vector_analysis.png")
438 plt.close()
439
440print("\n" + "=" * 70)
441print("DEMONSTRATION COMPLETE")
442print("=" * 70)