1"""
2Tensors - Operations, Einstein Summation, and Coordinate Transformations
3
4This script demonstrates:
5- Tensor operations with numpy
6- Index notation and Einstein summation (np.einsum)
7- Metric tensor
8- Christoffel symbols
9- Coordinate transformations
10- Tensor contraction
11- Tensor products
12"""
13
14import numpy as np
15
16try:
17 import matplotlib.pyplot as plt
18 HAS_MATPLOTLIB = True
19except ImportError:
20 HAS_MATPLOTLIB = False
21 print("Warning: matplotlib not available, skipping visualizations")
22
23
24def kronecker_delta(n):
25 """Kronecker delta tensor δ_ij"""
26 return np.eye(n)
27
28
29def levi_civita_3d():
30 """
31 Levi-Civita (permutation) tensor ε_ijk in 3D
32 ε_ijk = +1 for even permutations, -1 for odd, 0 otherwise
33 """
34 epsilon = np.zeros((3, 3, 3))
35 epsilon[0, 1, 2] = epsilon[1, 2, 0] = epsilon[2, 0, 1] = 1
36 epsilon[0, 2, 1] = epsilon[2, 1, 0] = epsilon[1, 0, 2] = -1
37 return epsilon
38
39
40def metric_tensor_euclidean_2d():
41 """Euclidean metric in 2D: g_ij = δ_ij"""
42 return np.eye(2)
43
44
45def metric_tensor_polar():
46 """
47 Metric tensor in polar coordinates (r, θ)
48 ds² = dr² + r²dθ²
49 g_ij = [[1, 0], [0, r²]]
50 """
51 def g(r):
52 return np.array([[1, 0], [0, r**2]])
53 return g
54
55
56def metric_tensor_spherical():
57 """
58 Metric tensor in spherical coordinates (r, θ, φ)
59 ds² = dr² + r²dθ² + r²sin²(θ)dφ²
60 """
61 def g(r, theta):
62 return np.array([
63 [1, 0, 0],
64 [0, r**2, 0],
65 [0, 0, r**2 * np.sin(theta)**2]
66 ])
67 return g
68
69
70def christoffel_symbols_polar(r):
71 """
72 Compute Christoffel symbols Γ^k_ij for polar coordinates
73 Γ^r_θθ = -r
74 Γ^θ_rθ = Γ^θ_θr = 1/r
75 """
76 Gamma = np.zeros((2, 2, 2))
77
78 # Γ^r_θθ = -r (component [0, 1, 1])
79 Gamma[0, 1, 1] = -r
80
81 # Γ^θ_rθ = Γ^θ_θr = 1/r (components [1, 0, 1] and [1, 1, 0])
82 Gamma[1, 0, 1] = 1 / r
83 Gamma[1, 1, 0] = 1 / r
84
85 return Gamma
86
87
88def christoffel_symbols_spherical(r, theta):
89 """
90 Compute Christoffel symbols Γ^k_ij for spherical coordinates
91 """
92 Gamma = np.zeros((3, 3, 3))
93
94 # Non-zero components
95 # Γ^r_θθ = -r
96 Gamma[0, 1, 1] = -r
97
98 # Γ^r_φφ = -r sin²(θ)
99 Gamma[0, 2, 2] = -r * np.sin(theta)**2
100
101 # Γ^θ_rθ = Γ^θ_θr = 1/r
102 Gamma[1, 0, 1] = 1 / r
103 Gamma[1, 1, 0] = 1 / r
104
105 # Γ^θ_φφ = -sin(θ)cos(θ)
106 Gamma[1, 2, 2] = -np.sin(theta) * np.cos(theta)
107
108 # Γ^φ_rφ = Γ^φ_φr = 1/r
109 Gamma[2, 0, 2] = 1 / r
110 Gamma[2, 2, 0] = 1 / r
111
112 # Γ^φ_θφ = Γ^φ_φθ = cot(θ)
113 if np.sin(theta) != 0:
114 Gamma[2, 1, 2] = np.cos(theta) / np.sin(theta)
115 Gamma[2, 2, 1] = np.cos(theta) / np.sin(theta)
116
117 return Gamma
118
119
120def coordinate_transform_cartesian_to_polar(x, y):
121 """Transform from Cartesian to polar coordinates"""
122 r = np.sqrt(x**2 + y**2)
123 theta = np.arctan2(y, x)
124 return r, theta
125
126
127def coordinate_transform_polar_to_cartesian(r, theta):
128 """Transform from polar to Cartesian coordinates"""
129 x = r * np.cos(theta)
130 y = r * np.sin(theta)
131 return x, y
132
133
134def jacobian_cartesian_to_polar(x, y):
135 """
136 Jacobian matrix for Cartesian to polar transformation
137 J_ij = ∂x^i_new / ∂x^j_old
138 """
139 r = np.sqrt(x**2 + y**2)
140 if r == 0:
141 return np.eye(2)
142
143 # dr/dx, dr/dy
144 # dθ/dx, dθ/dy
145 J = np.array([
146 [x/r, y/r],
147 [-y/(r**2), x/(r**2)]
148 ])
149 return J
150
151
152def transform_vector_contravariant(v, J):
153 """
154 Transform contravariant vector: v'^i = J^i_j v^j
155 """
156 return J @ v
157
158
159def transform_covector_covariant(w, J):
160 """
161 Transform covariant vector (covector): w'_i = (J^T)^(-1)_ij w_j
162 """
163 return np.linalg.inv(J.T) @ w
164
165
166# ============================================================================
167# MAIN DEMONSTRATIONS
168# ============================================================================
169
170print("=" * 70)
171print("TENSORS - OPERATIONS, EINSTEIN SUMMATION, TRANSFORMATIONS")
172print("=" * 70)
173
174# Test 1: Basic tensor operations
175print("\n1. BASIC TENSOR OPERATIONS")
176print("-" * 70)
177
178# Rank-2 tensors (matrices)
179A = np.array([[1, 2], [3, 4]])
180B = np.array([[5, 6], [7, 8]])
181
182print("Tensor A:")
183print(A)
184print("\nTensor B:")
185print(B)
186
187# Tensor addition
188print("\nA + B:")
189print(A + B)
190
191# Tensor product (outer product)
192C = np.outer(A.flatten(), B.flatten()).reshape(2, 2, 2, 2)
193print(f"\nTensor product A ⊗ B shape: {C.shape}")
194
195# Trace (contraction)
196trace_A = np.trace(A)
197print(f"\nTrace of A (contraction): {trace_A}")
198
199# Test 2: Einstein summation
200print("\n2. EINSTEIN SUMMATION CONVENTION (np.einsum)")
201print("-" * 70)
202
203# Matrix multiplication: C_ij = A_ik B_kj
204A = np.array([[1, 2], [3, 4]])
205B = np.array([[5, 6], [7, 8]])
206
207C_matmul = np.einsum('ik,kj->ij', A, B)
208C_direct = A @ B
209
210print("Matrix multiplication: C = A @ B")
211print("Using einsum:")
212print(C_matmul)
213print("Direct computation:")
214print(C_direct)
215print(f"Match: {np.allclose(C_matmul, C_direct)}")
216
217# Trace: tr(A) = A_ii
218trace_einsum = np.einsum('ii->', A)
219print(f"\nTrace using einsum: {trace_einsum}")
220print(f"Trace using np.trace: {np.trace(A)}")
221
222# Vector dot product: c = a_i b_i
223a = np.array([1, 2, 3])
224b = np.array([4, 5, 6])
225
226dot_einsum = np.einsum('i,i->', a, b)
227dot_direct = np.dot(a, b)
228
229print(f"\nDot product a·b:")
230print(f" Using einsum: {dot_einsum}")
231print(f" Using np.dot: {dot_direct}")
232
233# Outer product: C_ij = a_i b_j
234outer_einsum = np.einsum('i,j->ij', a, b)
235outer_direct = np.outer(a, b)
236
237print(f"\nOuter product a ⊗ b:")
238print("Using einsum:")
239print(outer_einsum)
240print("Using np.outer:")
241print(outer_direct)
242
243# Batch matrix multiplication: D_ijk = A_ij B_jk for multiple matrices
244n_batch = 3
245A_batch = np.random.rand(n_batch, 2, 3)
246B_batch = np.random.rand(n_batch, 3, 2)
247
248D_batch = np.einsum('nij,njk->nik', A_batch, B_batch)
249print(f"\nBatch matrix multiplication shape: {D_batch.shape}")
250
251# Test 3: Metric tensor
252print("\n3. METRIC TENSOR")
253print("-" * 70)
254
255# Euclidean 2D
256g_euclidean = metric_tensor_euclidean_2d()
257print("Euclidean metric (Cartesian):")
258print(g_euclidean)
259
260# Polar coordinates
261r = 2.0
262g_polar_func = metric_tensor_polar()
263g_polar = g_polar_func(r)
264print(f"\nPolar metric at r={r}:")
265print(g_polar)
266
267# Line element ds²
268dr = 0.1
269dtheta = 0.05
270ds_squared = g_polar[0, 0] * dr**2 + g_polar[1, 1] * dtheta**2
271print(f"\nLine element ds² = dr² + r²dθ²")
272print(f"For dr={dr}, dθ={dtheta}, r={r}:")
273print(f"ds² = {ds_squared:.6f}")
274print(f"ds = {np.sqrt(ds_squared):.6f}")
275
276# Test 4: Christoffel symbols
277print("\n4. CHRISTOFFEL SYMBOLS")
278print("-" * 70)
279
280r = 2.0
281Gamma_polar = christoffel_symbols_polar(r)
282
283print(f"Polar coordinates at r={r}:")
284print(f"Γ^r_θθ = {Gamma_polar[0, 1, 1]:.4f} (expected: {-r})")
285print(f"Γ^θ_rθ = {Gamma_polar[1, 0, 1]:.4f} (expected: {1/r:.4f})")
286print(f"Γ^θ_θr = {Gamma_polar[1, 1, 0]:.4f} (expected: {1/r:.4f})")
287
288# Spherical coordinates
289r, theta = 3.0, np.pi / 4
290Gamma_spherical = christoffel_symbols_spherical(r, theta)
291
292print(f"\nSpherical coordinates at r={r}, θ={np.degrees(theta):.1f}°:")
293print(f"Γ^r_θθ = {Gamma_spherical[0, 1, 1]:.4f} (expected: {-r})")
294print(f"Γ^r_φφ = {Gamma_spherical[0, 2, 2]:.4f} (expected: {-r*np.sin(theta)**2:.4f})")
295print(f"Γ^θ_rθ = {Gamma_spherical[1, 0, 1]:.4f} (expected: {1/r:.4f})")
296
297# Test 5: Coordinate transformations
298print("\n5. COORDINATE TRANSFORMATIONS")
299print("-" * 70)
300
301# Cartesian to polar
302x, y = 3.0, 4.0
303r, theta = coordinate_transform_cartesian_to_polar(x, y)
304
305print(f"Cartesian: (x, y) = ({x}, {y})")
306print(f"Polar: (r, θ) = ({r:.4f}, {np.degrees(theta):.2f}°)")
307
308# Back to Cartesian
309x_back, y_back = coordinate_transform_polar_to_cartesian(r, theta)
310print(f"Back to Cartesian: ({x_back:.4f}, {y_back:.4f})")
311
312# Jacobian
313J = jacobian_cartesian_to_polar(x, y)
314print(f"\nJacobian matrix (∂polar/∂cartesian):")
315print(J)
316
317# Verify det(J) = 1/r for area element transformation
318det_J = np.linalg.det(J)
319print(f"det(J) = {det_J:.6f}")
320print(f"Expected 1/r = {1/r:.6f}")
321
322# Test 6: Vector transformation
323print("\n6. VECTOR TRANSFORMATION")
324print("-" * 70)
325
326# Velocity vector in Cartesian coordinates
327v_cartesian = np.array([1.0, 2.0])
328print(f"Velocity in Cartesian: v = {v_cartesian}")
329
330# Transform to polar
331x, y = 3.0, 4.0
332J = jacobian_cartesian_to_polar(x, y)
333v_polar = transform_vector_contravariant(v_cartesian, J)
334
335print(f"Velocity in polar: v = {v_polar}")
336print(f" v^r = {v_polar[0]:.4f}")
337print(f" v^θ = {v_polar[1]:.4f}")
338
339# Transform back
340J_inverse = np.linalg.inv(J)
341v_back = transform_vector_contravariant(v_polar, J_inverse)
342print(f"Back to Cartesian: v = {v_back}")
343print(f"Match: {np.allclose(v_cartesian, v_back)}")
344
345# Test 7: Levi-Civita tensor and cross product
346print("\n7. LEVI-CIVITA TENSOR AND CROSS PRODUCT")
347print("-" * 70)
348
349epsilon = levi_civita_3d()
350print("Levi-Civita tensor ε_ijk (3D)")
351
352# Cross product using Einstein summation
353# (a × b)_i = ε_ijk a_j b_k
354a = np.array([1, 0, 0])
355b = np.array([0, 1, 0])
356
357cross_einsum = np.einsum('ijk,j,k->i', epsilon, a, b)
358cross_direct = np.cross(a, b)
359
360print(f"\na = {a}")
361print(f"b = {b}")
362print(f"a × b (einsum): {cross_einsum}")
363print(f"a × b (np.cross): {cross_direct}")
364
365# Verify ε_ijk ε_imn = δ_jm δ_kn - δ_jn δ_km
366product = np.einsum('ijk,imn->jkmn', epsilon, epsilon)
367delta = kronecker_delta(3)
368
369expected = (np.einsum('jm,kn->jkmn', delta, delta) -
370 np.einsum('jn,km->jkmn', delta, delta))
371
372print(f"\nVerify ε_ijk ε_imn = δ_jm δ_kn - δ_jn δ_km:")
373print(f"Match: {np.allclose(product, expected)}")
374
375# Test 8: Raising and lowering indices
376print("\n8. RAISING AND LOWERING INDICES")
377print("-" * 70)
378
379# In Euclidean space, g_ij = g^ij = δ_ij
380g = metric_tensor_euclidean_2d()
381g_inv = np.linalg.inv(g)
382
383print("Metric tensor g_ij:")
384print(g)
385print("\nInverse metric g^ij:")
386print(g_inv)
387
388# Covariant vector (covector)
389w_down = np.array([1, 2])
390print(f"\nCovariant vector w_i = {w_down}")
391
392# Raise index: w^i = g^ij w_j
393w_up = g_inv @ w_down
394print(f"Contravariant vector w^i = g^ij w_j = {w_up}")
395
396# In Euclidean space, they're the same
397print(f"Match (Euclidean): {np.allclose(w_down, w_up)}")
398
399# In polar coordinates
400r = 2.0
401g_polar_func = metric_tensor_polar()
402g_polar = g_polar_func(r)
403g_polar_inv = np.linalg.inv(g_polar)
404
405w_down_polar = np.array([1, 2])
406w_up_polar = g_polar_inv @ w_down_polar
407
408print(f"\nIn polar coordinates at r={r}:")
409print(f"Covariant w_i = {w_down_polar}")
410print(f"Contravariant w^i = {w_up_polar}")
411
412# Visualization
413if HAS_MATPLOTLIB:
414 print("\n" + "=" * 70)
415 print("VISUALIZATIONS")
416 print("=" * 70)
417
418 fig = plt.figure(figsize=(15, 10))
419
420 # Plot 1: Coordinate grid transformation
421 ax1 = plt.subplot(2, 3, 1)
422
423 # Cartesian grid
424 x = np.linspace(-3, 3, 7)
425 y = np.linspace(-3, 3, 7)
426
427 for xi in x:
428 ax1.plot([xi, xi], [-3, 3], 'b-', alpha=0.5, linewidth=0.8)
429 for yi in y:
430 ax1.plot([-3, 3], [yi, yi], 'r-', alpha=0.5, linewidth=0.8)
431
432 ax1.set_xlabel('x')
433 ax1.set_ylabel('y')
434 ax1.set_title('Cartesian Grid')
435 ax1.set_aspect('equal')
436 ax1.grid(True, alpha=0.3)
437
438 # Plot 2: Polar grid
439 ax2 = plt.subplot(2, 3, 2)
440
441 r_vals = np.linspace(0.5, 3, 6)
442 theta_vals = np.linspace(0, 2*np.pi, 13)
443
444 # Radial lines
445 for theta in theta_vals:
446 r_line = np.linspace(0, 3, 50)
447 x_line = r_line * np.cos(theta)
448 y_line = r_line * np.sin(theta)
449 ax2.plot(x_line, y_line, 'b-', alpha=0.5, linewidth=0.8)
450
451 # Circular lines
452 for r in r_vals:
453 theta_circle = np.linspace(0, 2*np.pi, 100)
454 x_circle = r * np.cos(theta_circle)
455 y_circle = r * np.sin(theta_circle)
456 ax2.plot(x_circle, y_circle, 'r-', alpha=0.5, linewidth=0.8)
457
458 ax2.set_xlabel('x')
459 ax2.set_ylabel('y')
460 ax2.set_title('Polar Grid in Cartesian Space')
461 ax2.set_aspect('equal')
462 ax2.grid(True, alpha=0.3)
463
464 # Plot 3: Vector transformation
465 ax3 = plt.subplot(2, 3, 3)
466
467 x0, y0 = 2.0, 1.5
468 v_cart = np.array([1.0, 0.5])
469
470 # Original vector
471 ax3.arrow(x0, y0, v_cart[0], v_cart[1], head_width=0.15, head_length=0.1,
472 fc='blue', ec='blue', linewidth=2, label='Cartesian')
473
474 # Transform to polar basis
475 J = jacobian_cartesian_to_polar(x0, y0)
476 v_polar = transform_vector_contravariant(v_cart, J)
477
478 # Polar basis vectors at (x0, y0)
479 r0, theta0 = coordinate_transform_cartesian_to_polar(x0, y0)
480
481 # e_r in Cartesian coordinates
482 e_r = np.array([np.cos(theta0), np.sin(theta0)])
483 # e_θ in Cartesian coordinates
484 e_theta = np.array([-np.sin(theta0), np.cos(theta0)])
485
486 # Reconstruct vector in Cartesian using polar components
487 v_reconstructed = v_polar[0] * e_r + v_polar[1] * e_theta
488
489 ax3.arrow(x0, y0, v_reconstructed[0], v_reconstructed[1],
490 head_width=0.15, head_length=0.1, fc='red', ec='red',
491 linewidth=2, alpha=0.5, linestyle='--', label='From polar')
492
493 # Show basis vectors
494 ax3.arrow(x0, y0, e_r[0], e_r[1], head_width=0.1, head_length=0.05,
495 fc='green', ec='green', alpha=0.5, label='e_r')
496 ax3.arrow(x0, y0, e_theta[0], e_theta[1], head_width=0.1, head_length=0.05,
497 fc='orange', ec='orange', alpha=0.5, label='e_θ')
498
499 ax3.plot(x0, y0, 'ko', markersize=8)
500 ax3.set_xlabel('x')
501 ax3.set_ylabel('y')
502 ax3.set_title('Vector Transformation')
503 ax3.legend()
504 ax3.grid(True, alpha=0.3)
505 ax3.set_aspect('equal')
506 ax3.set_xlim(0, 5)
507 ax3.set_ylim(0, 4)
508
509 # Plot 4: Metric tensor components in polar
510 ax4 = plt.subplot(2, 3, 4)
511
512 r_range = np.linspace(0.5, 5, 100)
513 g_rr = np.ones_like(r_range)
514 g_theta_theta = r_range**2
515
516 ax4.plot(r_range, g_rr, 'b-', linewidth=2, label='g_rr = 1')
517 ax4.plot(r_range, g_theta_theta, 'r-', linewidth=2, label='g_θθ = r²')
518 ax4.set_xlabel('r')
519 ax4.set_ylabel('Metric component')
520 ax4.set_title('Polar Metric Tensor Components')
521 ax4.legend()
522 ax4.grid(True, alpha=0.3)
523
524 # Plot 5: Christoffel symbols
525 ax5 = plt.subplot(2, 3, 5)
526
527 r_range = np.linspace(0.5, 5, 100)
528 Gamma_r_theta_theta = -r_range
529 Gamma_theta_r_theta = 1 / r_range
530
531 ax5.plot(r_range, Gamma_r_theta_theta, 'b-', linewidth=2, label='Γ^r_θθ = -r')
532 ax5.plot(r_range, Gamma_theta_r_theta, 'r-', linewidth=2, label='Γ^θ_rθ = 1/r')
533 ax5.set_xlabel('r')
534 ax5.set_ylabel('Christoffel symbol')
535 ax5.set_title('Polar Christoffel Symbols')
536 ax5.legend()
537 ax5.grid(True, alpha=0.3)
538
539 # Plot 6: Cross product with Levi-Civita
540 ax6 = plt.subplot(2, 3, 6, projection='3d')
541
542 # Vectors
543 a = np.array([1, 0, 0])
544 b = np.array([0, 1, 0])
545 c = np.cross(a, b)
546
547 origin = np.array([0, 0, 0])
548
549 ax6.quiver(origin[0], origin[1], origin[2], a[0], a[1], a[2],
550 color='red', arrow_length_ratio=0.2, linewidth=2, label='a')
551 ax6.quiver(origin[0], origin[1], origin[2], b[0], b[1], b[2],
552 color='blue', arrow_length_ratio=0.2, linewidth=2, label='b')
553 ax6.quiver(origin[0], origin[1], origin[2], c[0], c[1], c[2],
554 color='green', arrow_length_ratio=0.2, linewidth=2, label='a×b')
555
556 ax6.set_xlabel('x')
557 ax6.set_ylabel('y')
558 ax6.set_zlabel('z')
559 ax6.set_title('Cross Product: a × b')
560 ax6.legend()
561 ax6.set_xlim(-0.5, 1.5)
562 ax6.set_ylim(-0.5, 1.5)
563 ax6.set_zlim(-0.5, 1.5)
564
565 plt.tight_layout()
566 plt.savefig('/opt/projects/01_Personal/03_Study/examples/Mathematical_Methods/15_tensors.png', dpi=150)
567 print("Saved visualization: 15_tensors.png")
568 plt.close()
569
570print("\n" + "=" * 70)
571print("DEMONSTRATION COMPLETE")
572print("=" * 70)