3.6: Dynamics of Non-Linear Systems

Reading: Notes on Diffy Q’s Section 8.1, Section 8.2, and Section 8.3

Motion of a Pendulum

---

A mass \(m\) is attached to a point on the ceiling by a rod with length \(L\), and it is allowed to swing along a circlular path.

Image Credit: Chetvorno, Public domain, via Wikimedia Commons

If a pendulum is displaced from its equilibrium position, gravity will move the mass back towards its equilibrium position.

• Let \(\color{dodgerblue}{\theta}\) denote the angular position (in radians) of the mass.
  • Equilibrium position is \(\color{dodgerblue}{\theta =0}\).
• Let \(\color{tomato}{v = \dfrac{d\theta}{dt}}\) denote the angular velocity (radians per sec) of the mass.

Question 1:

---

For each of the initial conditions given below, answer the following questions.

• Describe how both the pendulum’s position and angular velocity change over time.
• Sketch a possible graph to illustrate the relation between \(\theta\) and \(v\) over time.
  • Plot values of \(\theta\) on the horizontal axis.
  • Plot values of \(v\) on the vertical axis.
• If you need a hint or want to check your work, see Hint for Question 1.

Question 1a:

---

Initial conditions \(\theta(0)=\frac{\pi}{6}\) and \(v(0)=0\).

Solution to Question 1a:

---

Question 1b:

---

Initial conditions \(\theta(0)= -\frac{5 \pi}{6}\) and \(v(0)=0\).

Solution to Question 1b:

---

Question 1c:

---

Initial conditions \(\theta(0)=0\) and \(v(0)=0.05\) radians per second (very small velocity).

Solution to Question 1c:

---

Question 1d:

---

Initial conditions \(\theta(0)=0\) and \(v(0)=50\) radians per second (very large velocity).

Solution to Question 1d:

---

Question 2:

---

How many equilibrium solutions does a pendulum have?

• Where are the equilibrium located in the \(\theta v\)-plane where we plotted solutions in Question 1?
• How would you classify the stablility of the equilibrium?

Solution to Question 2:

---

Equation of Motion for a Pendulum

---

There are two forces acting on bob of the pendulum:

1. The force due to gravity is given by \(F_g (\theta)= -mg\sin{\theta}\) where \(\color{dodgerblue}{m}\) denotes the mass of the bob and \(\color{mediumseagreen}{g}\) denotes the gravity constant.
2. The frictional force is given by \(F_f (\theta) = -b L \frac{d\theta}{dt}\), where \(\color{tomato}{b}\) denotes the damping coefficient and \(\color{mediumseagreen}{L}\) the length of the pendulum.

The force generated by the bob at the end of the pendulum is given by:

\[F = {\color{dodgerblue}{m}}a = {\color{dodgerblue}{m}} \frac{d^2s}{dt^2} = {\color{dodgerblue}{m}} \left( {\color{mediumseagreen}{L}} \frac{d^2 \theta}{dt^2} \right).\]

• \(s=\theta \cdot L\) denotes distance (arc length).

By balancing the forces, we get the following model for the motion of the simple pendulum:

\[\frac{d^2 \theta}{dt^2} + \frac{\color{tomato}{b}}{\color{dodgerblue}{m}} \frac{d\theta}{dt} + \frac{\color{mediumseagreen}{g}}{\color{mediumseagreen}{L}} \sin{\theta}=0.\]

• Below is a video deriving the equation of motion for the pendulum.

A Non-Linear Differential Equation

---

If we let \(v = \frac{d \theta}{dt}\) denote the velocity of the pendulum, set the length \(\color{mediumseagreen}{L = g= 9.8}\), mass \(\color{dodgerblue}{m=5}\) and damping coefficient \(\color{tomato}{b=1}\), the second-order differential equation above becomes

\[\frac{d^2 \theta}{dt^2} + 0.2 \frac{d \theta}{dt} + {\color{tomato}{\sin{\theta}}} = 0.\]

• Note this almost looks like a mass-spring system \(my'' + by' + ky = 0\).
• If the last term \(\color{tomato}{\sin{\theta}}\) was \(\theta\) instead, we would have the same dynamics.
• The \(\color{tomato}{\sin{\theta}}\) term means the differential equation is non-linear.

Expressing as a System of First Order Differential Equations

---

If we let \(\color{dodgerblue}{\dfrac{d \theta}{dt}=v}\), then we have \(\color{tomato}{\dfrac{d^2 \theta}{dt^2}=\dfrac{dv}{dt}}\), and the second order differential equation \({\color{tomato}{\dfrac{d^2 \theta}{dt^2}}} + 0.2 {\color{dodgerblue}{\dfrac{d \theta}{dt}}} + \sin{\theta} = 0\) can be rewritten as a first order differntial equation.

\[{\color{tomato}{\frac{dv}{dt}}} + 0.2 {\color{dodgerblue}{v}} + \sin{\theta} = 0.\]

Thus, the second order differential equation is equivalent to the non-linear system of two first order differential equations:

\[\begin{align} \color{dodgerblue}{\frac{d\theta}{dt}} & \color{dodgerblue}{=v} \\ {\color{tomato}{\frac{dv}{dt}}} &=-0.2{\color{dodgerblue}{v}}-\sin(\theta) \end{align}\]

Question 3:

---

The system of differential equations for the pendulum has how many equilibrium solutions?

• What are the equilibrium? Give answers as pairs \((\theta, v)\).
• Based on your intuition with pendulums, determine whether each equilibrium is stable, unstable, or semi-stable?
• How do your answers compare with your answers to Question 2?

\[\begin{align} \frac{d\theta}{dt} & =v \\ \frac{dv}{dt} &=-0.2v-\sin(\theta) \end{align}\]

Solution to Question 3:

---

Question 4:

---

You might recall that if \(\theta\) is small, \(\sin(\theta) \approx \theta\). Explain why this is true and then use this fact to approximate the above system with a linear system and classify the equilibrium solution at the origin.

Solution to Question 4:

---

Question 5:

---

If initially we push the mass with velocity \(v\) from its equilibrium position \(\theta=0\), approxiamte the range of initial velocities \(v\) that will result in the pendulum making exactly one complete rotation before eventually coming to rest.

Solution to Question 5:

---

If needed, change the initial conditions and view window for the phase plane plot!

# import the ode.tools module from github
!pip install git+https://github.com/CU-Denver-MathStats-OER/ODEs
from IPython.display import clear_output
clear_output()

from utils.ode_tools import plot_phase_sol  # Only need to import one time.
import numpy as np  # only need to import once

# Set viewing window

x = np.linspace(-10.0, 10.0, 21)  # theta is horizontal axis
y = np.linspace(-4.0, 4.0, 21)  # v is vertical axis

# Define system of odes
def f(Y, t):
    x, y = Y
    return [y , -0.2*y - np.sin(x)]  # enter f_1(theta, v) and f_2(theta ,v)

# Enter range of time
tspan = np.linspace(0, 50, 200) # range of time to visualize solution

# Enter initial values
x0 = 0  # initial value of theta is set to 0
y0 = 3  # initial value of v

# Plots a solution in phase plane
plot_phase_sol(x, y, f, tspan, x0, y0)

Question 6:

---

Based on the phase plane portraits you explored in Question 5, do you believe the equilibrium at \((\theta, v) = (\pi, 0)\) is stable, unstable, or semi-stable? Is this analysis consisten with your intuition?

Solution to Question 6:

---

Trust in Your Analysis!

---

Linearization and Linear Stability Analysis

---

We can perform linear stability analysis on a system of two or more variables, such as the one in the previous pendulum question. Consider a function \(f(x,y)\), then recall from multivariable calculus we find a formula for the linearization of \(f(x,y)\) around the point \((x_0,y_0)\) using the formula:

\[f(x,y) \approx L(x,y) = f(x_0,y_0) + f_x(x_0,y_0)({\color{dodgerblue}{x-x_0}}) + f_y(x_0,y_0)({\color{tomato}{y-y_0}}).\]

• Let \((x,y)=(x_0 + {\color{dodgerblue}{u}}, y_0 + {\color{tomato}{v}})\) denote a point close to \((x_0,y_0)\).
• Note that we have \({\color{dodgerblue}{u = x-x_0 = \Delta x}}\) and \({\color{tomato}{v = y - y_0 = \Delta y}}\).
• We can rewrite the linearization \(L(x,y)\) in terms of \(({\color{dodgerblue}{u}}, {\color{tomato}{v}})\).

\[f(x,y) \approx L({\color{dodgerblue}{u}}, {\color{tomato}{v}}) = f(x_0,y_0) + f_x(x_0,y_0){\color{dodgerblue}{u}} + f_y(x_0,y_0){\color{tomato}{v}}.\]

Linearizing a System at an Equilibrium

---

Let \((x_0, y_0)\) be an equilibrium of the system

\[\begin{align} \frac{dx}{dt} &= f(x,y) \\ \frac{dy}{dt} &= g(x,y) \end{align}.\]

• If we linearize at an equilibrium at \((x_0, y_0)\), then both \(f(x_0, y_0) =0\) and \(g(x_0, y_0) =0\).
• The resulting linearizations \(L_g\) and \(L_f\) for \(f\) and \(g\), respectively will be of the form

\[L_f(u,v) = {\color{dodgerblue}{f_x(x_0, y_0)}}u + {\color{tomato}{f_y(x_0, y_0)}}v = {\color{dodgerblue}{a}}u + {\color{tomato}{b}}v \qquad \mbox{and} \qquad L_g(u,v) = {\color{mediumseagreen}{g_x(x_0, y_0)}}u + {\color{mediumpurple}{g_y(x_0, y_0)}}v = {\color{mediumseagreen}{c}}u + {\color{mediumpurple}{d}}v\]

The linearized system of differential equations at equilibrium \((x_0,y_0)\) is therefore

\[\begin{align} \frac{dx}{dt} &= L_f(u,v) = {\color{dodgerblue}{a}}u+{\color{tomato}{b}}v\\ \frac{dy}{dt} &= L_g(u,v) = {\color{mediumseagreen}{c}}u + {\color{mediumpurple}{d}}v \end{align}.\]

Our derivatives are expressed with respect to the original variables \(x\) and \(y\), but the equations are now expressed with respect to variables \(u\) and \(v\). We have one last substitution to perform.

\[\mbox{if } u = x - x_0 \mbox{, then } \dfrac{du}{dt} = \dfrac{dx}{dt} \qquad \mbox{similarly} \qquad \mbox{if } v = y - y_0 \mbox{, then } \dfrac{dv}{dt} = \dfrac{dy}{dt}.\]

Thus, we obtain a linearized system of differential equations:

\[\begin{align} \frac{du}{dt} &= {\color{dodgerblue}{a}}u + {\color{tomato}{b}}v\\ \frac{dv}{dt} &= {\color{mediumseagreen}{c}}u + {\color{mediumpurple}{d}}v \end{align}.\]

• The behavior of the non-linear system near equilibrium \((x_0, y_0)\) is the same as the behavior of the linearized system near the origin \((0,0)\).
• The linearized system near \((x_0, y_0)\) cannot be used to analyze solutions that are not close \((x_0, y_0)\).

Expressing the Linearized System in Matrix Form

---

We can express the linearized system in matrix form

\[\begin{array}{l} & \dfrac{du}{dt} = {\color{dodgerblue}{a}}u + {\color{tomato}{b}}v\\ & \dfrac{dv}{dt} = {\color{mediumseagreen}{c}}u + {\color{mediumpurple}{d}}v \end{array} \qquad \longrightarrow \qquad \begin{bmatrix} \frac{du}{dt} \\ \frac{dv}{dt} \end{bmatrix} = \begin{bmatrix} {\color{dodgerblue}{a}} & {\color{tomato}{b}} \\ {\color{mediumseagreen}{c}} & {\color{mediumpurple}{d}} \end{bmatrix} \begin{bmatrix} u \\ v \end{bmatrix}.\]

The coefficients of the linearized system of differential equations are the partial derivatives of \(f\) and \(g\) evaluated at the equilibrium \((x_0, y_0)\):

\[{\color{dodgerblue}{a=f_x(x_0, y_0)}}, \qquad {\color{tomato}{b=f_y(x_0, y_0)}}, \qquad {\color{mediumseagreen}{c=g_x(x_0, y_0)}}, \qquad {\color{mediumpurple}{d=g_y(x_0, y_0)}}.\]

Therefore, the linearized system of differential equations can be written as

\[\begin{bmatrix} \frac{du}{dt} \\ \frac{dv}{dt} \end{bmatrix} = \begin{bmatrix} {\color{dodgerblue}{f_x(x_0, y_0)}} & {\color{tomato}{f_y(x_0, y_0)}} \\ {\color{mediumseagreen}{g_x(x_0, y_0)}} & {\color{mediumpurple}{g_y(x_0, y_0)}} \end{bmatrix} \begin{bmatrix} u \\ v \end{bmatrix}\]

The Jacobian Matrix

---

The matrix

\[J_{(x_0, y_0)} = \begin{bmatrix} {\color{dodgerblue}{f_x(x_0, y_0)}} & {\color{tomato}{f_y(x_0, y_0)}} \\ {\color{mediumseagreen}{g_x(x_0, y_0)}} & {\color{mediumpurple}{g_y(x_0, y_0)}} \end{bmatrix}\]

is called the Jacobian matrix.

Stability Analysis for Non-Linear Systems

---

Consider a non-linear system of first order differential equations

\[\begin{align} \frac{dx}{dt} &= f(x,y) \\ \frac{dy}{dt} &= g(x,y) \end{align}.\]

Step 1: Find all equilibrium points.

Step 2: Find expressions for Partial Derivatives \(f_x\), \(f_y\), \(g_x\), and \(g_y\).

Step 3: Find the Jacobian matrix for a selected equilibrium \((x_0, y_0)\).

Step 4: Find eigenvalues and eigenvectors for the Jacobian to sketch the phase plane near \((x_0, y_0)\).

Repeat for all equilibrium if desired.

Question 7:

---

Consider the system of differential equations

\[\begin{align} \frac{dx}{dt} &= 1-x^2 \\ \frac{dy}{dt} &= -3x -3y. \end{align}\]

Note the system has two equilibria, at \((1,-1)\) and \((-1,1)\). Linearize the system at \((1,-1)\), and determine whether the equilibrium is stable, unstable, or semi-stable.

Solution to Question 7:

---

Question 8:

---

Consider the system of differential equations

\[\begin{align} \frac{dx}{dt} &= 1-x^2 \\ \frac{dy}{dt} &= -3x -3y. \end{align}\]

Linearize the system at the other equilibrium \((-1,1)\), and determine whether that equilibrium is stable, unstable, or semi-stable.

Solution to Question 8:

---

Question 9:

---

Consider again

\[\begin{align} \frac{dx}{dt} &= 1-x^2 \\ \frac{dy}{dt} &= -3x -3y \end{align}\]

Combine your results from Question 7 and Question 8, to sketch a possible phase plane for the system of differential equations.

• Does an analysis of the system using nullclines corroborate your linear stability analysis?
• Sketch the phase plane portrait using the phase_portrait function.

Solution to Question 9:

---

from utils.ode_tools import phase_portrait  # only need to import one time
#import numpy as np  # only need to import one time

# Set viewing window
x = np.linspace(-5.0, 5.0, 21)  
y = np.linspace(-5.0, 5.0, 21)

def f(Y, t):
    x, y = Y
    return [??,
            ??]  

# Plots a phase portrait
phase_portrait(x, y, f)

Question 10:

---

Use linear stability analysis to classify the critical points you found in the pendulum system.

\[\begin{align} \frac{d\theta}{dt} &=v \\ \frac{dv}{dt} &= -0.2v - \sin(\theta) \end{align}\]

Solution to Question 10:

---

Question 11:

---

Use linear stability analysis to classify each of the critical points for the system of differential equations below.

\[\begin{align} \dfrac{dx}{dt} &= 4x-2xy\\ \dfrac{dy}{dt} &= 6y+2xy \end{align}\]

Solution to Question 11:

---

Appendix I: Plotting with phase_portrait and plot_phase_sol

---

In the code cells below:

• First load the ode_tools module from GitHub.
• Then from the ode_tools module we import phase_portrait and/or plot_phase_sol.

Loading ode_tools from GitHub

---

• Run the code cell below to load the most up to date modules stored in GitHub.
• You will only need to run this code cell one time during an active session.

!pip install git+https://github.com/CU-Denver-MathStats-OER/ODEs
from IPython.display import clear_output
clear_output()

Plot Phase Portrait with phase_portrait

---

Importing the phase_portrait Plotting Function

---

After you run the code cell above to load the ode_tools module from GitHub you are now ready to import the phase_portrait function.

from utils.ode_tools import phase_portrait  # Only need to import one time.

Defining the System of Differential Equations

---

We use \(x\) and \(y\) as the generic symbols for the two dependent variables.

\[\begin{align} \frac{dx}{dt} &= f_1(x, y)\\ \frac{dy}{dt} &= f_2(x, y) \end{align}\]

Below, we enter the system we analyzed in Question 3.

\[\begin{aligned} \frac{d\theta}{dt} &=v \\ \frac{dv}{dt} &= -0.2v - \sin{\theta} \end{aligned}\]

import numpy as np

# Set viewing window

x = np.linspace(-10.0, 10.0, 21)  # theta is horizontal axis
y = np.linspace(-2.0, 2.0, 21)  # v is vertical axis

def f(Y, t):
    x, y = Y
    return [y , -0.2*y - np.sin(x)]  # enter f_1(theta, v) and f_2(theta ,v)

Plotting with phase_portrait

---

# Plots a phase portrait
phase_portrait(x, y, f)

Plot Phase Portrait with plot_phase_sol

---

Importing the plot_phase_sol Plotting Function

---

After you have loaded the ode_tools module from GitHub you are now ready to import the phase_portrait_sol function.

from utils.ode_tools import plot_phase_sol  # Only need to import one time.

Defining the System of Differential Equations and Initial Conditions

---

Below, we enter the system we analyzed in Question 1a with initial conditions \(\theta(0) = \frac{\pi}{6}\) and \(v(0)=0\).

\[\begin{aligned} \frac{d\theta}{dt} &= v \\ \frac{dv}{dt} &=-0.2v - \sin{\theta} \end{aligned}\]

import numpy as np

# Set viewing window

x = np.linspace(-3, 3, 21)  # theta is horizontal axis
y = np.linspace(-1.0, 1.0, 21)  # v is vertical axis

def f(Y, t):
    x, y = Y
    return [y , -0.2*y - np.sin(x)]  # enter f_1(theta, v) and f_2(theta ,v)

# Enter range of time
tspan = np.linspace(0, 50, 200) # range of time to visualize solution

# Enter initial values
x0 = np.pi / 6  # initial value of theta
y0 = 0  # initial value of v

Plotting with phase_portrait_sol

---

# Plots a phase portrait
plot_phase_sol(x, y, f, tspan, x0, y0)

Appendix II: Hint for Question 1

---

Image credit: Jacopo Bertolotti, CC0, via Wikimedia Commons

Creative Commons License Information

---

Exploring Differential Equations by Adam Spiegler is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License.
Based on a work at https://github.com/CU-Denver-MathStats-OER/ODEs and original content created by Rasmussen, C., Keene, K. A., Dunmyre, J., & Fortune, N. (2018). Inquiry oriented differential equations: Course materials. Available at https://iode.sdsu.edu.