```python
import numpy as np 
import sympy as sp
import scipy as sc
from scipy.optimize import fsolve
from numpy import linalg as la
```

$PROBLEM$ $1$

$Part$ $A$


```python
A = np.array([[2,5,-3],[5,7,1], [-3,1,4]])
print(A)
print("\nSince it's 3x3 Matrix we have 3 degree characterstic polynomial")
λ = sp.symbols('λ')
display(Eq(-λ**3+13*λ**2-15*λ-139,0))
```

$Part$ $B$


```python
#Let's define the obove characterstic polynomial and compute it
def func(λ):
    return -λ**3+13*λ**2-15*λ-139

λ1=fsolve(func,1)
print("λ1 =",λ1)
λ2=fsolve(func,10)
print("λ2 =",λ2)
λ3=fsolve(func,5)
print("λ3 =",λ3)

```

$Part$ $C$


```python
eigVal,eigVect=la.eig(A)
print("Eigen Values\n",eigVal)
print("\nEigenvectors\n",eigVect)
```

$Part$ $D$


```python
#for a Matrix to be orthogonal Q^T=Q^-1
Q=eigVect
print(Q.T,"\n")
print(np.linalg.inv(Q))
print("\n Thus the EigenVectors are Orthogonal")
```

$PROBLEM$ $2$

$Part$ $A$


```python
M=np.array([[-1,8,2,10,45],[8,6,-11,-23,85],[2,-11,4,-3,-7],[10,-23,-3,2,2],[45,85,-7,2,2]])
for i in range(25):
 q,r=la.qr(M)
 M=r@q
 
print("\nEigenvalues by QR factorization\n")
print(np.rint(M))
```

    
    Eigenvalues by QR factorization
    
    [[105.  -6.   0.  -0.  -0.]
     [ -6. -91.   0.   0.   0.]
     [  0.   0. -24.   0.  -0.]
     [ -0.   0.   0.  18.  -0.]
     [  0.  -0.  -0.  -0.   5.]]


$Part$ $B$


```python
eigVal,eigVect=la.eig(M)
print("Eigen Values\n",np.around(eigVal,3))
print("\nEigenvectors\n",np.around(eigVect,2))

```

    Eigen Values
     [104.99  -91.557 -23.637  18.211   4.993]
    
    Eigenvectors
     [[-1.   -0.03 -0.    0.    0.  ]
     [ 0.03 -1.   -0.    0.    0.  ]
     [-0.    0.   -1.    0.   -0.  ]
     [ 0.   -0.    0.    1.    0.  ]
     [-0.    0.   -0.   -0.    1.  ]]


$Part$ $C$


```python
print(np.trace(M))
print(sum(eigVal))
```

    12.999999999999865
    12.999999999999819


$Part$ $D$


```python
print(la.det(M))
print(np.prod(eigVal))
```

    20660378.000000045
    20660378.00000005


$PROBLEM$ $3$

$Part$ $A$

Rewriting the Matrix as Eigen Problem

\begin{pmatrix}
\frac{k}{m_{1}}-w^{2} & -\frac{k}{m_{1}} & 0\\ 
 -\frac{k}{m_{2}}& \frac{2k}{m_{2}}-w^{2} & -\frac{k}{m_{2}} \\ 
 0 & -\frac{k}{m_{1}} & \frac{k}{m_{1}}-w^{2}
\end{pmatrix}


$Part$ $B$

$Substituting $ the $ values $ of $ k, $ m1, $ m2 $ in $ the $ equation

\begin{pmatrix}
3.06*10^{28}-\lambda & -3.06*10^{28} & 0\\ 
 -1.20*10^{29} & 2.41*10^{29}-\lambda & -1.20*10^{29} \\ 
 0 & -3.06*10^{28} & 3.06*10^{28}-\lambda
\end{pmatrix}



```python
V=np.array([[3.06*10**28,-3.06*10**28,0],[-1.20*10**29,2.41*10**29,-1.20*10**29],[0,-3.06*10**28,3.06*10**28]])
print(V)
eigValues,eigVectors=la.eig(V)
print("\nThe vibration Frequenices (w^2)1, (w^2)2, (w^2)3 are \n\n",eigValues)
```

    [[ 3.06e+28 -3.06e+28  0.00e+00]
     [-1.20e+29  2.41e+29 -1.20e+29]
     [ 0.00e+00 -3.06e+28  3.06e+28]]
    
    The vibration Frequenices (w^2)1, (w^2)2, (w^2)3 are 
    
     [2.71487288e+29 3.06000000e+28 1.12712460e+26]


$Part$ $C$


```python
print("\nEigen Vectors\n\n",eigVectors)
```

    
    Eigen Vectors
    
     [[-1.25028831e-01 -7.07106781e-01  5.78059140e-01]
     [ 9.84243660e-01 -1.84988936e-17  5.75929909e-01]
     [-1.25028831e-01  7.07106781e-01  5.78059140e-01]]


$PROBLEM$ $4$

$Part$ $A$


```python
S=np.array([[20,8,-18],[8,28,12],[-18,12,4]])
# print(S)
eigVal,eigVect=la.eig(S)
print("\nPrincipal Stresses in Mega Pascals are\n\n",eigVal)
```

    
    Principal Stresses in Mega Pascals are
    
     [-12.79567553  31.67988775  33.11578778]


$Part$ $B$


```python
#To find the angles we need the normal(eigenvectors) vectors frist
#Also python normalizes the normal vectors automatically
#angles are measured from the +ve axis CCW direction.
N=eigVect
print("Normal Vectors to the Prinipal Plane\n",N)

print("\nThe Maximum principal stress is 33.11 MPa")
print("It's Unit Normal to plane is ", N[:,2])

print("\nX-axis in degrees ",np.degrees(np.arccos(N[0,2]))," or radians",np.arccos(N[0,2]))
print("Y-axis in degrees ",np.degrees(np.arccos(N[1,2]))," or radians",np.arccos(N[1,2]))
print("Z-axis in degrees ",np.degrees(np.arccos(N[2,2]))," or radians",np.arccos(N[2,2]))
```

    Normal Vectors to the Prinipal Plane
     [[-0.51481828  0.81232573  0.27402382]
     [ 0.33329183 -0.10484587  0.93697593]
     [-0.78985992 -0.57370224  0.21676496]]
    
    The Maximum principal stress is 33.11 MPa
    It's Unit Normal to plane is  [0.27402382 0.93697593 0.21676496]
    
    X-axis in degrees  74.09615133093509  or radians 1.2932218037807959
    Y-axis in degrees  20.450248042933367  or radians 0.35692416119871395
    Z-axis in degrees  77.48090601002706  or radians 1.3522969173032346


$Part$ $C$


```python
#similary the Minimum Princial stress is -12.79Mpa
#Also python normalizes the normal vectors automatically
#angles are measured from the +ve axis CCW direction.
N=eigVect
print("Normal Vectors to the Prinipal Plane\n",N)

print("\nThe Minimum principal stress is -12.79 MPa")
print("It's Unit Normal to plane is ", N[:,0])

print("\nX-axis in degrees ",np.degrees(np.arccos(N[0,0]))," or radians",np.arccos(N[0,0]))
print("Y-axis in degrees ",np.degrees(np.arccos(N[1,0]))," or radians",np.arccos(N[1,0]))
print("Z-axis in degrees ",np.degrees(np.arccos(N[2,0]))," or radians",np.arccos(N[2,0]))
```

    Normal Vectors to the Prinipal Plane
     [[-0.51481828  0.81232573  0.27402382]
     [ 0.33329183 -0.10484587  0.93697593]
     [-0.78985992 -0.57370224  0.21676496]]
    
    The Minimum principal stress is -12.79 MPa
    It's Unit Normal to plane is  [-0.51481828  0.33329183 -0.78985992]
    
    X-axis in degrees  120.9853092679676  or radians 2.111591993269645
    Y-axis in degrees  70.53130164099609  or radians 1.2310034393526612
    Z-axis in degrees  142.17242293287993  or radians 2.4813768857166476


f3661279a1b40fbc89d96c635c80195dd8e4a1d7
ipynb
Jupyter Notebook
assignments/assignment03/assignment03_solution.ipynb
eyobghiday/computational-mechanics
c64aebca5d26fdf5fc93977a5ec3dbaa4e715a48
MIT
2021-06-27T06:21:23.000Z
2022-03-17T09:48:45.000Z
assignments/assignment03/assignment03_solution.ipynb
eyobghiday/computational-mechanics
c64aebca5d26fdf5fc93977a5ec3dbaa4e715a48
MIT
assignments/assignment03/assignment03_solution.ipynb
eyobghiday/computational-mechanics
c64aebca5d26fdf5fc93977a5ec3dbaa4e715a48
MIT
Qwen/Qwen-72B
1. YES
2. YES