https://lcvmwww.epfl.ch/~lcvm/dna_teaching_05_06/exercises/ex5.pdf

$$E(\theta, \phi, \lambda) = \frac12 \theta^2 + \frac12(\phi - \theta)^2 + \lambda(\cos\theta + \cos\phi)$$


```python
from math import sin, cos
def energy(theta, phi, lam):
    return 0.5 * theta ** 2 + 0.5 * (phi - theta)**2 + lam * (cos(theta) + cos(phi))
```

The equilibrium condition is given by 
$$0 = \frac{\partial E}{\partial \theta} = \theta + \theta - \phi - \lambda\sin\theta = 2\theta-\phi- \lambda\sin\theta,\quad
  0 = \frac{\partial E}{\partial \phi} = \phi - \theta - \lambda\sin\phi$$



```python
def F(theta, phi, lam):
    return 2*theta - phi - lam * sin(theta), phi - theta - lam * sin(phi)
```

The Jacobian of F is given by 
$$J = \begin{pmatrix}2-\lambda\cos(\theta) & -1\\-1 & 1-\lambda\cos\phi\end{pmatrix}$$


```python
import numpy as np
def J(theta, phi, lam):
    return np.array([[2-lam*cos(theta), -1], [-1, 1-lam*cos(phi)]])
J(0,0,0)
```




    array([[ 2., -1.],
           [-1.,  1.]])



In the case of the straight rod (\theta = \phi = 0) J(0,0,\lambda) is singular when $$\lambda^2 - 3\lambda + 1 = 0$$
thus when $\lambda = \frac12(3\pm\sqrt5)$


```python
J(0,0,0.5*(3-5**.5))
```




    array([[ 1.61803399, -1.        ],
           [-1.        ,  0.61803399]])




```python
J(0,0,0.5*(3+5**.5))
```




    array([[-0.61803399, -1.        ],
           [-1.        , -1.61803399]])



The null spaces are spanned by $(1, \frac12(\sqrt5+1))$ and $(\frac12(\sqrt5+1), -1)$ respectively


```python
J(0,0,0.5*(3-5**.5)) @ np.array([1,0.5*(5**0.5+1)]), J(0,0,0.5*(3+5**.5)) @ np.array([0.5*(5**0.5+1), -1])
```




    (array([0.00000000e+00, 2.22044605e-16]),
     array([-2.22044605e-16,  0.00000000e+00]))



## Stability
The hessian is identical to J in this case. Its eigenvalues are given by
$$\mu_1\mu_2 = \det J(0,0,\lambda) = \lambda^2 - 3\lambda+1\quad\text{and}$$
$$\mu_1 + \mu_2 = \operatorname{tr}J(0,0,\lambda) = 3-2\lambda$$
Thus for $\lambda < \frac12(3-\sqrt5)$ and $\lambda > \frac12(3+\sqrt5)$ the eigenvalues have the same sign, and in between they have opposite signs. For $\lambda < \frac12(3-\sqrt5)$ their sum is positive thus the two eigenvalues are positive and for $\lambda > \frac12(3+\sqrt5)$ their sum is negative and thus they are both negative.


```python
np.linalg.eig(J(0,0,0.3))
```




    (array([2.31803399, 0.08196601]),
     array([[ 0.85065081,  0.52573111],
            [-0.52573111,  0.85065081]]))




```python
np.linalg.eig(J(0,0,1.5))
```




    (array([ 1.11803399, -1.11803399]),
     array([[ 0.85065081,  0.52573111],
            [-0.52573111,  0.85065081]]))




```python
np.linalg.eig(J(0,0,2.7))
```




    (array([-0.08196601, -2.31803399]),
     array([[ 0.85065081,  0.52573111],
            [-0.52573111,  0.85065081]]))




```python
np.linalg.eig(J(0,0,0.5*(3-5**.5)))
```




    (array([2.23606798e+00, 1.11022302e-16]),
     array([[ 0.85065081,  0.52573111],
            [-0.52573111,  0.85065081]]))



# Bifurcation Shape


```python
from sympy import *
eps, phi, the, lam, E = symbols('ε φ θ λ E')
```


```python
phi_ = symbols(['φ_%d' % i for i in range(4)])
the_ = symbols(['θ_%d' % i for i in range(4)])
lam_ = symbols(['λ_%d' % i for i in range(4)])
phi_eps = sum(eps**i * phi_[i] for i in range(4))
the_eps = sum(eps**i * the_[i] for i in range(4))
lam_eps = sum(eps**i * lam_[i] for i in range(4))
the_eps
```




$\displaystyle ε^{3} θ_{3} + ε^{2} θ_{2} + ε θ_{1} + θ_{0}$




```python
E = (the**2 + (phi-the)**2)/2 + lam*(cos(phi)+cos(the))
```


```python
F1, F2 = diff(E, the), diff(E, phi)
```


```python
F1_eps = F1.subs([(the, the_eps), (phi, phi_eps), (lam, lam_eps)])
```


```python
F1_eps.subs(eps,0)
```




$\displaystyle 2 θ_{0} - λ_{0} \sin{\left(θ_{0} \right)} - φ_{0}$




```python
diff(F1_eps, eps).subs(eps, 0)
```




$\displaystyle - θ_{1} λ_{0} \cos{\left(θ_{0} \right)} + 2 θ_{1} - λ_{1} \sin{\left(θ_{0} \right)} - φ_{1}$




```python
conditions = [
    diff(diff(E, x).subs([(the, the_eps), (phi, phi_eps), (lam, lam_eps)]), eps, i).subs(eps, 0) * factorial(i)
    for i in range(4)
    for x in (phi, the)
]; conditions
```




    [-θ_0 - λ_0*sin(φ_0) + φ_0,
     2*θ_0 - λ_0*sin(θ_0) - φ_0,
     -θ_1 - λ_0*φ_1*cos(φ_0) - λ_1*sin(φ_0) + φ_1,
     -θ_1*λ_0*cos(θ_0) + 2*θ_1 - λ_1*sin(θ_0) - φ_1,
     -4*θ_2 + 2*λ_0*φ_1**2*sin(φ_0) - 4*λ_0*φ_2*cos(φ_0) - 4*λ_1*φ_1*cos(φ_0) - 4*λ_2*sin(φ_0) + 4*φ_2,
     2*θ_1**2*λ_0*sin(θ_0) - 4*θ_1*λ_1*cos(θ_0) - 4*θ_2*λ_0*cos(θ_0) + 8*θ_2 - 4*λ_2*sin(θ_0) - 4*φ_2,
     -36*θ_3 + 6*λ_0*φ_1**3*cos(φ_0) + 36*λ_0*φ_1*φ_2*sin(φ_0) - 36*λ_0*φ_3*cos(φ_0) + 18*λ_1*φ_1**2*sin(φ_0) - 36*λ_1*φ_2*cos(φ_0) - 36*λ_2*φ_1*cos(φ_0) - 36*λ_3*sin(φ_0) + 36*φ_3,
     6*θ_1**3*λ_0*cos(θ_0) + 18*θ_1**2*λ_1*sin(θ_0) + 36*θ_1*θ_2*λ_0*sin(θ_0) - 36*θ_1*λ_2*cos(θ_0) - 36*θ_2*λ_1*cos(θ_0) - 36*θ_3*λ_0*cos(θ_0) + 72*θ_3 - 36*λ_3*sin(θ_0) - 36*φ_3]




```python
cond_1 = [c.subs([(the_[0], 0), (phi_[0],0), (lam_[0], (3-sqrt(5))/2)]) for c in conditions]; cond_1
```




    [0,
     0,
     -θ_1 - φ_1*(3/2 - sqrt(5)/2) + φ_1,
     -θ_1*(3/2 - sqrt(5)/2) + 2*θ_1 - φ_1,
     -4*θ_2 - 4*λ_1*φ_1 - 4*φ_2*(3/2 - sqrt(5)/2) + 4*φ_2,
     -4*θ_1*λ_1 - 4*θ_2*(3/2 - sqrt(5)/2) + 8*θ_2 - 4*φ_2,
     -36*θ_3 - 36*λ_1*φ_2 - 36*λ_2*φ_1 + 6*φ_1**3*(3/2 - sqrt(5)/2) - 36*φ_3*(3/2 - sqrt(5)/2) + 36*φ_3,
     6*θ_1**3*(3/2 - sqrt(5)/2) - 36*θ_1*λ_2 - 36*θ_2*λ_1 - 36*θ_3*(3/2 - sqrt(5)/2) + 72*θ_3 - 36*φ_3]




```python
from sympy.solvers.solveset import linsolve
linsolve(cond_1[:4], (the_[1], phi_[1]))

```




$\displaystyle \left\{\left( \frac{φ_{1} \left(-1 + \sqrt{5}\right)}{2}, \  φ_{1}\right)\right\}$




```python
cond_2 = [c.subs([(phi_[1],1), (the_[1], (sqrt(5)-1)/2)]).simplify() for c in cond_1]
#cond_2[4] /= 2
#cond_2[5] /= (1+sqrt(5))
#cond_2[5] = cond_2[5].simplify()
cond_2
```




    [0,
     0,
     0,
     0,
     -4*θ_2 - 4*λ_1 - 2*φ_2 + 2*sqrt(5)*φ_2,
     2*θ_2 + 2*sqrt(5)*θ_2 - 2*sqrt(5)*λ_1 + 2*λ_1 - 4*φ_2,
     -36*θ_3 - 36*λ_1*φ_2 - 36*λ_2 - 18*φ_3 + 18*sqrt(5)*φ_3 - 3*sqrt(5) + 9,
     -36*θ_2*λ_1 + 18*θ_3 + 18*sqrt(5)*θ_3 - 18*sqrt(5)*λ_2 + 18*λ_2 - 36*φ_3 - 33 + 15*sqrt(5)]




```python
linsolve(cond_2[:6], (the_[2], phi_[2], lam_[1]))
```




$\displaystyle \left\{\left( \frac{φ_{2} \left(-1 + \sqrt{5}\right)}{2}, \  φ_{2}, \  0\right)\right\}$




```python
cond_3 = [c.subs([(phi_[2],1), (the_[2], (sqrt(5)-1)/2), (lam_[1], 0)]).simplify() for c in cond_2]
cond_3
```




    [0,
     0,
     0,
     0,
     0,
     0,
     -36*θ_3 - 36*λ_2 - 18*φ_3 + 18*sqrt(5)*φ_3 - 3*sqrt(5) + 9,
     18*θ_3 + 18*sqrt(5)*θ_3 - 18*sqrt(5)*λ_2 + 18*λ_2 - 36*φ_3 - 33 + 15*sqrt(5)]




```python
linsolve(cond_3, (the_[3], phi_[3], lam_[2]))
```




$\displaystyle \left\{\left( - \frac{φ_{3} \left(1 - \sqrt{5}\right)}{2} - \frac{1}{4} + \frac{7 \sqrt{5}}{60}, \  φ_{3}, \  \frac{1}{2} - \frac{\sqrt{5}}{5}\right)\right\}$




```python

```


```python

```


```python

```


```python

```


```python

```

5235b8745a3f693d6b60420bfdfb0717752dc243
ipynb
Jupyter Notebook
TwoBars.ipynb
leovt/cosserat
7a03c2768365f39ae6a5a0592a8539a5c7c269d9
MIT
TwoBars.ipynb
leovt/cosserat
7a03c2768365f39ae6a5a0592a8539a5c7c269d9
MIT
TwoBars.ipynb
leovt/cosserat
7a03c2768365f39ae6a5a0592a8539a5c7c269d9
MIT
Qwen/Qwen-72B
1. YES
2. YES