Solve
\begin{cases} \dot{x} &= \sigma \cdot (y - x) \\ \dot{y} &= x \cdot (r - z) - y \\ \dot{z} &= x y - b z \\ \end{cases}
for 0 < t < 40 with initial conditions
\begin{cases} x(0) = -11 \\ y(0) = -16 \\ z(0) = 22.5 \\ \end{cases}
and \sigma=10, r=28 and b=8/3, which are the classical parameters that generate the butterfly as presented by Edward Lorenz back in his seminal 1963 paper Deterministic non-periodic flow. This example’s input file ressembles the parameters, inital conditions and differential equations of the problem as naturally as possible with an ASCII file.
PHASE_SPACE x y z # Lorenz attractor’s phase space is x-y-z
end_time = 40 # we go from t=0 to 40 non-dimensional units
= 10 # the original parameters from the 1963 paper
sigma = 28
r = 8/3
b
= -11 # initial conditions
x_0 = -16
y_0 = 22.5
z_0
# the dynamical system's equations written as naturally as possible
= sigma*(y - x)
x_dot = x*(r - z) - y
y_dot = x*y - b*z
z_dot
PRINT t x y z # four-column plain-ASCII output
$ feenox lorenz.fee > lorenz.dat
$ gnuplot lorenz.gp
$ python3 lorenz.py
$ sh lorenz2x3d.sh < lorenz.dat > lorenz.html
Figure 1: The Lorenz attractor computed with FeenoX plotted with two different tools. a — Gnuplot, b — Matplotlib
Find the location of the two bobs vs time in the double pendulum in fig. 2.
Use these four different approaches:
Hamiltonian formulation with numerical derivatives
\begin{aligned} \mathcal{H}(\theta_1, \theta_2, p_1, p_2) =& - \frac{\ell_2^2 m_2 p_1^2 - 2 \ell_1 \ell_2 m_2 \cos(\theta_1-\theta_2) p_1 p_2 + \ell_1^2 (m_1+m_2) p_2^2} {\ell_1^2 \ell_2^2 m_2 \left[-2m_1-m_2+m_2\cos\Big(2(\theta_1-\theta_2)\Big)\right]} \\ & - \Big[ m_1 g \ell_1 \cos \theta_1 + m_2 g (\ell_1 \cos \theta_1 + \ell_2 \cos \theta_2) \Big] \end{aligned}
\begin{cases} \displaystyle \dot{\theta}_1 &= \displaystyle +\frac{\partial \mathcal{H}}{\partial p_1} \\ \displaystyle \dot{\theta}_2 &= \displaystyle +\frac{\partial \mathcal{H}}{\partial p_2} \\ \displaystyle \dot{p}_1 &= \displaystyle -\frac{\partial \mathcal{H}}{\partial \theta_1} \\ \displaystyle \dot{p}_2 &= \displaystyle -\frac{\partial \mathcal{H}}{\partial \theta_2} \\ \end{cases}
# the double pendulum solved by the hamiltonian formulation
# and numerically computing its derivatives
PHASE_SPACE theta1 theta2 p1 p2
VAR theta1' theta2' p1' p2'
(theta1,theta2,p1,p2) = \
H- (m1*g*l1*cos(theta1) + m2*g*(l1*cos(theta1) \
+ l2*cos(theta2))) \
- (l2^2*m2*p1^2 - 2*l1*l2*m2*cos(theta1-theta2)*p1*p2 + \
^2*(m1+m2)*p2^2)/(l1^2*l2^2*m2 \
l1* (-2*m1-m2+m2*cos(2*(theta1-theta2))))
= +derivative(H(theta1,theta2,p1',p2), p1', p1)
theta1_dot .= +derivative(H(theta1,theta2,p1,p2'), p2', p2)
theta2_dot .
= -derivative(H(theta1',theta2,p1,p2), theta1', theta1)
p1_dot .= -derivative(H(theta1,theta2',p1,p2), theta2', theta2) p2_dot .
Hamiltonian formulation with analytical derivatives
\begin{cases} \dot{\theta}_1 &= \displaystyle \frac{p_1 \ell_2 - p_2 \ell_1 \cos(\theta_1-\theta_2)}{\ell_1^2 \ell_2 [m_1 + m_2 \sin^2(\theta_1-\theta_2)]} \\ \dot{\theta}_2 &= \displaystyle \frac{p_2 (m_1+m_2)/m_2 \ell_1 - p_1 \ell_2 \cos(\theta_1-\theta_2)}{\ell_1 \ell_2^2 [m_1 + m_2 \sin^2(\theta_1-\theta_2)]} \\ \dot{p_1} &= \displaystyle -(m_1+m_2) g \ell_1 \sin(\theta_1) - c_1 + c_2 \\ \dot{p_2} &= \displaystyle -m_2 g \ell_2 \sin(\theta_2) + c_1 - c_2 \end{cases} where the expressions c_1 and c_2 are
\begin{aligned} c1 &= \frac{p_1 p_2 \sin(\theta_1-\theta_2)}{\ell_1 \ell_2 \Big[m_1+m_2 \sin(\theta_1-\theta_2)^2\Big]} \\ c2 &= \frac{\Big[ p_1^2 m_2 \ell_2^2 - 2 p_1 p_2 m_2 \ell_1 \ell_2 \cos(\theta_1-\theta_2) + p_2^2 (m_1+m_2) \ell_1^2)\Big] \sin(2 (\theta_1-\theta_2)}{ 2 \ell_1^2 \ell_2^2 \left[m_1+m_2 \sin^2(\theta_1-\theta_2)\right]^2} \end{aligned}
# the double pendulum solved by the hamiltonian formulation
# and analytically computing its derivatives
PHASE_SPACE theta1 theta2 p1 p2 c1 c2
= (p1*l2 - p2*l1*cos(theta1-theta2))/(l1^2*l2*(m1 + m2*sin(theta1-theta2)^2))
theta1_dot .= (p2*(m1+m2)/m2*l1 - p1*l2*cos(theta1-theta2))/(l1*l2^2*(m1 + m2*sin(theta1-theta2)^2))
theta2_dot .
= -(m1+m2)*g*l1*sin(theta1) - c1 + c2
p1_dot .= -m2*g*l2*sin(theta2) + c1 - c2
p2_dot .
= p1*p2*sin(theta1-theta2)/(l1*l2*(m1+m2*sin(theta1-theta2)^2))
c1 .= { (p1^2*m2*l2^2 - 2*p1*p2*m2*l1*l2*cos(theta1-theta2)
c2 .+ p2^2*(m1+m2)*l1^2)*sin(2*(theta1-theta2))/
(2*l1^2*l2^2*(m1+m2*sin(theta1-theta2)^2)^2) }
Lagrangian formulation with numerical derivatives
\begin{aligned} \mathcal{L}(\theta_1, \theta_2, \dot{\theta}_1, \dot{\theta}_2) =& \frac{1}{2} m_1 \ell_1^2 \dot{\theta}_1^2 + \frac{1}{2} m_2 \left[\ell_1^2 \dot{\theta}_1^2 + \ell_2^2 \dot{\theta}_2^2 + 2 \ell_1 \ell_2 \dot{\theta}_1 \dot{\theta}_2 \cos(\theta_1-\theta_2)\right] + \\ & m_1 g \ell_1\cos \theta_1 + m_2 g \left(\ell_1\cos \theta_1 + \ell_2 \cos \theta_2 \right) \end{aligned}
\begin{cases} \displaystyle \frac{\partial}{\partial t}\left(\frac{\partial \mathcal{L}}{\partial \dot{\theta}_1}\right) &= \displaystyle \frac{\partial \mathcal{L}}{\partial \theta_1} \\ \displaystyle \frac{\partial}{\partial t}\left(\frac{\partial \mathcal{L}}{\partial \dot{\theta}_2}\right) &= \displaystyle \frac{\partial \mathcal{L}}{\partial \theta_2} \end{cases}
# the double pendulum solved by the lagrangian formulation
# and numerically computing its derivatives
PHASE_SPACE theta1 theta2 dL_dthetadot1 dL_dthetadot2
VAR theta1' theta2' theta1_dot' theta2_dot'
(theta1,theta2,theta1_dot,theta2_dot) = {
L# kinetic energy of m1
1/2*m1*l1^2*theta1_dot^2 +
# kinetic energy of m2
1/2*m2*(l1^2*theta1_dot^2 + l2^2*theta2_dot^2 + 2*l1*l2*theta1_dot*theta2_dot*cos(theta1-theta2))
+ (
# potential energy of m1
*g * l1*cos(theta1) +
m1# potential energy of m2
*g * (l1*cos(theta1) + l2*cos(theta2))
m2) }
# there is nothing wrong with numerical derivatives, is there?
= derivative(L(theta1, theta2, theta1_dot', theta2_dot), theta1_dot', theta1_dot)
dL_dthetadot1 .= derivative(L(theta1, theta2, theta1_dot, theta2_dot'), theta2_dot', theta2_dot)
dL_dthetadot2 .
= derivative(L(theta1', theta2,theta1_dot, theta2_dot), theta1', theta1)
dL_dthetadot1_dot .= derivative(L(theta1, theta2',theta1_dot, theta2_dot), theta2', theta2) dL_dthetadot2_dot .
Lagrangian formulation with analytical derivatives
\begin{cases} 0 &= (m_1+m_2) \ell_1^2 \ddot{\theta}_1 + m_2 \ell_1 \ell_2 \ddot{\theta}_2 \cos(\theta_1-\theta_2) + m_2 \ell_1 \ell_2 \dot{\theta}_2^2 \sin(\theta_1-\theta_2) + \ell_1 (m_1+m_2) g \sin(\theta_1) \\ 0 &= m_2 \ell_2^2 \ddot{\theta}_2 + m_2 \ell_1 \ell_2 \ddot{\theta}_1 \cos(\theta_1-\theta_2) - m_2 \ell_1 \ell_2 \dot{\theta}_1^2 \sin(\theta_1-\theta_2) + \ell_2 m_2 g \sin(\theta_2) \end{cases}
# the double pendulum solved by the lagrangian formulation
# and analytically computing its derivatives
PHASE_SPACE theta1 theta2 omega1 omega2
# reduction to a first-order system
= theta1_dot
omega1 .= theta2_dot
omega2 .
# lagrange equations
0 .= (m1+m2)*l1^2*omega1_dot + m2*l1*l2*omega2_dot*cos(theta1-theta2) + \
*l1*l2*omega2^2*sin(theta1-theta2) + l1*(m1+m2)*g*sin(theta1)
m2
0 .= m2*l2^2*omega2_dot + m2*l1*l2*omega1_dot*cos(theta1-theta2) - \
*l1*l2*omega1^2*sin(theta1-theta2) + l2*m2*g*sin(theta2) m2
The combination Hamilton/Lagrange and numerical/analytical is given
in the command line as arguments $1
and $2
respectively.
# the double pendulum solved using different formulations
# parameters
end_time = 10
= 0.3
m1 = 0.2
m2 = 0.3
l1 = 0.25
l2 = 9.8
g
# inital conditions
= pi/2
theta1_0 = pi
theta2_0
# include the selected formulation
DEFAULT_ARGUMENT_VALUE 1 hamilton
DEFAULT_ARGUMENT_VALUE 2 numerical
INCLUDE double-$1-$2.fee
# output the results vs. time
PRINT t theta1 theta2 theta1_dot theta2_dot \
*sin(theta1) -l1*cos(theta1) \
l1*sin(theta1)+l2*sin(theta2) -l1*cos(theta1)-l2*cos(theta2) l1
$ for guy in hamilton lagrange; do
for form in numerical analytical; do
feenox double.fee $guy $form > double-${guy}-${form}.tsv;
m4 -Dguy=$guy -Dform=$form -Dtype=$RANDOM double.gp.m4 | gnuplot -;
done;
done
$
Figure 3: Position of the double pendulum’s m_2 solved with four (slightly) different formulations. a — Hamilton numerical, b — Hamilton analytical, c — Lagrange numerical, d — Lagrange analytical