```python
import sympy
sympy.init_printing()
import verify_3
```

# Section 1 - Linear Alebra Types

You can define both a vector and a matrix by creating an instance of Matrix from an array


```python
def demo_1a():
    
    vector = sympy.Matrix([1,2,3])
    matrix = sympy.Matrix([[1,2],[3,5]])
    display(vector)
    display(matrix)
demo_1a()
```


$\displaystyle \left[\begin{matrix}1\\2\\3\end{matrix}\right]$



$\displaystyle \left[\begin{matrix}1 & 2\\3 & 5\end{matrix}\right]$


Transpositions


```python
def demo_1b():
    
    vector = sympy.Matrix([1,2,3])
    matrix = sympy.Matrix([[1,2],[3,5]])
    display(vector)
    display(vector.transpose())
    display(matrix)
    display(matrix.transpose())
demo_1b()
```


$\displaystyle \left[\begin{matrix}1\\2\\3\end{matrix}\right]$



$\displaystyle \left[\begin{matrix}1 & 2 & 3\end{matrix}\right]$



$\displaystyle \left[\begin{matrix}1 & 2\\3 & 5\end{matrix}\right]$



$\displaystyle \left[\begin{matrix}1 & 3\\2 & 5\end{matrix}\right]$


Access individual entries


```python
def demo_1c():
    
    matrix = sympy.Matrix([[1,2],[3,5]])
    display(matrix)
    display(matrix[3]) # Read
    display(matrix[1,0])
    matrix[3] = 7 # Write
    display(matrix)
demo_1c()
```

__Exercise a__

Create the column vector $\left[x,y\right]^T$


```python
def exercise_1a():
    
    x = sympy.Symbol('x')
    y = sympy.Symbol('y')
    
    # Enter answer here
    answer = 0
    display(answer)
    print(verify_3.verify_1a(answer))
exercise_1a()
```

Slightly more interesting example, Jacobian matrix for polar coordinates


```python
def demo_1d():
    
    r = sympy.Symbol('r', positive=True)
    q = sympy.Symbol('theta', positive=True)
    x = r*sympy.cos(q)
    y = r*sympy.sin(q)
    jac = sympy.Matrix([[var1.diff(var2) for var1 in [x,y]] for var2 in [r,q]])
    display(jac)

demo_1d()
```


$\displaystyle \left[\begin{matrix}\cos{\left(\theta \right)} & \sin{\left(\theta \right)}\\- r \sin{\left(\theta \right)} & r \cos{\left(\theta \right)}\end{matrix}\right]$


# Section 2 - Multiplication

To multiply matrices, use the operator *


```python
def demo_2a():
    
    matrix_a = sympy.Matrix([[1,1],[1,1]])
    matrix_b = sympy.Matrix([[1,-1],[-1,1]])
    display(matrix_a)
    display(matrix_b)
    display(matrix_a*matrix_b)
demo_2a()
```


$\displaystyle \left[\begin{matrix}1 & 1\\1 & 1\end{matrix}\right]$



$\displaystyle \left[\begin{matrix}1 & -1\\-1 & 1\end{matrix}\right]$



$\displaystyle \left[\begin{matrix}0 & 0\\0 & 0\end{matrix}\right]$


Similarly, matrices and vectors


```python
def demo_2b():
    
    matrix = sympy.Matrix([[1,2],[3,4]])
    vector = sympy.Matrix([1,1])
    display(matrix)
    display(vector)
    display(matrix*vector)
demo_2b()
```


$\displaystyle \left[\begin{matrix}1 & 2\\3 & 4\end{matrix}\right]$



$\displaystyle \left[\begin{matrix}1\\1\end{matrix}\right]$



$\displaystyle \left[\begin{matrix}3\\7\end{matrix}\right]$


When multiplying vectors, and order and orientation is important


```python
def demo_1c():
    
    vector_a = sympy.Matrix([1,2])
    vector_b = sympy.Matrix([3,4])
    display(vector_a)
    display(vector_b)
    display(vector_a.transpose()*vector_b)
    display(vector_a*vector_b.transpose())
demo_1c()
```


$\displaystyle \left[\begin{matrix}1\\2\end{matrix}\right]$



$\displaystyle \left[\begin{matrix}3\\4\end{matrix}\right]$



$\displaystyle \left[\begin{matrix}11\end{matrix}\right]$



$\displaystyle \left[\begin{matrix}3 & 4\\6 & 8\end{matrix}\right]$


__Exercise a__

Calculate the product $R \cdot b$ where $R$ is the 2D rotation matrix and $b = [1,0]$


```python
def exercise_2a():
    
    q = sympy.Symbol('theta', positive=True)# Rotation angle
    b = sympy.Matrix([1,0])
    R = sympy.Matrix([[sympy.cos(q), -sympy.sin(q)],[sympy.sin(q), sympy.cos(q)]])
    display(R,b)
    
    # Enter answer here
    answer = 0
    display(answer)
    print(verify_3.verify_2a(answer))
exercise_2a()
```

# Section 3 - Inversion and Null Spaces

To get the matrix use .inv()


```python
def demo_3a():
    
    matrix = sympy.Matrix([[1,2],[3,4]])
    display(matrix)
    display(matrix.inv())
demo_3a()
```


$\displaystyle \left[\begin{matrix}1 & 2\\3 & 4\end{matrix}\right]$



$\displaystyle \left[\begin{matrix}-2 & 1\\\frac{3}{2} & - \frac{1}{2}\end{matrix}\right]$


__Exercise a__

Solve the following system of equations $A \cdot x = b$


```python
def exercise_3a():
    
    A = sympy.Matrix([[1,2],[3,4]])
    b = sympy.Matrix([5,6])
    
    sol = 0
    print(verify_3.verify_3a(sol))
exercise_3a()
```

    False


If the matrix is not invertible, we can find its null space


```python
def demo_3b():
    
    matrix = sympy.Matrix([[1,1],[1,1]])
    display(matrix)
    display(matrix.nullspace())
demo_3b()
```


$\displaystyle \left[\begin{matrix}1 & 1\\1 & 1\end{matrix}\right]$



$\displaystyle \left[ \left[\begin{matrix}-1\\1\end{matrix}\right]\right]$


# Section 4 - Determinant

Taking the determinant


```python
def demo_4a():
    
    matrix = sympy.Matrix([[1,2],[3,4]])
    display(matrix)
    display(matrix.det())
demo_4a()
```

__Exercise a__

Find x for which the following matrix is singular


```python
def exercise_4a():
    
    x = sympy.Symbol('x', positive=True)
    matrix = sympy.Matrix([[1,2-x],[3,4]])
    
    answer = 0
    print(verify_3.verify_4a(answer))
    
exercise_4a()
```

    False


# Section 5 - Eigenvalue system

To find the eigenvalues


```python
def demo_5a():
    
    matrix = sympy.Matrix([[1,2],[3,4]])
    display(matrix.eigenvals())

demo_5a()
```

The first number is the eigenvalue, and the second number is the multiplicity.

To get the eigenvectors


```python
def demo_5b():
    
    matrix = sympy.Matrix([[1,2],[3,4]])
    display(matrix.eigenvects())

demo_5b()
```


$\displaystyle \left[ \left( \frac{5}{2} - \frac{\sqrt{33}}{2}, \  1, \  \left[ \left[\begin{matrix}- \frac{2}{- \frac{3}{2} + \frac{\sqrt{33}}{2}}\\1\end{matrix}\right]\right]\right), \  \left( \frac{5}{2} + \frac{\sqrt{33}}{2}, \  1, \  \left[ \left[\begin{matrix}- \frac{2}{- \frac{\sqrt{33}}{2} - \frac{3}{2}}\\1\end{matrix}\right]\right]\right)\right]$


In some cases it is also useful to find the characteristic polynomial


```python
def demo_5c():
    
    x = sympy.Symbol('x')
    matrix = sympy.Matrix([[1,2],[3,4]])
    display(matrix.charpoly('x'))
demo_5c()
```

__Exercise a__

A series of aligned point masses, with springs connecting the two adjacent pairs of masses. Their equations of motion can be described by 

$\ddot{\vec{x}} = \vec{M} \cdot \vec{x}$. Find the eigensystem (eigenvectors + eigenvalues) for the matrix $M$. What does each of the modes represent?


```python
def exercise_5a():
    
    M = sympy.Matrix([[-1,1,0],
                      [1,-2,1],
                      [0,1,-1]])
    
    display(M.eigenvects())
    
    answer = 0
    print(verify_3.verify_5a(answer))
exercise_5a()
```


$\displaystyle \left[ \left( -3, \  1, \  \left[ \left[\begin{matrix}1\\-2\\1\end{matrix}\right]\right]\right), \  \left( -1, \  1, \  \left[ \left[\begin{matrix}-1\\0\\1\end{matrix}\right]\right]\right), \  \left( 0, \  1, \  \left[ \left[\begin{matrix}1\\1\\1\end{matrix}\right]\right]\right)\right]$


    False


# Section 6 - Matrix functions

sympy interprets functions on matrices as matrix functions (rather than operations on the 


```python
def demo_6a():
    
    matrix = sympy.Matrix([[0,-1],[1,0]])
    display(matrix)
    display(sympy.exp(matrix))
    
demo_6a()
```


$\displaystyle \left[\begin{matrix}0 & -1\\1 & 0\end{matrix}\right]$



$\displaystyle \left[\begin{matrix}\cos{\left(1 \right)} & - \sin{\left(1 \right)}\\\sin{\left(1 \right)} & \cos{\left(1 \right)}\end{matrix}\right]$



```python
def demo_6b():
    
    matrix = sympy.Matrix([[0,-1],[1,0]])
    display(matrix)
    display(matrix**2)
    
demo_6b()
```


$\displaystyle \left[\begin{matrix}0 & -1\\1 & 0\end{matrix}\right]$



$\displaystyle \left[\begin{matrix}-1 & 0\\0 & -1\end{matrix}\right]$


__Exercise a__

In this problem we will explore the generator to the Lorentz boost. Calculate $\exp \left(\psi Z\right)$ ,where $Z$ is the anti diagonal unit matrix and $\tanh \psi = \beta$ is the rapidity.


```python
def exercise_6a():
    
    psi = sympy.Symbol('psi', positive=True) # Rapidity
    gen = sympy.Matrix([[0,-1],[-1,0]])
    
    # Enter anser here
    temp = sympy.exp(psi*gen)
    temp.simplify()
    display(temp)
    answer = sympy.Rational(0,1)
    print(verify_3.verify_6a(answer))
    
exercise_6a()
```


$\displaystyle \left[\begin{matrix}\frac{e^{\psi}}{2} + \frac{e^{- \psi}}{2} & - \frac{e^{\psi}}{2} + \frac{e^{- \psi}}{2}\\- \frac{e^{\psi}}{2} + \frac{e^{- \psi}}{2} & \frac{e^{\psi}}{2} + \frac{e^{- \psi}}{2}\end{matrix}\right]$


    False


# Section 7 - Worked example

## Tidal Field

In this section we derive the tidal field of a point mass. Close to some distant point the potential behaves roughly as $\varphi \approx \frac{1}{2} \vec{x} T \vec{x}$ where $T$ is the tidal tensor.


```python
def demo_7a():
    
    x = sympy.Symbol('x', real=True)
    y = sympy.Symbol('y', real=True)
    z = sympy.Symbol('z', real=True)
    G = sympy.Symbol('G', positive=True) # Gravitation constant
    M = sympy.Symbol('M', positive=True) # Mass
    a = sympy.Symbol('a', positive=True) # Distance
    
    potential = G*M/sympy.sqrt(x**2+y**2+z**2)
    
    display(potential)
    
    hessian = sympy.Matrix([[potential.diff(var1).diff(var2).simplify()
                             for var1 in [x,y,z]] for var2 in [x,y,z]])
    
    display(hessian)
    
demo_7a()
```

## Two Body Problem


```python
def demo_7b():
    
    r = sympy.Function('r', positive=True)
    s = sympy.Function('s', positive=True)
    qf = sympy.Function('theta', real=True)
    q = sympy.Symbol('theta', real=True)
    t = sympy.Symbol('t', real=True)
    pos = sympy.Matrix([r(t)*sympy.cos(qf(t)),r(t)*sympy.sin(qf(t))])
    G = sympy.Symbol('G', positive=True) # Gravitation constant
    M = sympy.Symbol('M', positive=True) # Mass
    l = sympy.Symbol('l', positive=True) # Angular momentum
    u = sympy.Symbol('u', real=True) # Energy
    xi = sympy.Symbol('xi', positive=True) # Auxiliary variable
    energy = -G*M/r(t)+(r(t).diff(t)**2+qf(t).diff(t)**2*r(t)**2)/2-u
    eom = pos.diff(t,2)+G*M*pos/(pos.transpose()*pos)**sympy.Rational(3,2)
    
    v = pos.diff(t)
    LRL = l*sympy.Matrix([v[1],-v[0]])-G*M*pos/sympy.sqrt((pos.transpose()*pos)[0])
    LRL = sympy.Matrix([LRL[0].simplify(), LRL[1].simplify()])
    LRL = LRL.subs(r(t),xi).subs(xi,r(t))
    
    print('conservation of angular momentum')
    temp = eom.transpose()*sympy.Matrix([sympy.sin(qf(t)), -sympy.cos(qf(t))])
    temp = temp[0].simplify()
    display(temp)
    temp = (r(t)**2*qf(t).diff(t)).diff(t).subs(sympy.solve(temp,qf(t).diff(t,2),dict=True)[0])
    display(temp)
    
    print('conservation of the LRL vector')
    temp = LRL
    display(LRL)
    temp = temp.diff(t)    
    temp = temp.subs(sympy.solve(eom,[r(t).diff(t,2), qf(t).diff(t,2)],dict=True)[0])
    temp = temp.subs(qf(t).diff(t),l/r(t)**2)
    temp = sympy.Matrix([temp[0].simplify(), temp[1].simplify()])
    display(temp)
    
    print('trajectory')
    temp = energy
    temp = temp.subs(r(t).diff(t), r(q).diff(q)*qf(t).diff(t))
    temp = temp.subs(qf(t).diff(t), l/r(q)**2)
    temp = temp.subs(r(t), r(q))
    energy_r_q = temp
    temp = temp.subs(r(q),1/s(q))
    temp = temp.doit()
    temp = temp.diff(q)
    temp = temp/s(q).diff(q)
    temp = temp.simplify()
    temp = sympy.dsolve(temp,s(q))
    temp = 1/sympy.solve(temp,s(q))[0]
    temp = temp.subs(sympy.solve(temp.diff(q).subs(q,0),sympy.Symbol('C1'),dict=True)[0])
    display(temp)
    temp = temp.subs(sympy.solve(energy_r_q.subs(r(q),temp).doit(),sympy.Symbol('C2'),dict=True)[1])
    display(temp)
    
demo_7b()
```

## Sound Waves


```python
def demo_7c():
    
    x = sympy.Symbol('x', real=True) # Position
    t = sympy.Symbol('t', real=True) # Time
    rho = sympy.Symbol('rho') # Density
    v = sympy.Symbol('v') # Velocity
    p = sympy.Symbol('p') # Pressure
    s = sympy.Symbol('s') # Entropy
    J = sympy.Symbol('J') # Mass flux
    xi = sympy.Symbol('xi') # Auxiliary variable
    gamma = sympy.Symbol('gamma', positive=True) # Adiabatic index
    
    # Hydrodynamic equations
    continuity_equation = sympy.Derivative(rho,t)+sympy.Derivative(J,x)
    momentum_conservatoin = sympy.Derivative(v,t)+v*sympy.Derivative(v,x)+sympy.Derivative(p,x)/rho
    entropy_conservation = sympy.Derivative(s,t)+v*sympy.Derivative(s,x)
    hydro_eqns = [continuity_equation,
                  momentum_conservatoin,
                  entropy_conservation]
    print('equations of motion')
    display(hydro_eqns)
    
    # Perturbation
    rho_0 = sympy.Symbol('rho_0', positive=True) # Ambient density
    p_0 = sympy.Symbol('p_0', positive=True) # Ambient pressure
    delta_rho = sympy.Symbol(r'\delta \rho', positive=True) # Density fluctuation amplitude
    delta_p = sympy.Symbol(r'\delta p', positive=True) # Pressure fluctuation amplitude
    delta_v = sympy.Symbol(r'\delta v', positive=True) # Velocity fluctuation amplitude
    omega = sympy.Symbol('omega', real=True) # Frequency
    k = sympy.Symbol('k', positive=True) # Wavenumber
    epsilon = sympy.Symbol('epsilon') # Perturbative parameter
    temp = hydro_eqns
    mode = sympy.exp(sympy.I*(omega*t-k*x))
    temp = [itm.subs({rho:rho_0 + epsilon*delta_rho*mode,
                      p:p_0 + epsilon*delta_p*mode,
                      v:epsilon*delta_v*mode,
                      J:(rho_0 + epsilon*delta_rho*mode)*(epsilon*delta_v*mode), # sympy spoonfeed
                      s:sympy.log(p_0 + epsilon*delta_p*mode) - gamma*sympy.log(rho_0 + epsilon*delta_rho*mode)}) 
            for itm in temp]
    temp = [itm.doit() for itm in temp]
    temp = [itm.diff(epsilon).subs(epsilon,0) for itm in temp]
    temp = [itm/mode for itm in temp]
    temp = [sympy.expand(itm) for itm in temp]
    temp = sympy.Matrix([[itm.diff(var) for itm in temp] for var in [delta_rho, delta_p, delta_v]])
    print('Solutions are only possible if this matrix is singular')
    display(temp)
    print('The determinant is')
    temp = temp.det().simplify()
    display(temp)
    print('And the frequencies are')
    display(sympy.solve(temp,omega))
    print('These solutions represent a left going wave, a right going wave and a density perturbation that does not move')
    
demo_7c()
```


```python

```

28c23d1036fa8cf975b61f84f0ac7ca2c9d9cd9b
ipynb
Jupyter Notebook
lesson_3.ipynb
bolverk/cyborg_math
298224dadd4218ebcc266b0a135e15595a359343
MIT
2020-02-25T22:29:21.000Z
2022-03-01T21:49:08.000Z
lesson_3.ipynb
bolverk/cyborg_math
298224dadd4218ebcc266b0a135e15595a359343
MIT
lesson_3.ipynb
bolverk/cyborg_math
298224dadd4218ebcc266b0a135e15595a359343
MIT
2020-05-14T21:02:56.000Z
2021-10-05T19:12:11.000Z
Qwen/Qwen-72B
1. YES
2. YES