# References

1. Joseph C. Kolecki, (September 2002) An Introduction to Tensors for Students of Physics and Engineering National Aeronautics and Space Administration, Glenn Research Center, Cleveland, Ohio 44135
2. Daniel A Fleisch - A student's guide to vectors and tensors (2012, Cambridge University Press)
3. [H.K._Dass]_Advanced_Engineering_Mathematics(z-lib.org) (21st Revised Edition)
4. K.A. Stroud, Dexter J. Booth, Advanced Engineering Mathematics

## Imports


```python
import numpy as np
import math
import sympy
import sympy as sp
from sympy import Eq, IndexedBase, symbols, Idx, Indexed, Sum, S, N
from sympy.functions.special.tensor_functions import KroneckerDelta
from sympy.vector import Vector, CoordSys3D, AxisOrienter, BodyOrienter, Del, curl, divergence, gradient, is_conservative, is_solenoidal, scalar_potential, Point, scalar_potential_difference, Del, express, matrix_to_vector
import matplotlib.pyplot as plt
from sympy.physics.vector import ReferenceFrame
import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
from IPython.display import display
from IPython.display import display_latex
from sympy import latex
from sympy import Array, Matrix, transpose, zeros, diff, Function, Derivative, cos, sin, sqrt, solve, linsolve, acos, atan, asin

from sympy import symbols
from sympy.plotting import plot
```

## Setup


```python
sympy.init_printing()
np.set_printoptions(precision=3)
```


```python
a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x, y, z = sympy.symbols('a:z')
a_1, a_2, a_3, b_1, b_2, b_3, c_1, c_2, c_3, d_1, d_2, d_3 = symbols('a_1 a_2 a_3 b_1 b_2 b_3 c_1 c_2 c_3 d_1 d_2 d_3')
u_1, u_2, u_3, v_1, v_2, v_3, h_1, h_2, h_3 = symbols('u_1 u_2 u_3 v_1 v_2 v_3 h_1 h_2 h_3')
alpha, beta, gamma, theta, phi, rho, omega = sympy.symbols('alpha beta gamma theta phi rho omega')
Alpha, Beta, Gamma, Theta, Phi, Rho, Omega = sympy.symbols('Alpha Beta Gamma Theta Phi Rho Omega')
```

## $LaTeX$ Cheat Sheet

### Symbols

$\mathbb{R}$  
$\bar{R}$  
$\frac{\partial^2 R}{\partial t}$  
$\overline{R}$  
$\vec{R}$  
$\hat{x}$  
$\widehat{x}$  
$\therefore$  
$\because$  
$\Rightarrow$  
$\rightarrow$  
$\implies$  
$\iff$  
$\leftrightarrow$  
$\Leftrightarrow$  
$\equiv$  
$\subset$  
$\supset$  
$\wedge$  
$\parallel$  
$\oplus$  
$\ominus$  
$\otimes$  
$\times$  
$\top$  
$\bot$  
$\forall$  
$\exists$  
$\exists!$  

## Functions


```python
# Usage: display_equation('u_x', x)
def display_equation(idx, symObj):
    if(isinstance(idx, str)):
        eqn = '\\[' + idx + ' = ' + latex(symObj) + '\\]'
        display_latex(eqn, raw=True)
    else:
        eqn = '\\[' + latex(idx) + ' = ' + latex(symObj) + '\\]'
        display_latex(eqn, raw=True)
    return
```


```python
# Usage: display_full_latex('u_x')
def display_full_latex(idx):
    if(isinstance(idx, str)):
        eqn = '\\[' + idx + '\\]'
        display_latex(eqn, raw=True)
    else:
        eqn = '\\[' + latex(idx) + '\\]'
        display_latex(eqn, raw=True)
    return
```


```python
def vplot2d(O, a, origin_label='O', tip_label='P'):
    dx = a[0][0] - O[0][0]
    dy = a[1][0] - O[1][0]

    head_length = 0.3
    vec_ab = [dx,dy]

    vec_ab_magnitude = math.sqrt(dx**2+dy**2)

    dx = dx / vec_ab_magnitude
    dy = dy / vec_ab_magnitude

    vec_ab_magnitude = vec_ab_magnitude - head_length

    ax = plt.axes()
    ax.set_aspect('equal', 'box')
    #ax.axis([-range, range, -range, range]) Sub for xlim and ylim

    ax.arrow(O[0][0], O[1][0], vec_ab_magnitude*dx, vec_ab_magnitude*dy, head_width=0.1, head_length=head_length, fc='lightblue', ec='black')
    plt.scatter(O[0][0],O[1][0],color='black')
    plt.scatter(a[0][0],a[1][0],color='black')

    ax.annotate(origin_label, (O[0][0]-0.4,O[1][0]),fontsize=14)
    ax.annotate(tip_label, (a[0][0]+0.3,a[1][0]),fontsize=14)
    return
```


```python
def rotation(theta):
    R = zeros(3, 3)
    R[0,0] = cos(theta)
    R[0,1] = -sin(theta)
    R[0,2] = 0
    R[1,0] = sin(theta)
    R[1,1] = cos(theta)
    R[1,2] = 0
    R[2,0] = 0
    R[2,1] = 0
    R[2,2] = 1
    return R
    return R
```


```python
def inRadians(theta):
    return (sympy.pi / 180) * theta

def inDegrees(theta):
    return theta * (180 / sympy.pi)
```

### Coordinate Transformation


```python
x_prime, y_prime, z_prime = sympy.symbols('x\' y\' z\'')
l_1, l_2, l_3, m_1, m_2, m_3, n_1, n_2, n_3 = sympy.symbols('l_1 l_2 l_3 m_1, m_2 m_3 n_1 n_2 n_3')
x, y, z = sympy.symbols('x y z')
```


```python
x_prime, y_prime, z_prime = sp.symbols('x\' y\' z\'')
l_1, l_2, l_3, m_1, m_2, m_3, n_1, n_2, n_3 = sp.symbols('l_1 l_2 l_3 m_1, m_2 m_3 n_1 n_2 n_3')
x, y, z = sp.symbols('x y z')

# Equation 1
x_prime = l_1 * x + m_1 * y + n_1 * z
y_prime = l_2 * x + m_2 * y + n_2 * z
z_prime = l_3 * x + m_3 * y + n_3 * z

display_equation('x\'', x_prime)
display_equation('y\'', y_prime)
display_equation('z\'', z_prime)
```


\[x' = l_{1} x + m_{1} y + n_{1} z\]



\[y' = l_{2} x + m_{2} y + n_{2} z\]



\[z' = l_{3} x + m_{3} y + n_{3} z\]



```python
# Redefining the primes so equation 1 doesn't get substituted in
x_prime, y_prime, z_prime = sp.symbols('x\' y\' z\'', Integer=True)

x = l_1 * x_prime + m_1 * y_prime + n_1 * z_prime
y = l_2 * x_prime + m_2 * y_prime + n_2 * z_prime
z = l_3 * x_prime + m_3 * y_prime + n_3 * z_prime

display_equation('x\'', x)
display_equation('y\'', y)
display_equation('z\'', z)
```


\[x' = l_{1} x' + m_{1} y' + n_{1} z'\]



\[y' = l_{2} x' + m_{2} y' + n_{2} z'\]



\[z' = l_{3} x' + m_{3} y' + n_{3} z'\]


|  | $x$ | $y$ | $z$ |
|---|---|---|---|
| $x'$ | $l_1$ | $m_1$ | $n_1$ |
| $y'$ | $l_2$ | $m_2$ | $n_2$ |
| $z'$ | $l_3$ | $m_3$ | $n_3$ |

### OR

$\bar{x_1} = l_{11}x_1 + l_{21}x_2 + l_{31}x_3$  
$\bar{x_2} = l_{13}x_1 + l_{22}x_2 + l_{32}x_3$  
$\bar{x_3} = l_{13}x_1 + l_{23}x_2 + l_{33}x_3$  

$x_1 = l_{11}\bar{x_1} + l_{12}\bar{x_2} + l_{13}\bar{x_3}$  
$x_2 = l_{21}\bar{x_1} + l_{22}\bar{x_2} + l_{23}\bar{x_3}$  
$x_3 = l_{31}\bar{x_1} + l_{32}\bar{x_2} + l_{33}\bar{x_3}$  

### With column by row indexing

|  | $x_1$ | $x_2$ | $x_3$ |
|---|---|---|---|
| $\bar{x_1}$ | $l_{11}$ | $l_{21}$ | $l_{31}$ |
| $\bar{x_2}$ | $l_{12}$ | $l_{22}$ | $l_{32}$ |
| $\bar{x_3}$ | $l_{13}$ | $l_{23}$ | $l_{33}$ |

$\bar{x_j} = l_{ji}x_i$  
$x_i = l_{ji}\bar{x_j}$  
Where  
$i$ = column  
$j$ = row 

#### With row by column indexing (More common matrix indexing, used here)

|  | $\bar{x_1}$ | $\bar{x_2}$ | $\bar{x_3}$ |
|---|---|---|---|
| $x_1$ | $l_{11}$ | $l_{12}$ | $l_{13}$ |
| $x_2$ | $l_{21}$ | $l_{22}$ | $l_{23}$ |
| $x_3$ | $l_{31}$ | $l_{32}$ | $l_{33}$ |

$x_i = l_{ij}\bar{x_j}$  
$\bar{x_j} = l_{ij}x_i$  
Where  
$i$ = row  
$j$ = column 

#### Relation between the Direction Cosines of Three Mutually Perpendicular Straight Lines

1. The dot product between any two is 0.
2. The dot product of any vector on itself is 1.

In Einstein summation notation


```python
i, k = sympy.symbols('i k')
display_equation('I_{ij}I_{kj}', sp.Piecewise((1, sp.Eq(i, k)), (0, sp.Ne(i, k))))
display_equation('\delta_{kj}', sp.Piecewise((1, sp.Eq(i, k)), (0, sp.Ne(i, k))))
display_full_latex('\\therefore')
display_equation('I_{ij}I_{kj}', KroneckerDelta(i, k))
```


\[I_{ij}I_{kj} = \begin{cases} 1 & \text{for}\: i = k \\0 & \text{otherwise} \end{cases}\]



\[\delta_{kj} = \begin{cases} 1 & \text{for}\: i = k \\0 & \text{otherwise} \end{cases}\]



\[\therefore\]



\[I_{ij}I_{kj} = \delta_{i k}\]



```python
i, j, k = sp.symbols('i j k')
l = sp.IndexedBase('l')
sp.x_bar_j = l[i, j] * x[i]
display_equation('\\bar{x_j}', x_bar_j)
f = x_bar_j * l[j, k]
display_equation('l_{jk}\\bar{x_j}', f)
f = f.subs(l[i, j] * l[j, k], KroneckerDelta(i, k))
display_equation('l_{jk}\\bar{x_j}', f)
f = f.subs(i, k)
display_equation('l_{jk}\\bar{x_j}', f)
```

$\implies x_k = l_{jk}\bar{x_j}$

Note: $I_{ij}I_{kj}$ is actually the dot product of the ith row by the kth row, while $I_{ij}I_{jk}$ is the product of the jth row & the jth column. In other words, 
* $I_{ij}I_{kj}$ is the dot product of two $\parallel$ vectors, divided by the product of the norms of the vectors.
* $I_{ij}I_{jk}$ is the dot product of two $\bot$ vectors, divided by the product of the norms of the vectors.

$v_i = l_iv \tag{1a}$
$\bar{v_j} = v\bar{l_j} \tag{1b}$

$\bar{l_j} = l_{ij}l_i \tag{2a}$
$l_i = l_{ij}\bar{l_j} \tag{2b}$


```python
i, j, v_i, v_j, v_bar_i, v_bar_j, l_bar_i, l_bar_j = sp.symbols('i j v_i v_j \\bar{v_i} \\bar{v_j} \\bar{l_i} \\bar{l_j}')
l_i, l_j = sp.symbols('l_i l_j')
v, l = sp.symbols('v l')
l_ij = sp.symbols('l_{ij}')

f = v * l_bar_j # From eq 1b
display_equation('\\bar{v_j}', f)
print('Substitute eq 2a')
f = f.subs(l_bar_j, l_ij * l_i) # From eq 2a
display_equation('\\bar{v_j}', f)
print('Substitute eq 1a')
f = f.subs(v, v_i / l_i) # From eq 1a
display_equation('\\bar{v_j}', f)
```

### Rank of a Tensor

| Tensor | Symbol | Rank |
|---|---|---|
| Scalar | $A$ | zero |
| Contravariant Tensor | $B^i$ | 1 |
| Covariant Tensor | $C_k$ | 1 |
| Covariant Tensor | $D_{ij}$ | 2 |
| Mixed Tensor | $E^{il}_{jkl}$ | 3 |

In an n-dimensional space, the number of components of a tensor of rank $r$ is $n^r$.

**First Order**  
$a_i = l_{ij}\bar{a}_j$  
$\bar{a}_p = l_{ip}a_i$

**Second Order**  
$\bar{a}_{pq} = l_{ip}l_{jp}a_{ij}$  

**Any Order**  
$\bar{a}_{pqrs}... = l_{ip}l_{jp}l_{kr}l_{ls}...a_{ijkl}$

#### Indexed variables


```python
Delta = symbols('\\Delta')
n, m = symbols('n m', integer=True)
i = Idx('i', (1, n))
j = Idx('j', (1, n))
k = Idx('k', (0, n))
A = IndexedBase('A', shape=(m, n))
a = IndexedBase('a', shape=(m, n))
e = Eq(a[i,j] *  A[k,j], Delta * KroneckerDelta(i, k))
display(e)
print('Using k = i')
display(e.subs(k, i))
```


```python
a, b, c, d, e, f = symbols('a:f')
u1, u2, u3, v1, v2, v3 = symbols('u_1 u_2 u_3 v_1 v_2 v_3')
E = ReferenceFrame('E')
u = u1 * E.x + u2 * E.y + u3 * E.z
v = v1 * E.x + v2 * E.y + v3 * E.z
```

#### Dot Product


```python
u & v
```

#### Cross Product


```python
u ^ v
```

#### The Three Scalar Invariants of a Second Order Tensor


```python
n = 3
i = Idx('i', (1, n))
j = Idx('j', (1, n))
p = Idx('p', (1, n))
q = Idx('q', (1, n))
a = IndexedBase('a', shape=(n, n))
```

1. $a_{ii}$


```python
Sum(a[i,i], (i, 1, n)).doit()
```

2. $\frac{1}{2} (a_{ii}a_{jj} - a_{ij}a_{ji})$


```python
S('1/2') * Sum(a[i,i] * a[i,i] - a[i,j] * a[j,i] , (i, 1, n), (j, 1, n)).doit()
```

3. $|a_{ij}|$

#### Singular Tensor  
$|a_{ij}| = 0$

#### Tensor Product


```python
from sympy.abc import x,y,z,t
n = 3
a = IndexedBase('a', shape=(n))
b = IndexedBase('b', shape=(n))
```


```python
from sympy.physics.quantum import TensorProduct
n = 2
A = zeros(1, n)
B = zeros(1, n)
#A = Array([a[1], a[2]])
#B = Array([b[1], b[2]])
#display(A.rank())
#display(B.rank())
#tensorproduct(A, B).rank()

for i in range(n):
    A[i] = a[i + 1]
    B[i] = b[i + 1]

display(A.T)
display(B.T)
TensorProduct(A.T, B.T)
```


```python
P = np.array([[2, 4],[5, 6]])
B = np.array([[1, 2],[3, 5]])
C = np.array([[1, 2],[3, 5]])
```


```python
3*P - 2*P
```


```python
A = np.array([[2, 3, 1],[4, 2, 1], [3, 1, 1]])
B = np.array([2, 1, 3]).T
A@B
```


```python
E = CoordSys3D('E')
U = u_1 * E.i + u_2 * E.j + u_3 * E.k
V = v_1 * E.i + v_2 * E.j + v_3 * E.k
H = h_1 * E.i + h_2 * E.j + h_3 * E.k
```


```python
U.outer(V)
```


```python
R = dyad & H
R.simplify()
```


```python
V & H
```


```python
U * (V & H)
```


```python
u, v, w = sympy.symbols('u v w')
#E.x = Function('x')(u, v, w)
#E.y = Function('y')(u, v, w)
#E.z = Function('z')(u, v, w)
x, y, z = [E.x, E.y, E.z]
r = x * E.i + y * E.j + z * E.k
#r = E.x * E.i + E.y * E.j + E.z * E.k
r
```


```python
u, v = [x, y] # Assumed
#u, v = [x, x] # Assumed
#u, v = [y, x] # Assumed All three to prove a point from [2]
#u = Function('u')(x, y, z)
#v = Function('v')(x, y, z)
#w = Function('w')(x, y, z)
```


```python
delop = Del()
a = delop(u)
a
```


```python
diff(r, u)
```


```python
res = diff(r, u).doit() & delop(u)
res
```


```python
res.doit()
```


```python
res = diff(r, u).doit() & delop(v)
res
```


```python
res.doit()
```


```python

```


```python

```


```python
x, y, xd, xdd, yd, ydd = sympy.symbols('x y \dot{x} \ddot{x} \dot{y} \ddot{y}')
x1, x2, xp1, xp2 = sympy.symbols('x_1 x_2 x^{\prime}_1 x^{\prime}_2')
Ap1, Ap2 = sympy.symbols('A^{\prime1} A^{\prime2}')
r, theta = sympy.symbols('r theta')
```


```python
x1 = x
x2 = y
xp1 = r
xp2 = theta
x = r*cos(theta)
y = r*sin(theta)
```


```python
eq = Eq(r * sin(theta), x)
solve(eq, theta)
```


```python

```


```python

```


```python

```


```python
E = CoordSys3D('E')
Ap = sympy.symbols('A^\prime')
EA = E.orient_new_axis('E_A', theta, E.k)
EAp = EA.orient_new_axis('EA^\prime', alpha, EA.k)
A = EA.i + EA.j + EA.k
A = express(A, E)
Ap = EAp.i + EAp.j + EAp.k
Ap = express(Ap, E)
EA.rotation_matrix(EAp)

#EA.position_wrt(E)
#EA.to_matrix(E)
#EA.origin.express_coordinates(E)
#type(EA.origin)
#A.__dir__()

#EA.rotation_matrix(E)
```


```python
N = CoordSys3D('N')
A = 5 * N.i + 3 * N.j
ang = inRadians(150)
Rot = rotation(ang)
res = Rot @ A.to_matrix(N)
res = res.evalf(2)
matrix_to_vector(res, N)
```


```python
E = CoordSys3D(E)
etra1, etra2, eco1, eco2 = sympy.symbols('\vec{E}^1 \vec{E}^2 \vec{E}_1 \vec{E}_2')
A, Atra1, Atra2, Ax, Ay = sympy.symbols('\vec{A} A^1 A^2 A_x A_y')

eco1 = E.i + 3 * E.j
eco2 = 4 * E.i
A = 7 * E.i + 2 * E.j
eqns = [Atra1 * eco1.components[E.i] + Atra2 * eco2.components[E.i] - A.components[E.i], Atra1 * eco1.components[E.j] + Atra2 * 0 - A.components[E.j] ]
m = sympy.linear_eq_to_matrix(eqns, [Atra1, Atra2])
res = linsolve(eqns, Atra1, Atra2)
(Atra1, Atra2) = next(iter(res))
angle = acos((eco1 & eco2) / (eco1.magnitude() * eco2.magnitude()))
angle = inDegrees(angle)
angle = angle.evalf()
t = 90 - angle
etra1mag = 1/(eco1.magnitude() * cos(inRadians(t))).evalf()
etra2mag = 1/(eco2.magnitude() * cos(inRadians(t))).evalf()
etra1 = etra1mag * E.j
etra2 = etra2mag * cos(t) * E.i + etra2mag * sin(t) * E.j
A_angle = atan(A.components[E.j] / A.components[E.i])
eco1A_ang = acos((eco1 & A) / (eco1.magnitude() * A.magnitude()))
eco2A_ang = acos((eco2 & A) / (eco2.magnitude() * A.magnitude()))
l1 = A.magnitude() * cos(eco1A_ang)
Acon1 = l1/(cos(t) * etra1mag)
Acon1 = Acon1.evalf()
l2 = A.magnitude() * cos(eco2A_ang)
Acon2 = l2/(cos(t) * etra2mag).evalf()
```


```python
Atra1, Atra2
```


```python
m
```


```python

t
```


```python
etra1mag
```


```python
etra2mag
```


```python
#l = A.magnitude() * cos(A_angle)
A_angle.evalf()
```


```python
inDegrees(e1A_ang).evalf()
```


```python
Acon1
```


```python
Acon2
```


```python
l1.evalf(), l2
```


```python
etra2
```


```python
A & eco1, A & eco2
```


```python
A & etra1, A & etra2
```


```python

```


```python

```


```python
x, y = symbols('x y')
u = [4, 8, 16, 32]
v = [-16, -12, -8, -4, 0, 4, 8]

curves = None

for i in range(len(u) - 1):
    pu = plot(u[i]/x, show=False, xlim=(0, 4), ylim=(0, 20)) #label=str(u[i]), legend=True
    
    if (curves == None):
        curves = pu
    else:
        curves.append(pu[0])
    
for i in range(len(v) - 1):
    pv = plot(x**2 - v[i], show=False, xlim=(0, 4), ylim=(0, 20)) #label=str(v[i]), legend=True
    
    if (curves == None):
        curves = pv
    else:
        curves.append(pv[0])
    
curves.show()
```


```python

```

843fcf42b93417fa680b6218a0cc027578d3d08c
ipynb
Jupyter Notebook
Tensor Analysis.ipynb
Valentine-Efagene/Jupyter-Notebooks
91a1d98354a270d214316eba21e4a435b3e17f5d
MIT
Tensor Analysis.ipynb
Valentine-Efagene/Jupyter-Notebooks
91a1d98354a270d214316eba21e4a435b3e17f5d
MIT
Tensor Analysis.ipynb
Valentine-Efagene/Jupyter-Notebooks
91a1d98354a270d214316eba21e4a435b3e17f5d
MIT
Qwen/Qwen-72B
1. YES
2. YES