Assume that $$ \begin{cases} \sigma > 0 \\ b > 0 \\ r > 0 \end{cases} $$
Equation to solve is
$$ \frac{ d \overline{R} }{ dt } = \begin{bmatrix} 0 \\ 0 \\ 0 \end{bmatrix} $$which means
$$ \begin{bmatrix} 0 \\ 0 \\ 0 \end{bmatrix} = \begin{bmatrix} \sigma \cdot \left( y - x \right) \\ -y + x \cdot \left( r - z \right) \\ -b \cdot z + x \cdot y \\ \end{bmatrix} \Longrightarrow \begin{cases} \sigma \cdot \left( y - x \right) &= 0 \\ -y + x \cdot \left( r - z \right) &= 0 \\ -b \cdot z + x \cdot y &= 0 \\ \end{cases} $$It's obvious that $x = y = p$, where $p$ is some new variable $$ \begin{cases} x = y &= p \\ p \cdot \left( r - z - 1 \right) &= 0 \\ -b \cdot z + p^2 &= 0 \\ \end{cases} $$
When $p = 0$ $$ \begin{cases} x &= 0 \\ y &= 0 \\ z &= 0 \end{cases} $$
Otherwise $r - z - 1 = 0$ $$ \begin{cases} z &= r - 1 \\ p^2 &= b \cdot z \end{cases} \Longrightarrow \begin{cases} z &= r - 1 \\ x &= \pm \sqrt{b \cdot z} \\ y &= \pm \sqrt{b \cdot z} \\ \end{cases} $$
Finally $$ \begin{cases} z &= r - 1 \\ x &= \pm \sqrt{b \cdot \left( r - 1 \right)} \\ y &= \pm \sqrt{b \cdot \left( r - 1 \right)} \end{cases} $$
We cannot solve the equation because of $x \cdot y$ and $x \cdot z$ so it needs to be linearized
Just kick those products $$ \frac{d\overline{R}}{dt} = \begin{bmatrix} \sigma \cdot \left( y - x \right) \\ -y + x \cdot \left( r - z \right) \\ -b \cdot z + x \cdot y \\ \end{bmatrix} $$
In matrix form $$ \frac{d\overline{R}}{dt} = \begin{bmatrix}
- \sigma, &\sigma, &0 \\
r, &-1, &0 \\
0, &0, &-b
\end{bmatrix}
\cdot
\overline{R}
\frac{d\overline{R}}{dt} = \begin{bmatrix}
- \sigma, &\sigma, &0 \\
r, &-1, &0 \\
0, &0, &-b
\end{bmatrix}
\cdot
\overline{R}
\overline{R} = \overline{R}_0 \cdot \exp{\left{ t \cdot \begin{bmatrix}
- \sigma, &\sigma, &0 \\
r, &-1, &0 \\
0, &0, &-b
\end{bmatrix}\right\}}
$$
Obvious that one sulution is $$ \lambda_1 = -b $$
It's easy to solve $$ \lambda^2 + \lambda \cdot \left( \sigma + 1 \right)+ \sigma \cdot \left( 1 - r \right) = 0 $$
$$ D = \sigma^2 + 2 \cdot \sigma + 1 - 4 \cdot \sigma + 4 \cdot \sigma \cdot r = \left( \sigma - 1 \right)^2 + 4 \cdot \sigma \cdot r $$$$ \lambda_{2, 3} = \frac{-\sigma - 1 \pm \sqrt{\left( \sigma - 1 \right)^2 + 4 \cdot \sigma \cdot r}}{2} $$Conition of a fixpoint $$ \pm \sqrt{\left( \sigma - 1 \right)^2 + 4 \cdot \sigma \cdot r} < \sigma + 1 $$
$$ r < \frac{\left( \sigma - 1 \right)^2}{4 \cdot \sigma} $$Square both parts
$$ \left( \sigma - 1 \right)^2 + 4 \cdot \sigma \cdot r < \left( \sigma + 1 \right)^2 $$$$ 4 \cdot \sigma \cdot r < 4 \cdot \sigma $$$$ r < 1 $$Fixpoints are $$ \overline{R}_0 = \begin{bmatrix} \pm \sqrt{b \cdot \left( r - 1 \right)} \\ \pm \sqrt{b \cdot \left( r - 1 \right)} \\ r + 1 \end{bmatrix} $$
Say that $\sqrt{b \cdot \left( r - 1 \right)}$ are multiplied by some constants $c \in \left\{ -1, 1 \right\}$ to ease further calculations
Also define $k = \sqrt{b \cdot \left( r + 1 \right)}$
$$ \overline{R}_0 = \begin{bmatrix} - c \cdot k \\ - c \cdot k \\ r - 1 \end{bmatrix} $$Move the center of our coordinates system to these points $$ \begin{bmatrix} x_c \\ y_c \\ z_c \end{bmatrix} = \overline{C} = \overline{R} - \overline{R}_0 = \begin{bmatrix} x + c \cdot k \\ y + c \cdot k \\ z - r + 1 \end{bmatrix} $$
We have a polynomial $$ P\left( \lambda \right) = \lambda^3
+ \lambda^2 \cdot \left( 1 + b + \sigma \right)
+ \lambda \cdot b \cdot \left( \sigma + r \right)
+ 2 \cdot b \cdot \left( r -1 \right)
= 0 $$
Since $b$ and $\sigma$ are positive and $r > 1$ out of critical value for $\left[0, 0, 0\right]$ $$ \begin{cases} \frac{P\left( \lambda \right)}{d\lambda} > 0,\qquad \lambda \ge 0 \\ P\left( \lambda \right) = 2 \cdot b \cdot \left( r -1 \right) > 0 \end{cases} \Longrightarrow P\left( \lambda \right) > 0,\qquad \lambda \in \mathbb{R}^+ $$
This means that non-negative real $\lambda$ is not a solution, thus $\forall r > 1: \lambda_1 < 0$
Other two eigenvalues should occur as a complex conjugate pair
We need to explore the critical value where $Re\left( \lambda \right) = 0$ $$ P\left( \pm i \cdot \Lambda \right) = \pm i \cdot \Lambda^3
- \Lambda^2 \cdot \left( 1 + b + \sigma \right)
\pm i \cdot \Lambda \cdot b \cdot \left( \sigma + r \right)
+ 2 \cdot b \cdot \sigma \cdot \left( r -1 \right)
= 0 $$
It doesn't matter whether $\Lambda$ is positive or negative $$ P\left( i \cdot \Lambda \right) = \left[ 2 \cdot b \cdot \sigma \cdot \left( r -1 \right) - \Lambda^2 \cdot \left( 1 + b + \sigma \right) \right]
+ i \cdot \Lambda \cdot \left[ b \cdot \left( \sigma + r \right) - \Lambda^2 \right]
\frac{2 \cdot b \cdot \sigma \cdot \left( r -1 \right)}{1 + b + \sigma} = b \cdot \left( \sigma + r \right) $$
$$ 2 \cdot b \cdot \sigma \cdot r - 2 \cdot b \cdot \sigma = b \cdot r \cdot \left( 1 + b + \sigma \right) + b \cdot \sigma \cdot \left( 1 + b + \sigma \right) $$$$ r \cdot \left( 2 \cdot \sigma - 1 - b - \sigma \right) = \sigma \cdot \left( 1 + b + \sigma + 2 \right) $$Finally $$ r_c = \frac{\sigma \cdot \left( \sigma + b + 3 \right)}{\sigma - b - 1} $$
Need to calculate $$ \frac{1}{V} \cdot \frac{dV}{dt} = ?, \qquad V = x \cdot y \cdot z $$
According to Liouville theorem $$ \frac{1}{V} \cdot \frac{dV}{dt} = \overline{\nabla} \cdot V = - \left( \sigma + 1 + b \right) $$
Solution $$ V\left( t \right) = \exp{\left\{V_0 - t\cdot\left( \sigma + 1 + b \right) \right\}} $$
This means that phase volume shrinks exponentially and the system is dissipative
In [1]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import pylab
pylab.rcParams['figure.figsize'] = (20.0, 20.0)
In [2]:
T = 250.
dt = 1E-2
STEPS = int(T/dt)
In [3]:
S = 10.
B = 8./3
# INITIAL = (2., -1., 0.)
INITIAL = (19.03492016970706, 22.626227200095837, 36.651063105952005)
In [11]:
R_C = S * (S + B + 3) / (S - B - 1)
R_L = R_C * .5
R_G = R_C * 2
In [5]:
def lorenz(x, y, z, r=R_C, s=S, b=B):
x_dot = s*(y - x)
y_dot = r*x - y - x*z
z_dot = x*y - b*z
return x_dot, y_dot, z_dot
In [6]:
# Need one more for the initial values
xs = np.empty((STEPS + 1,))
ys = np.empty((STEPS + 1,))
zs = np.empty((STEPS + 1,))
# Setting initial values
xs[0], ys[0], zs[0] = INITIAL
In [7]:
# Stepping through "time".
for i in range(STEPS):
# Derivatives of the X, Y, Z state
x_dot, y_dot, z_dot = lorenz(xs[i], ys[i], zs[i], r=R_C)
xs[i + 1] = xs[i] + (x_dot * dt)
ys[i + 1] = ys[i] + (y_dot * dt)
zs[i + 1] = zs[i] + (z_dot * dt)
In [8]:
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(xs, ys, zs)
ax.set_xlabel("X Axis")
ax.set_ylabel("Y Axis")
ax.set_zlabel("Z Axis")
ax.set_title("Lorenz Attractor")
plt.show()
In [9]:
plt.figure(1)
ax = plt.subplot(211)
plt.plot(xs, ys)
plt.title('XY phase portrait')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax = plt.subplot(212)
plt.plot(xs, zs)
plt.title('XZ phase portrait')
ax.set_xlabel('X')
ax.set_ylabel('Z')
plt.show()
In [10]:
INITIAL_VOLUME = INITIAL[0] * INITIAL[1] * INITIAL[2]
# INITIAL_VOLUME = 1
TIME = np.empty((STEPS+1,))
real_phase_volume = np.zeros((STEPS+1,), 'f')
time = 0
for t in range(STEPS+1):
time = t*dt
TIME[t] = time
real_phase_volume[t] = INITIAL_VOLUME*np.exp(-time*(S + 1 + B))
calculated_phase_volume = xs * ys * zs
plt.plot(TIME, calculated_phase_volume)
plt.plot(TIME, real_phase_volume)
plt.show()