System of differential equations for a pendulum is derived using Newtonian and Lagrangian mechanics. The system is then linearized about three different equilibium points. Transfer function and state space model are developed for the linearized systems. Simulation code is also provided.

The pendulum is made of a point mass attached to the end of a massless bar. The other end of the bar is attached to a fixed ground. The mass and the bar rotate about this fixed point. Gravity acts downwards.

The point mass has mass *m* and the bar is of length *l*.
The angle the bar makes
with the downward direction is *θ*. When the mass is at the bottom-most
position, *θ* = 0 and *θ* increases as the pendulum rotates
counter-clockwise. There is also air resistance which acts in the direction
opposite to the velocity and whose magnitude scales with the speed by *μ*.
An external torque *τ* is applied to this pendulum. *τ* is positive
counter-clockwise.

Variables and constants of concern are listed and organized as follows:

(input) | τ(t) | = | external torque applied to pendulum | |

(output) | θ(t) | = | angle of pendulum | unknown |

m | = | mass | ||

l | = | arm length | ||

μ | = | air friction coefficient |

Differential equation describing the dynamics of the pendulum is derived using Newtonian mechanics. Since there is only one unknown, one equation is needed.

Torques around the axis of rotation can be visualized using the following free body diagram:

Using Newton's 2nd law:

where *I* is the inertia and is *ml*^{2} for a point mass
rotating *l* away from the center of rotation.
*α* is angular acceleration.

Substituting all of the variables:

Solving for angular acceleration, the differential equation describing the motion of the pendulum can be written as:

The above equation has sin*θ* so the system is nonlinear.

Differential equation describing the dynamics of the pendlum is derived using Lagrangian mechanics.

Kinetic energy is

where speed *v* = *l**θdot*.

Potential energy is

where the height *h* is 0 at the axis of rotation and positive upwards
so *h* = -*l*cos*θ*.

Velocity-proportional friction forces can be accounted in Lagrangian mechanics using Rayleigh dissipation function:

The lagrangian is

and the equation of motion with external force is solved with Lagrange's equation,

This is solved using Python's symbolic math library:

from sympy import * t, m, g, l, mu, tau = symbols('t m g l mu tau') theta = symbols(r'\theta', cls=Function) theta = theta(t) theta_d = diff(theta, t) theta_dd = diff(theta_d, t) v = l * theta_d h = -l * cos(theta) T = 1/2 * m * v ** 2 V = m * g * h L = T - V D = 1/2 * mu * theta_d ** 2 # Lagrangian LE = diff(diff(L, theta_d), t) - diff(L, theta) + diff(D, theta_d) - tau # Set equal to zero and solve for theta_dd rhs = solve(LE, theta_dd) pprint(rhs)

Output is:

which is identical to the result obtained with Newton's method.

The equation of motion derived in the previous sections is a second order differential equation. This is converted to a system of first order differential equations by letting

The state space model is then,

Equilibrium points are states vec*x _{eq}* and
input

This refers to the state when the pendulum is not moving nor accelerating. Letting

then

Doing Taylor expansion about vec*δx*=0 and *δu*=0 and ignoring second
and higher order terms, the nonlinear state space model can be linearized as

where *A*, *B*, *C*, *D* matrices are jacobians evaluated
at the equilibrium point:

In the case of the pendulum, the jacobians are:

Three arbitrary equilibrium points are selected for analysis:

This occurs when *θ* = 0 and *τ* = 0 (when no external torque is applied).
Letting

The linearized state space at this equilibrium is

The eigenvalues of the system matrix *A* provides insights into the stability of
the system. Analytical solution is given as:

If we let *g* = 9.81, *m* = 0.5, *l* = 0.2, *μ* = 0.15,
the eigenvalues of the
system matrix *A* are -3.75 ± 5.91502i. Since all real parts are negative,
the system is stable. A phase plane also shows that this is a stable node:

Next, it is possible to compute the transfer function using the following formula:

Solving this gives,

The above can be converted to differential equation as

which is a simple harmonic oscillator if there is no friction (*μ* = 0) and
no external torque (*τ* = 0).

The above calculations were done using Python's symbolic math library:

from sympy import * m, g, l, mu, u, s, t = symbols('m g l mu u s t') x = MatrixSymbol('x', 2, 1) F = Matrix([ [x[1]], [-g/l * sin(x[0]) - mu/(m * l**2) * x[1] + 1/(m * l**2) * u] ]) G = Matrix([ [x[0]] ]) x_eq = Matrix([ [0], # 0 or pi [0] ]) u_eq = 0 A = F.jacobian(x).subs(x, x_eq) B = F.jacobian([u]).subs(u, u_eq) C = G.jacobian(x).subs(x, x_eq) D = G.jacobian([u]).subs(u, u_eq) transfer_function = C * (s * eye(2) - A) ** -1 * B + D A_subs = A.subs([(g, 9.81), (m, 0.5), (l, 0.2), (mu, 0.15)]) pprint(A) pprint(B) pprint(C) pprint(D) pprint(A.eigenvals()) pprint(A_subs) pprint(A_subs.eigenvals()) pprint(transfer_function)

Phase portrait was plotted as below:

import numpy as np import matplotlib.pyplot as plt v = 3 w = np.pi/4 Y, X = np.mgrid[-v:v:100j, -w:w:100j] A = np.array([ [0, 1], [-49.05, -7.5] ]) s, v = np.linalg.eig(A) U = A[0,0]*X + A[0,1]*Y V = A[1,0]*X + A[1,1]*Y if s.imag[0] == 0 and s.imag[1] == 0: t = np.arange(-w, w, 0.01) plt.plot(t, (v[1,0]/v[0,0])*t) plt.plot(t, (v[1,1]/v[0,1])*t) plt.streamplot(X, Y, U, V) plt.xlabel(r'$\delta \theta$') plt.ylabel(r'$\delta \dot{\theta}$') plt.show()

This occurs when *θ* = *π* and *τ* = 0
(when no external torque is applied). Letting

The linearized state space at this equilibium is

The eigenvalues of the system matrix *A* is:

If we let *g* = 9.81, *m* = 0.5, *l* = 0.2, *μ* = 0.15,
the eigenvalues of *A* are
-11.69434 and 4.19434. Since there exists an eigenvalue with positive real part,
it is unstable. There is one positive and one negative eigenvalue so this is a
saddle point and this can also be seen from the phase portrait:

Transfer function is,

Differential equation is,

The pendulum can be in equilibium in a tilted position if there is an external torque holding it in place. This happens at:

The external torque that just balances the torque produced by the gravity can be obtained using a free body diagram:

Another way to obtain the torque at this equilibium is by solving it algebraically:

It is also possible to find this equilibium point by using numerical method.
If we let *g* = 9.81, *m* = 0.5, *l* = 0.2, *μ* = 0.15,
an optimization solver is used to find the equilibium torque:

from scipy import optimize import numpy as np # Constants g = 9.81 l = 0.2 m = 0.5 mu = 0.15 # Desired angle at equilibrium theta_eq = np.pi/4 # Nonlinear pendulum differential equation: xdot = f(x, u) def f(x, u): return [ x[1], -g / l * np.sin(x[0]) - mu/(m * l**2) * x[1] + 1/(m * l**2) * u ] # Optimization function def fun(z): x = z[:2] u = z[2] xdot = f(x, u) theta_err = x[0] - theta_eq return [ xdot[0], xdot[1], theta_err ] # Initial guess x0 = np.array([theta_eq, 0]) u0 = np.array([0]) # Find equilibrium z0 = np.concatenate((x0, u0)) result = optimize.root(fun, z0) print(result.x[:2]) # x at equilibium print(result.x[2]) # u at equilibrium

Above code returns the equilibium torque to be 0.69367 which is identical to the analytical solution:

Now, linear state space at this equilibium is determined. Letting

The linearized state space at this equilibium is

The linearized dynamics is only valid for small angles around the equilibium point. The full nonlinear dynamics can be simulated using numerical integration. Python's ODE solver is used:

import numpy as np from scipy.integrate import solve_ivp import matplotlib.pyplot as plt # Constants g = 9.81 l = 0.2 m = 0.5 mu = 0.02 def xdot(t, x): u = 0 # No input torque return [ x[1], -g / l * np.sin(x[0]) - mu/(m * l**2) * x[1] + 1/(m * l**2) * u ] x0 = [0, 17] # Initial state [theta0, theta_dot0] t0 = 0 tf = 6 dt = 0.01 # Simulate sol = solve_ivp(xdot, [t0, tf], x0, t_eval=np.arange(t0, tf, dt)) # Plot as function of time fig, axs = plt.subplots(2) axs[0].plot(sol.t, sol.y[0]) axs[0].set_ylabel(r'$\theta$') axs[0].grid() axs[1].plot(sol.t, sol.y[1]) axs[1].set_ylabel(r'$\dot{\theta}$') axs[1].grid() plt.show()

Output:

In this simulation, the pendulum is initially positioned at the bottom with a counter-clockwise angular speed. The pendulum passes through the top and comes back down. After oscillating a little, it slowly comes to rest due to air resistance. There is no input torque. The animation at the top of this page was created with this method.

The following versions were used to run the code:

- Python 3.11.1
- sympy 1.12
- numpy 1.26.4
- scipy 1.10.1
- matplotlib 3.7.1