A dynamical system can be expressed in the form
$$ \frac{\mathrm{d}}{\mathrm{d}t} \phi(t)= \mathcal{F}\big(\phi(t),\,t\big) \tag{1}. $$Above, $\phi(t)$ includes all variables needed to describe the state of the system. (Hereafter, dots above variables denote time-differentiation.)
For a pendulum that evolves freely in three dimensions, $\phi$ includes the position and the velocity of its free end. In this case, $\phi$ has six components: $$ \phi = [x(t), y(t), z(t), \dot{x}(t), \dot{y}(t), \dot{z}(t)]. $$
If $\mathcal{F}$ in Eq. (1) does not depend explicitly on time $t$ we say that the dynamical system is autonomous. Otherwise, the dynamical system is non-autonomous.
Given equation (1) and the state of the system at some instance $t=t_0$, e.g., $\phi(t_0)=\phi_0$, we can predict the state of the system $\phi(t)$ for all times $t$.
Any $n$-th order differential equation can be rewritten as a system of $n$ differential equations all being first order. In that sense, all dynamical systems are first-order in time (but of course their dimensionality may vary).
Show that any non-autonomous system of $n$ differential equations corresponds to an autonomous system of $n+1$ differential equations.
Hint: You can start by an example: Show that the forced 1D harmonic oscillator
\begin{align*} \ddot{x} + 2\gamma\dot{x}+\omega^2 x = f_0\cos(2\pi t), \end{align*}can be rewritten as a system of 3 equations that do not depend explicitly on time. Then generalise this to what the exercise wants.
Hint: You may need to change/introduce a new variable.
Show that any nonlinear dynamical system can be rewritten as a linear dynamical system of higher dimensionality.
Take as an example the system:
\begin{align*} \dot{x} = -x^2. \end{align*}How big is the dimension of the equivalent linear system for the example above?
Food for thought: If any nonlinear system can be transformed into a linear one then what's the point of distinguish between linear and nonlinear dynamical systems?
Fixed points or equilibrium points are special values of $\phi=\phi^e$ which satisfy:
$$ \mathcal{F}(\phi^e) = 0. \tag{2} $$They are called equilibria or fixed points because if the system finds itself in one of these states then it will stay there for eternity.
The character of a fixed point is studied by assuming small deviations about that fixed point,
$$ \phi = \phi^e + \phi',\quad\text{where}\quad\|\phi'\|\ll 1. \tag{3} $$Determining how these small deviations will evolve characterises the fixed point. We do that by writing out the linear dynamical system that governs $\phi'$. By inserting Eq. (3) in Eq. (1) and after linearization (i.e., discard of the terms that involve quantites which are quadratic or higher order in components of $\phi'$), we get:
$$ \dot{\phi'} = \mathcal{F}(\phi^e + \phi') \approx \mathcal{A}(\phi^e)\,\phi',\tag{4} $$where $\mathcal{A}$ is a linear operator that depends on the fixed point $\phi^e$.
Consider a pendulum in one dimension. We learn that the angle $\theta$ from rest is governed by
$$ \ddot{\theta} + \frac{g}{\ell} \sin{\theta} = 0. \tag{5} $$Note that this is a second-order dynamical law. But if we define $\omega(t) \equiv \dot{\theta}$ and $\phi=[ \theta, \omega ]^\mathrm{T}$ then we can rewrite the second-order equation of motion as:
$$ \begin{bmatrix}\dot{\theta}(t)\\ \dot{\omega}(t) \end{bmatrix} = \begin{bmatrix}\omega(t)\\-\frac{g}{\ell} \sin[\theta(t)]\end{bmatrix}. \tag{6} $$What are the pendulum's fixed points? $$ \phi^e_1 = \begin{bmatrix}0\\0\end{bmatrix}\quad\textrm{and}\quad \phi^e_2 = \begin{bmatrix}\pi\\0\end{bmatrix} \tag{7} $$
We say that $\phi^e$ is *attracting* or *asymptotically stable* (or *modally stable*) if all trajectories that start near $\phi^e$ approach it as $t\to \infty$.
In other words, any trajectory $\phi^e(t)$ that starts within a "distance" $\delta$ of $\phi^e$ is guaranteed to converge to $\phi^e$ *eventually*.
Asymptotic stability concerns what happens at $t\to\infty$. Eigenanalysis of a dynamical system reveals the $t\to\infty$ behavior of it.
How do we study the asymptotic/modal stability of these fixed point?
where operator $\mathcal{A}$ depends on the equilibrium state $\phi^e$. 4. Search for eigensolutions $\phi' = \widetilde{\phi}\,e^{\sigma t}$. This transforms (6) in an eigenvalue problem with $\sigma$ the eigenvalue and $\widetilde{\phi}$ the eigenvector/eigenfunction: $$ \mathcal{A}\,\widetilde{\phi} = \sigma \, \widetilde{\phi}. \tag{9} $$ 5. Eigenanalysis of $\mathcal{A}$ determines its spectrum (spectrum is the set of all eigenvalues). 6. If we have at least one eigenvalue $\sigma$ with $\mathrm{Re}(\sigma)>0$ then we have modal instability. Otherwise the equilibrium $\phi^e$ is modally stable.
A bit more on point 6: We say that if we have an eigenvalue $\sigma$ with $\mathrm{Re}(\sigma)>0$ then we have modal instability. That is because if we have a perturbation which initially looks like the eigenfunction of that fixed point, $\phi'(t=0) = \widetilde{\phi}$, then that perturbation will grow away from $\phi^e$ at an exponential rate.
Confirm the above statement and then generalise it. What happens in the above system if initially we start of with any random perturbation (not with the "special" perturbation $\phi'(t=0) = \widetilde{\phi}$)? Will the system again move away from $\phi^e$ at an exponential rate?
For our pendulum example we have that for $\phi^e_1$ \begin{align*} \mathcal{A}_1 = \begin{pmatrix}0 & 1 \\ -\frac{g}{\ell} & 0\end{pmatrix} \end{align*} while for $\phi^e_2$: \begin{align*} \mathcal{A}_2 = \begin{pmatrix}0 & 1 \\ +\frac{g}{\ell} & 0\end{pmatrix} \end{align*}
Eigenanalysis of $\mathcal{A}_1$ gives that its eigenvalues are $\sigma = +i\sqrt{\frac{g}{\ell}}, -i\sqrt{\frac{g}{\ell}}$.
Eigenanalysis of $\mathcal{A}_2$ gives that its eigenvalues are $\sigma = +\sqrt{\frac{g}{\ell}}, -\sqrt{\frac{g}{\ell}}$.
None of the eigenvalues of $\mathcal{A}_1$ have positive real part $\Rightarrow$ $\phi^e_1$ is modally stable.
There is exist one eigenvalue of $\mathcal{A}_2$ with $\mathrm{Re}(\sigma)>0$ $\Rightarrow$ $\phi^e_2$ is modally unstable.
Fill in the missing dots in the pendulum example. First make sure that you can reproduce the above. Then proceed finding the eigenvectors that correspond to those eigenvalues. How would the fastest growing mode look like?
A take-away message from the lectures is that actually obtaining analytic solutions for a dynamical system is tedious (sometimes even impossible) and, furthermore, it provides very little insight.
Instead, one can get a very good idea of how the dynamical system will behave by drawing the so-called phase portrait.
The phase portrait is constructed by finding attractors (e.g., fixed points or limit cycles) and the corresponding stability characteristic of these attractors. Oftentimes it helps to draw the nullclines (locus of points where each component of $\dot{\phi}$ vanishes).
In most cases finding analytical solutions in closed form for a dynamical system is hard or impossible. That is, we want a method for constructing approximations of the solution for:
\begin{align*} \dot{\phi}(t) &= \mathcal{F}\big(\phi(t)\,,\, t\big), \tag{10a}\\ \phi(t_0) & =\phi_0. \tag{10b} \end{align*}How do we proceed?
Let $\phi(t)$ the solution of Eqs. (10) we want to approximate. Now use a Taylor-series expansion about $t=t_0$:
\begin{align*} \phi(t) &= \phi(t_0) + (t-t_0)\dot{\phi}\big|_{t=t_0} + \tfrac1{2!}(t-t_0)^2\ddot{\phi}\big|_{t=t_0} + \dots \\ &= \phi(t_0) + (t-t_0) \mathcal{F}(\phi_0,t_0) + \tfrac1{2!}(t-t_0)^2 \frac{\mathrm{d}\mathcal{F}}{\mathrm{d}t}\big|_{t=t_0} + \dots \tag{11}\\ \end{align*}If $t$ is close to $t_0$ then we can neglect terms of order $(t-t_0)^2$ and higher. This gives the so-called Euler approximation.
Assume a time-discretization with time-step $\delta$, i.e.,
$$ t_n = t_0 + n\,\delta,\ n=0,1,\dots. \tag{12} $$Let's denote the values of the approximations of $\phi(t)$ at times above as $\phi_n \equiv \phi(t_0+n\delta)$. The Forward Euler approximation is then
$$ \phi_{n+1} = \phi_n + \delta\,\mathcal{F}\big(\phi_n, t_n\big). \tag{13} $$Let's try to integrate
\begin{align*} \dot{x} &= -x + \sin t, \\ x(0) &= x_0. \end{align*}In this case, we could actually find a closed-form solution of the system:
$$ x(t) = e^{-t}\,x_0 + \frac1{2}(e^{-t} + \sin t - \cos t). $$IMPORTANT: Always start with an example that you know what to expect!
def RHS(x, t):
xdot = -x + sin(t)
return xdot
def Euler(xn, tn, dt):
xnext = xn + dt*RHS(xn, tn)
return xnext
x0 = 2.0 # initial position
tfin, nsteps = 12.0, 25 # final time, number of time-steps
t = np.linspace(0, tfin, nsteps); dt = t[1]-t[0]
x_euler = 0*t; x_euler[0] = x0 #initialize array and set initial condition
for j in np.arange(0, nsteps-1):
xn, tn = x_euler[j], t[j]
x_euler[j+1] = Euler(xn, tn, dt)
t_an = np.linspace(0, 12.0, 100*nsteps) # a very dense time array
x_an = x0*exp(-t_an) + 0.5*(exp(-t_an) + sin(t_an) - cos(t_an))
plt.figure(figsize=(12, 4))
plt.plot(t_an, x_an, label="truth", color="dodgerblue")
plt.plot(t, x_euler, '^', label="Euler", color="darkgreen")
plt.legend(fontsize=20); plt.xlabel(r"$t$"); plt.ylabel(r"$x(t)$");
In Forward-Euler timestepping scheme we ignore all terms of $\mathcal{O}(\delta^2)$.
@ every time-step $\longrightarrow$ error $\sim \delta^2$
$\dfrac{t_{\rm final}}{\delta}$ time-steps $\longrightarrow$ total error $\sim \dfrac{t_{\rm final}}{\delta} \delta^2 \sim \delta$
The Runge-Kutta 4th order time-stepping scheme is: $$ \phi_{n+1} = \phi_n + \dfrac{\delta}{6}\big(k_1 + 2k_2 + 2k_3 + k4\big), \tag{14} $$ where \begin{align} k_1 & = \mathcal{F}\big(\phi_n\,,\, t_n\big), \tag{15a}\\ k_2 & = \mathcal{F}\big(\phi_n + \tfrac1{2} k_1 \delta \,,\, t_n + \tfrac1{2}\delta\big), \tag{15b}\\ k_3 & = \mathcal{F}\big(\phi_n + \tfrac1{2} k_2 \delta \,,\, t_n + \tfrac1{2}\delta\big), \tag{15c}\\ k_4 & = \mathcal{F}\big(\phi_n + k_3 \delta \,,\, t_n + \delta\big). \tag{15d}\\ \end{align}
RK4 @ every time-step $\longrightarrow$ error $\sim \delta^5$
total error $\longrightarrow$ $\mathcal{O}(\delta^4)$
def RK4(xn, tn, dt):
k1 = RHS(xn, tn)
k2 = RHS(xn+k1*dt/2, tn+dt/2)
k3 = RHS(xn+k2*dt/2, tn+dt/2)
k4 = RHS(xn+k3*dt , tn+dt)
xnext = xn + (dt/6)*(k1+2*k2+2*k3+k4)
return xnext
x_rk4 = 0*t; x_rk4[0] = x0 #initialize array and set initial condition
for j in np.arange(0, nsteps-1):
xn, tn = x_rk4[j], t[j]
x_rk4[j+1] = RK4(xn, tn, dt)
plt.figure(figsize=(12, 4))
plt.plot(t_an, x_an, label="truth", color="dodgerblue")
plt.plot(t, x_euler, '^', label="Euler", color="darkgreen")
plt.plot(t, x_rk4, 's', label="RK4", color="salmon")
plt.legend(fontsize=20); plt.xlabel(r"$t$"); plt.ylabel(r"$x(t)$");
Consider the following:
\begin{align*} \dot{x} &= y,\\ \dot{y} &= -x,\\ x(0)&=x_0,\\ y(0)&=y_0. \end{align*}(What's do these equations describe?)
def RHS(x, t):
xdot = np.zeros((2,))
xdot[0] = x[1]
xdot[1] = -x[0]
return xdot
The solution is
\begin{align*} x(t) &= x_0 \cos t + y_0 \sin t,\\ y(t) &= -x_0 \sin t + y_0 \cos t. \end{align*}x0 = np.zeros((2,))
x0[0], x0[1] = 1.0, 0.0
tfin, nsteps = 3*pi, 40 # final time, number of time-steps
t = np.linspace(0, tfin, nsteps); dt = t[1]-t[0]
x_euler = np.zeros((2, len(t))); x_euler[:, 0] = x0
x_rk4 = np.zeros((2, len(t))); x_rk4[:, 0] = x0
for j in np.arange(0, nsteps-1):
xn, tn = x_euler[:, j], t[j]
x_euler[:, j+1] = Euler(xn, tn, dt)
xn, tn = x_rk4[:, j], t[j]
x_rk4[:, j+1] = RK4(xn, tn, dt)
t_an = np.linspace(0, tfin, 100*nsteps)
x_an = np.zeros((2, len(t_an)))
x_an[0, :], x_an[1, :] = cos(t_an), -sin(t_an)
plt.figure(figsize=(12, 4))
plt.plot(t_an, x_an[0, :], label="truth", color="dodgerblue")
plt.plot(t, x_rk4[0, :], 's', label="RK4", color="salmon")
plt.plot(t, x_euler[0, :], '^', label="Euler", color="darkgreen")
plt.legend(fontsize=20)
plt.xlabel(r"$t$"); plt.ylabel(r"$x(t)$");
What do you expect to see?
plt.figure(figsize=(6, 6))
x, xdot = x_an[0, :], x_an[1, :]
plt.plot(x, xdot, label="truth", color="dodgerblue")
x, xdot = x_rk4[0, :], x_rk4[1, :]
plt.plot(x, xdot, "s--", label="RK4", color="salmon")
x, xdot = x_euler[0, :], x_euler[1, :]
plt.plot(x, xdot, "^-.", label="Euler", color="darkgreen")
plt.xlabel(r"$x$"); plt.ylabel(r"$\dot{x}$");
plt.axis('square');
Why Euler method systematically spirals out?
Consider the damped harmonic oscillator:
$$ \ddot{x} + 2\gamma \dot{x} + x = 0,\;\; x(0)=x_0,\;\; \dot{x}(0)=\upsilon_0. $$Take $\gamma=0.1$ and use numerical parameters as above tfin=3pi
and nsteps=40
. Integrate using Euler and RK4. What do you conclude about the stability character of the fixed point $x=\dot{x}=0$ just from the results of numerical integrations?
How can you make sure that any conclusions you deduce from numerics are actually what the physics imply and not a nuance of the numerics?