In earlier chapters, many models were obtained directly from force balances, circuit laws, or free-body diagrams. The Lagrangian approach provides a unified alternative: instead of summing forces component-by-component in Cartesian coordinates, we choose physically meaningful coordinates \( q_1,\dots,q_n \), compute kinetic and potential energies, and derive equations from an energy balance-compatible principle.

The core object is the \( Lagrangian \) \( L(q,\dot{q},t)=T(q,\dot{q},t)-V(q,t) \), where \( T \) is kinetic energy and \( V \) is potential energy. For dissipative systems we also use a Rayleigh dissipation function \( \mathcal{D}(\dot{q}) \).

\[ \frac{d}{dt}\left(\frac{\partial L}{\partial \dot{q}_i}\right) - \frac{\partial L}{\partial q_i} + \frac{\partial \mathcal{D} }{\partial \dot{q}_i} = Q_i, \qquad i=1,\dots,n \]

Here \( Q_i \) is the non-conservative generalized force associated with coordinate \( q_i \). This form is especially useful for systems with coupled translational and rotational motion, geometric constraints, and nonlinear coordinate mappings.

flowchart TD
  A["Physical system and assumptions"] --> B["Choose generalized coordinates q"]
  B --> C["Kinematics: positions and velocities"]
  C --> D["Compute energies T and V"]
  D --> E["Add dissipation D and generalized forces Q"]
  E --> F["Apply Euler-Lagrange equations"]
  F --> G["Obtain nonlinear ODEs / state-space model"]
  G --> H["Simulate, linearize, validate"]
        

A mechanical system with many particle coordinates can often be described by a smaller set of independent variables called \( generalized\ coordinates \). If the configuration can be parameterized by \( q \in \mathbb{R}^n \), then each particle position is a smooth function of \( q \) and (possibly) time:

\[ \mathbf{r}_k = \mathbf{r}_k(q,t), \qquad \dot{\mathbf{r} }_k = \sum_{i=1}^{n}\frac{\partial \mathbf{r}_k}{\partial q_i}\dot{q}_i + \frac{\partial \mathbf{r}_k}{\partial t}. \]

If holonomic constraints are present, they can be written as \( \phi_\alpha(q,t)=0 \) for \( \alpha=1,\dots,m \). Using an independent coordinate set removes these constraints explicitly. If we retain redundant coordinates, constraint forces appear through Lagrange multipliers.

\[ \frac{d}{dt}\left(\frac{\partial L}{\partial \dot{q} }\right) - \frac{\partial L}{\partial q} = Q + A(q,t)^T\lambda, \qquad \phi(q,t)=0, \]

where \( A(q,t)=\partial \phi / \partial q \) is the constraint Jacobian and \( \lambda \) collects the multiplier variables. This matrix form is important in multi-body and constrained mechatronic systems.

Let \( \delta q \) be a virtual displacement compatible with the constraints at a fixed time. If a force \( \mathbf{F}_k \) acts on particle \( k \), the virtual work is

\[ \delta W = \sum_{k}\mathbf{F}_k \cdot \delta \mathbf{r}_k = \sum_{k}\mathbf{F}_k \cdot \left(\sum_{i=1}^{n} \frac{\partial \mathbf{r}_k}{\partial q_i}\delta q_i\right) = \sum_{i=1}^{n} Q_i \delta q_i, \]

which defines the generalized force components:

\[ Q_i = \sum_k \mathbf{F}_k \cdot \frac{\partial \mathbf{r}_k}{\partial q_i}. \]

For conservative forces derived from a potential, generalized forces are already represented in \( V(q) \). Non-conservative contributions (actuator forces, viscous damping, friction approximations, torques, and inputs) should be added through \( Q_i \) and/or \( \mathcal{D}(\dot{q}) \).

If the dissipative force is linear viscous and can be represented by a positive semidefinite matrix \( R \), we may choose

\[ \mathcal{D}(\dot{q})=\frac{1}{2}\dot{q}^T R \dot{q}, \qquad \frac{\partial \mathcal{D} }{\partial \dot{q} } = R\dot{q}. \]

We now derive the governing equations using Hamilton’s variational principle in a form suitable for system dynamics. Consider the action functional

\[ \mathcal{S}[q] = \int_{t_0}^{t_1} L(q,\dot{q},t)\,dt. \]

For variations \( \delta q(t) \) satisfying \( \delta q(t_0)=\delta q(t_1)=0 \), stationarity of the action plus virtual work of non-conservative forces gives

\[ \delta \mathcal{S} + \int_{t_0}^{t_1}\sum_{i=1}^{n} Q_i \delta q_i \, dt = 0. \]

Expand the variation:

\[ \delta \mathcal{S} = \int_{t_0}^{t_1} \sum_{i=1}^{n} \left( \frac{\partial L}{\partial q_i}\delta q_i + \frac{\partial L}{\partial \dot{q}_i}\delta \dot{q}_i \right)dt. \]

Integrate the second term by parts:

\[ \int_{t_0}^{t_1}\frac{\partial L}{\partial \dot{q}_i}\delta \dot{q}_i\,dt = \left[\frac{\partial L}{\partial \dot{q}_i}\delta q_i\right]_{t_0}^{t_1} - \int_{t_0}^{t_1}\frac{d}{dt}\left(\frac{\partial L}{\partial \dot{q}_i}\right)\delta q_i\,dt. \]

The boundary term is zero because the endpoint variations vanish. Therefore

\[ \delta \mathcal{S} = \int_{t_0}^{t_1}\sum_{i=1}^{n} \left( \frac{\partial L}{\partial q_i} - \frac{d}{dt}\left(\frac{\partial L}{\partial \dot{q}_i}\right) \right)\delta q_i\,dt. \]

Substituting the generalized work term and using the arbitrariness of \( \delta q_i \) yields

\[ \frac{d}{dt}\left(\frac{\partial L}{\partial \dot{q}_i}\right) - \frac{\partial L}{\partial q_i} = Q_i. \]

If Rayleigh dissipation is used, include \( \partial \mathcal{D}/\partial \dot{q}_i \) on the left-hand side, which is equivalent to moving viscous generalized forces into a structured dissipation term.

flowchart TD
  A["Action S = integral L dt"] --> B["Take variation delta S"]
  B --> C["Use delta q(t0)=delta q(t1)=0"]
  C --> D["Integrate by parts"]
  D --> E["Add virtual work term integral Q*delta q dt"]
  E --> F["Collect coefficients of delta q_i"]
  F --> G["Euler-Lagrange equations"]
        

For standard mechanical systems with generalized coordinates \( q \in \mathbb{R}^n \), the Euler–Lagrange equations can be reorganized as

\[ M(q)\ddot{q} + h(q,\dot{q}) + g(q) + R\dot{q} = B(q)u, \]

where \( M(q) \) is the mass/inertia matrix, \( h(q,\dot{q}) \) contains velocity-coupling (centrifugal/Coriolis-type) terms, \( g(q)=\partial V/\partial q \), and \( R\dot{q} \) represents linear viscous dissipation. In constrained form:

\[ M(q)\ddot{q} + h(q,\dot{q}) + g(q) = B(q)u + A(q,t)^T\lambda, \qquad \phi(q,t)=0. \]

Under physically meaningful parameters (positive masses and inertias), the matrix \( M(q) \) is symmetric and positive definite in admissible coordinates, which guarantees local solvability for \( \ddot{q} \) in unconstrained systems.

The total mechanical energy \( E(q,\dot{q})=T(q,\dot{q})+V(q) \) satisfies the power balance

\[ \frac{dE}{dt} = \dot{q}^T Q - \dot{q}^T R\dot{q} + \frac{\partial L}{\partial t}. \]

For autonomous systems with no explicit time dependence and no external actuation \( Q=0 \), energy decreases monotonically if \( R \) is positive semidefinite, and is conserved if \( R=0 \).

We model a cart of mass \( M \) connected to a spring \( k \) and damper \( c \), carrying a pendulum bob of mass \( m \) and center-of-mass distance \( l \). The pendulum may also have rotational inertia \( J \) about the pivot and viscous joint damping \( b \). A horizontal actuator force \( u \) acts on the cart.

Choose generalized coordinates \( q_1=x \) (cart displacement) and \( q_2=\theta \) (pendulum angle from the downward vertical). The pendulum mass center is

\[ x_p = x + l\sin\theta, \qquad y_p = -l\cos\theta. \]

Differentiating gives

\[ \dot{x}_p = \dot{x} + l\cos\theta\,\dot{\theta}, \qquad \dot{y}_p = l\sin\theta\,\dot{\theta}. \]

Therefore the kinetic energy is

\[ \begin{aligned} T &= \frac{1}{2}M\dot{x}^2 + \frac{1}{2}m\left(\dot{x}_p^2+\dot{y}_p^2\right) + \frac{1}{2}J\dot{\theta}^2 \\ &= \frac{1}{2}(M+m)\dot{x}^2 + ml\cos\theta\,\dot{x}\dot{\theta} + \frac{1}{2}(J+ml^2)\dot{\theta}^2. \end{aligned} \]

Potential energy and Rayleigh dissipation are

\[ V = \frac{1}{2}kx^2 + mgl(1-\cos\theta), \qquad \mathcal{D} = \frac{1}{2}c\dot{x}^2 + \frac{1}{2}b\dot{\theta}^2. \]

With generalized non-conservative input \( Q=[u,\;0]^T \), the Euler–Lagrange equations become

\[ (M+m)\ddot{x} + ml\cos\theta\,\ddot{\theta} - ml\sin\theta\,\dot{\theta}^2 + c\dot{x} + kx = u \]

\[ ml\cos\theta\,\ddot{x} + (J+ml^2)\ddot{\theta} + b\dot{\theta} + mgl\sin\theta = 0. \]

In compact matrix form:

\[ \underbrace{\begin{bmatrix} M+m & ml\cos\theta \\ ml\cos\theta & J+ml^2 \end{bmatrix} }_{M(q)} \begin{bmatrix}\ddot{x}\\ \ddot{\theta}\end{bmatrix} + \underbrace{\begin{bmatrix} -ml\sin\theta\,\dot{\theta}^2 \\ 0 \end{bmatrix} }_{h(q,\dot{q})} + \underbrace{\begin{bmatrix} kx \\ mgl\sin\theta \end{bmatrix} }_{g(q)} + \underbrace{\begin{bmatrix} c\dot{x} \\ b\dot{\theta} \end{bmatrix} }_{R\dot{q} } = \begin{bmatrix} u \\ 0 \end{bmatrix}. \]

The determinant of the inertia matrix is

\[ \Delta = (M+m)(J+ml^2) - m^2l^2\cos^2\theta, \]

which remains positive for practical parameter ranges with positive masses and inertia. Hence \( \ddot{x} \) and \( \ddot{\theta} \) can be solved explicitly.

For small angles around the hanging equilibrium \( \theta \approx 0 \), using \( \sin\theta \approx \theta \) and \( \cos\theta \approx 1 \), the linearized equations are

\[ (M+m)\ddot{x} + ml\ddot{\theta} + c\dot{x} + kx = u, \]

\[ ml\ddot{x} + (J+ml^2)\ddot{\theta} + b\dot{\theta} + mgl\,\theta = 0. \]

The Python implementation below verifies the Euler–Lagrange equations symbolically (using sympy) and simulates the nonlinear system (using scipy.integrate.solve_ivp). It also exports state trajectories and energy to a CSV file.

Code: Chapter18_Lesson2.py

# Chapter18_Lesson2.py
# Lagrangian modeling of a cart-pendulum system using generalized coordinates q=[x, theta]
# Includes (1) symbolic derivation check with SymPy, and (2) numerical simulation with SciPy.

import numpy as np
from math import sin, cos
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

try:
    import sympy as sp
    SYMPY_AVAILABLE = True
except Exception:
    SYMPY_AVAILABLE = False


def symbolic_derivation():
    if not SYMPY_AVAILABLE:
        print("SymPy not available; skipping symbolic derivation.")
        return

    t = sp.symbols('t', real=True)
    M, m, l, J, g, k, c, b, u = sp.symbols('M m l J g k c b u', positive=True, real=True)
    x = sp.Function('x')(t)
    th = sp.Function('th')(t)
    xd = sp.diff(x, t)
    thd = sp.diff(th, t)

    T = sp.Rational(1, 2) * (M + m) * xd**2 + m * l * sp.cos(th) * xd * thd + sp.Rational(1, 2) * (J + m * l**2) * thd**2
    V = sp.Rational(1, 2) * k * x**2 + m * g * l * (1 - sp.cos(th))
    D = sp.Rational(1, 2) * c * xd**2 + sp.Rational(1, 2) * b * thd**2
    L = T - V

    q = [x, th]
    qd = [xd, thd]
    Q = [u, 0]

    eom = []
    for qi, qdi, Qi in zip(q, qd, Q):
        expr = sp.diff(sp.diff(L, qdi), t) - sp.diff(L, qi) + sp.diff(D, qdi) - Qi
        eom.append(sp.simplify(sp.expand(expr)))

    print("\nSymbolic Euler-Lagrange equations (left-hand sides = 0):")
    for i, eq in enumerate(eom, start=1):
        print(f"Eq{i}:", sp.factor(eq))


def cart_pendulum_rhs(t, state, p):
    x, th, xd, thd = state
    M = p["M"]
    m = p["m"]
    l = p["l"]
    J = p["J"]
    g = p["g"]
    k = p["k"]
    c = p["c"]
    b = p["b"]

    # Example actuation: smooth pulse on the cart
    u = 2.0 if 0.5 <= t <= 2.5 else 0.0

    # M(q) qdd = rhs
    M11 = M + m
    M12 = m * l * cos(th)
    M21 = M12
    M22 = J + m * l * l

    rhs1 = u - c * xd - k * x + m * l * sin(th) * thd * thd
    rhs2 = -b * thd - m * g * l * sin(th)

    det = M11 * M22 - M12 * M21
    xdd = (rhs1 * M22 - rhs2 * M12) / det
    thdd = (M11 * rhs2 - M21 * rhs1) / det
    return np.array([xd, thd, xdd, thdd], dtype=float)


def total_energy(state, p):
    x, th, xd, thd = state
    M = p["M"]
    m = p["m"]
    l = p["l"]
    J = p["J"]
    g = p["g"]
    k = p["k"]

    T = 0.5 * (M + m) * xd * xd + m * l * np.cos(th) * xd * thd + 0.5 * (J + m * l * l) * thd * thd
    V = 0.5 * k * x * x + m * g * l * (1.0 - np.cos(th))
    return T + V


def simulate():
    p = dict(M=1.0, m=0.25, l=0.5, J=0.02, g=9.81, k=8.0, c=0.35, b=0.05)
    x0 = np.array([0.05, 0.35, 0.0, 0.0], dtype=float)

    sol = solve_ivp(
        fun=lambda t, y: cart_pendulum_rhs(t, y, p),
        t_span=(0.0, 10.0),
        y0=x0,
        method="RK45",
        max_step=0.01,
        rtol=1e-8,
        atol=1e-9,
    )

    E = np.array([total_energy(sol.y[:, i], p) for i in range(sol.y.shape[1])])

    fig1 = plt.figure(figsize=(8, 4))
    plt.plot(sol.t, sol.y[0], label="x (m)")
    plt.plot(sol.t, sol.y[1], label="theta (rad)")
    plt.xlabel("Time (s)")
    plt.ylabel("States")
    plt.title("Cart-Pendulum Generalized Coordinates")
    plt.grid(True)
    plt.legend()
    fig1.tight_layout()
    fig1.savefig("Chapter18_Lesson2_states.png", dpi=150)

    fig2 = plt.figure(figsize=(8, 4))
    plt.plot(sol.t, E)
    plt.xlabel("Time (s)")
    plt.ylabel("Energy (J)")
    plt.title("Total Mechanical Energy (with damping and input pulse)")
    plt.grid(True)
    fig2.tight_layout()
    fig2.savefig("Chapter18_Lesson2_energy.png", dpi=150)

    out = np.column_stack([sol.t, sol.y.T, E])
    np.savetxt(
        "Chapter18_Lesson2_simulation.csv",
        out,
        delimiter=",",
        header="t,x,theta,xdot,thetadot,energy",
        comments=""
    )
    print("Saved: Chapter18_Lesson2_states.png, Chapter18_Lesson2_energy.png, Chapter18_Lesson2_simulation.csv")


if __name__ == "__main__":
    symbolic_derivation()
    simulate()
      

This C++17 code implements the same nonlinear model using a hand-written fourth-order Runge–Kutta integrator. It is useful for embedded-oriented workflows and for understanding solver structure.

Code: Chapter18_Lesson2.cpp

// Chapter18_Lesson2.cpp
// Lagrangian cart-pendulum simulation (generalized coordinates q = [x, theta])
// Standalone C++17 code with RK4 integration and CSV output.

#include <cmath>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <array>
#include <string>

struct Params {
    double M{1.0};
    double m{0.25};
    double l{0.5};
    double J{0.02};
    double g{9.81};
    double k{8.0};
    double c{0.35};
    double b{0.05};
};

using State = std::array<double, 4>; // [x, theta, xdot, thetadot]

double control_input(double t) {
    return (t >= 0.5 && t <= 2.5) ? 2.0 : 0.0;
}

State rhs(double t, const State& s, const Params& p) {
    const double x = s[0];
    const double th = s[1];
    const double xd = s[2];
    const double thd = s[3];

    const double M11 = p.M + p.m;
    const double M12 = p.m * p.l * std::cos(th);
    const double M21 = M12;
    const double M22 = p.J + p.m * p.l * p.l;

    const double u = control_input(t);
    const double rhs1 = u - p.c * xd - p.k * x + p.m * p.l * std::sin(th) * thd * thd;
    const double rhs2 = -p.b * thd - p.m * p.g * p.l * std::sin(th);

    const double det = M11 * M22 - M12 * M21;
    const double xdd = (rhs1 * M22 - rhs2 * M12) / det;
    const double thdd = (M11 * rhs2 - M21 * rhs1) / det;

    return {xd, thd, xdd, thdd};
}

State add_scaled(const State& a, const State& b, double scale) {
    return {a[0] + scale * b[0], a[1] + scale * b[1], a[2] + scale * b[2], a[3] + scale * b[3]};
}

double total_energy(const State& s, const Params& p) {
    const double x = s[0];
    const double th = s[1];
    const double xd = s[2];
    const double thd = s[3];

    const double T = 0.5 * (p.M + p.m) * xd * xd
                   + p.m * p.l * std::cos(th) * xd * thd
                   + 0.5 * (p.J + p.m * p.l * p.l) * thd * thd;
    const double V = 0.5 * p.k * x * x + p.m * p.g * p.l * (1.0 - std::cos(th));
    return T + V;
}

int main() {
    Params p;
    State s{0.05, 0.35, 0.0, 0.0};

    const double t0 = 0.0;
    const double tf = 10.0;
    const double h = 0.001;
    const int nSteps = static_cast<int>((tf - t0) / h);

    std::ofstream file("Chapter18_Lesson2_cpp_output.csv");
    file << "t,x,theta,xdot,thetadot,energy\n";
    file << std::fixed << std::setprecision(8);

    double t = t0;
    for (int k = 0; k <= nSteps; ++k) {
        file << t << "," << s[0] << "," << s[1] << "," << s[2] << "," << s[3] << "," << total_energy(s, p) << "\n";

        State k1 = rhs(t, s, p);
        State k2 = rhs(t + 0.5 * h, add_scaled(s, k1, 0.5 * h), p);
        State k3 = rhs(t + 0.5 * h, add_scaled(s, k2, 0.5 * h), p);
        State k4 = rhs(t + h, add_scaled(s, k3, h), p);

        for (int i = 0; i < 4; ++i) {
            s[i] += (h / 6.0) * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]);
        }
        t += h;
    }

    file.close();

    std::cout << "Simulation finished. Output saved to Chapter18_Lesson2_cpp_output.csv\n";
    std::cout << "Compile example: g++ -O2 -std=c++17 Chapter18_Lesson2.cpp -o Chapter18_Lesson2\n";
    return 0;
}
      

The Java version mirrors the C++ solver and demonstrates how Lagrangian system simulation can be implemented in a platform-neutral environment without external numerical libraries.

Code: Chapter18_Lesson2.java

// Chapter18_Lesson2.java
// Lagrangian cart-pendulum simulation (generalized coordinates q = [x, theta])
// Pure Java implementation with RK4 integration and CSV export.

import java.io.FileWriter;
import java.io.IOException;
import java.io.PrintWriter;

public class Chapter18_Lesson2 {

    static class Params {
        double M = 1.0;
        double m = 0.25;
        double l = 0.5;
        double J = 0.02;
        double g = 9.81;
        double k = 8.0;
        double c = 0.35;
        double b = 0.05;
    }

    static double controlInput(double t) {
        return (t >= 0.5 && t <= 2.5) ? 2.0 : 0.0;
    }

    static double[] rhs(double t, double[] s, Params p) {
        double x = s[0];
        double th = s[1];
        double xd = s[2];
        double thd = s[3];

        double M11 = p.M + p.m;
        double M12 = p.m * p.l * Math.cos(th);
        double M21 = M12;
        double M22 = p.J + p.m * p.l * p.l;

        double u = controlInput(t);
        double rhs1 = u - p.c * xd - p.k * x + p.m * p.l * Math.sin(th) * thd * thd;
        double rhs2 = -p.b * thd - p.m * p.g * p.l * Math.sin(th);

        double det = M11 * M22 - M12 * M21;
        double xdd = (rhs1 * M22 - rhs2 * M12) / det;
        double thdd = (M11 * rhs2 - M21 * rhs1) / det;

        return new double[]{xd, thd, xdd, thdd};
    }

    static double[] addScaled(double[] a, double[] b, double scale) {
        return new double[]{
            a[0] + scale * b[0],
            a[1] + scale * b[1],
            a[2] + scale * b[2],
            a[3] + scale * b[3]
        };
    }

    static double totalEnergy(double[] s, Params p) {
        double x = s[0];
        double th = s[1];
        double xd = s[2];
        double thd = s[3];

        double T = 0.5 * (p.M + p.m) * xd * xd
                 + p.m * p.l * Math.cos(th) * xd * thd
                 + 0.5 * (p.J + p.m * p.l * p.l) * thd * thd;
        double V = 0.5 * p.k * x * x + p.m * p.g * p.l * (1.0 - Math.cos(th));
        return T + V;
    }

    public static void main(String[] args) throws IOException {
        Params p = new Params();
        double[] s = new double[]{0.05, 0.35, 0.0, 0.0};

        double t0 = 0.0;
        double tf = 10.0;
        double h = 0.001;
        int nSteps = (int)((tf - t0) / h);

        try (PrintWriter out = new PrintWriter(new FileWriter("Chapter18_Lesson2_java_output.csv"))) {
            out.println("t,x,theta,xdot,thetadot,energy");

            double t = t0;
            for (int k = 0; k <= nSteps; k++) {
                out.printf(java.util.Locale.US, "%.8f,%.8f,%.8f,%.8f,%.8f,%.8f%n",
                        t, s[0], s[1], s[2], s[3], totalEnergy(s, p));

                double[] k1 = rhs(t, s, p);
                double[] k2 = rhs(t + 0.5 * h, addScaled(s, k1, 0.5 * h), p);
                double[] k3 = rhs(t + 0.5 * h, addScaled(s, k2, 0.5 * h), p);
                double[] k4 = rhs(t + h, addScaled(s, k3, h), p);

                for (int i = 0; i < 4; i++) {
                    s[i] += (h / 6.0) * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]);
                }
                t += h;
            }
        }

        System.out.println("Simulation finished. Output saved to Chapter18_Lesson2_java_output.csv");
        System.out.println("Compile: javac Chapter18_Lesson2.java");
        System.out.println("Run: java Chapter18_Lesson2");
    }
}
      

The MATLAB code uses ode45 and exports results to CSV. For Simulink, the same dynamics can be implemented using a MATLAB Function block that computes \( \dot{x} \) from the state and input, connected to an Integrator chain for the state vector.

Code: Chapter18_Lesson2.m

% Chapter18_Lesson2.m
% Lagrangian modeling and simulation of a cart-pendulum system
% Generalized coordinates: q = [x; theta]
% Uses ode45 and exports simulation results.

clear; clc;

p.M = 1.0;
p.m = 0.25;
p.l = 0.5;
p.J = 0.02;
p.g = 9.81;
p.k = 8.0;
p.c = 0.35;
p.b = 0.05;

x0 = [0.05; 0.35; 0.0; 0.0]; % [x; theta; xdot; thetadot]

tspan = [0 10];
opts = odeset('RelTol',1e-8,'AbsTol',1e-9);
[t, X] = ode45(@(t,x) cartPendulumRHS(t, x, p), tspan, x0, opts);

E = zeros(length(t),1);
for i = 1:length(t)
    E(i) = totalEnergy(X(i,:).', p);
end

figure;
plot(t, X(:,1), 'LineWidth', 1.2); hold on;
plot(t, X(:,2), 'LineWidth', 1.2);
grid on; xlabel('Time (s)'); ylabel('States');
legend('x (m)','theta (rad)');
title('Cart-Pendulum Generalized Coordinates');

figure;
plot(t, E, 'LineWidth', 1.2);
grid on; xlabel('Time (s)'); ylabel('Energy (J)');
title('Total Mechanical Energy');

T = table(t, X(:,1), X(:,2), X(:,3), X(:,4), E, ...
    'VariableNames', {'t','x','theta','xdot','thetadot','energy'});
writetable(T, 'Chapter18_Lesson2_matlab_output.csv');

disp('Saved Chapter18_Lesson2_matlab_output.csv');

% Optional symbolic verification (requires Symbolic Math Toolbox)
% syms x(t) th(t) M m l J g k c b u real
% xd = diff(x,t); thd = diff(th,t);
% T = 1/2*(M+m)*xd^2 + m*l*cos(th)*xd*thd + 1/2*(J+m*l^2)*thd^2;
% V = 1/2*k*x^2 + m*g*l*(1-cos(th));
% D = 1/2*c*xd^2 + 1/2*b*thd^2;
% L = T - V;
% Eqx  = simplify(diff(diff(L, xd), t) - diff(L, x) + diff(D, xd) - u);
% Eqth = simplify(diff(diff(L, thd), t) - diff(L, th) + diff(D, thd));
% pretty(Eqx), pretty(Eqth)

function dx = cartPendulumRHS(t, x, p)
    q1 = x(1); th = x(2); q1d = x(3); thd = x(4);

    if t >= 0.5 && t <= 2.5
        u = 2.0;
    else
        u = 0.0;
    end

    M11 = p.M + p.m;
    M12 = p.m * p.l * cos(th);
    M21 = M12;
    M22 = p.J + p.m * p.l^2;

    rhs1 = u - p.c*q1d - p.k*q1 + p.m*p.l*sin(th)*thd^2;
    rhs2 = -p.b*thd - p.m*p.g*p.l*sin(th);

    detM = M11*M22 - M12*M21;
    q1dd = (rhs1*M22 - rhs2*M12) / detM;
    thdd = (M11*rhs2 - M21*rhs1) / detM;

    dx = [q1d; thd; q1dd; thdd];
end

function E = totalEnergy(x, p)
    q1 = x(1); th = x(2); q1d = x(3); thd = x(4);
    T = 0.5*(p.M + p.m)*q1d^2 + p.m*p.l*cos(th)*q1d*thd + 0.5*(p.J + p.m*p.l^2)*thd^2;
    V = 0.5*p.k*q1^2 + p.m*p.g*p.l*(1 - cos(th));
    E = T + V;
end
      

Simulink mapping (conceptual): Use state vector \( [x,\theta,\dot{x},\dot{\theta}]^T \), compute \( \ddot{x}, \ddot{\theta} \) inside a MATLAB Function block from the inertia matrix and right-hand side terms, then integrate to produce \( \dot{x}, \dot{\theta}, x, \theta \). Add scopes for state and energy monitoring.

The Mathematica code derives the Euler–Lagrange expressions symbolically and solves the nonlinear model with NDSolve. It also exports plots and time-series data.

Code: Chapter18_Lesson2.nb

(* Chapter18_Lesson2.nb *)
(* Wolfram Language notebook-style script for Lagrangian modeling of a cart-pendulum *)
(* Generalized coordinates: q = {x[t], th[t]} *)

ClearAll["Global`*"];

params = <|
  "M" -> 1.0, "m" -> 0.25, "l" -> 0.5, "J" -> 0.02,
  "g" -> 9.81, "k" -> 8.0, "c" -> 0.35, "b" -> 0.05
|>;

x = x[t];
th = th[t];

Texpr = 1/2 (M + m) D[x, t]^2 + m l Cos[th] D[x, t] D[th, t] + 1/2 (J + m l^2) D[th, t]^2;
Vexpr = 1/2 k x^2 + m g l (1 - Cos[th]);
Dexpr = 1/2 c D[x, t]^2 + 1/2 b D[th, t]^2;
Lexpr = Texpr - Vexpr;

qList = {x, th};
qdList = D[qList, t];
Qnc = {u[t], 0};

eulerLagrange = Table[
  FullSimplify[D[D[Lexpr, qdList[[i]]], t] - D[Lexpr, qList[[i]]] + D[Dexpr, qdList[[i]]] - Qnc[[i]]],
  {i, 1, 2}
];

Print["Euler-Lagrange equations (left-hand side = 0):"];
Print /@ eulerLagrange;

u[t_] := Piecewise[{ {2.0, 0.5 <= t <= 2.5} }, 0.0];

eqns = {
  (M + m) x''[t] + m l Cos[th[t]] th''[t] - m l Sin[th[t]] th'[t]^2 + c x'[t] + k x[t] == u[t],
  m l Cos[th[t]] x''[t] + (J + m l^2) th''[t] + b th'[t] + m g l Sin[th[t]] == 0
} /. {M -> params["M"], m -> params["m"], l -> params["l"], J -> params["J"],
       g -> params["g"], k -> params["k"], c -> params["c"], b -> params["b"]};

ics = {
  x[0] == 0.05, th[0] == 0.35,
  x'[0] == 0.0, th'[0] == 0.0
};

sol = NDSolve[Join[eqns, ics], {x, th}, {t, 0, 10}, MaxStepFraction -> 1/100][[1]];

energy[t_] := Evaluate[
  (Texpr + Vexpr) /. {M -> params["M"], m -> params["m"], l -> params["l"], J -> params["J"],
                      g -> params["g"], k -> params["k"]} /. sol
];

p1 = Plot[Evaluate[{x[t], th[t]} /. sol], {t, 0, 10},
  PlotLegends -> {"x(t)", "th(t)"},
  GridLines -> Automatic,
  PlotLabel -> "Generalized Coordinates"];

p2 = Plot[energy[t], {t, 0, 10},
  GridLines -> Automatic,
  PlotLabel -> "Total Mechanical Energy"];

Export["Chapter18_Lesson2_mathematica_states.png", p1];
Export["Chapter18_Lesson2_mathematica_energy.png", p2];

data = Table[
  {tt, x[tt] /. sol, th[tt] /. sol, (x'[tt] /. sol), (th'[tt] /. sol), energy[tt]},
  {tt, 0, 10, 0.01}
];

Export["Chapter18_Lesson2_mathematica_output.csv",
  Prepend[data, {"t","x","theta","xdot","thetadot","energy"}]
];

Print["Saved Mathematica plots and CSV output."];
      

Problem 1 (Single Pendulum via Lagrange): A simple pendulum of mass \( m \) and length \( l \) moves in a vertical plane with generalized coordinate \( \theta \) measured from the downward vertical. Derive the equation of motion assuming no damping and no actuation.

Solution: The kinetic and potential energies are

\[ T=\frac{1}{2}ml^2\dot{\theta}^2, \qquad V=mgl(1-\cos\theta). \]

Thus \( L=T-V \). Applying Euler–Lagrange:

\[ \frac{d}{dt}\left(\frac{\partial L}{\partial \dot{\theta} }\right)-\frac{\partial L}{\partial \theta}=0 \quad \Rightarrow \quad ml^2\ddot{\theta}+mgl\sin\theta=0. \]

Dividing by \( ml^2 \) gives \( \ddot{\theta} + (g/l)\sin\theta = 0 \).

Problem 2 (Generalized Force from Cartesian Force): A particle moves on a line with generalized coordinate \( q \), and the Cartesian position is \( x = a q^2 \). A force \( F_x \) acts along the x-axis. Find the generalized force \( Q_q \).

Solution: By virtual work, \( \delta W = F_x \delta x = Q_q \delta q \) and \( \delta x = \frac{\partial x}{\partial q}\delta q = 2aq\,\delta q \). Therefore

\[ Q_q = F_x \frac{\partial x}{\partial q} = 2aqF_x. \]

Problem 3 (Mass Matrix Positivity): For the cart–pendulum system in this lesson, show that the inertia matrix is symmetric. Under what condition is it invertible?

Solution: The inertia matrix is

\[ M(q)=\begin{bmatrix} M+m & ml\cos\theta \\ ml\cos\theta & J+ml^2 \end{bmatrix}, \]

which is symmetric because the off-diagonal entries are equal. It is invertible when \( \det(M(q))=\Delta \neq 0 \). A sufficient practical condition is positive physical parameters and non-degenerate geometry, which yields \( \Delta > 0 \) over the operating range.

Problem 4 (Linearization about the Hanging Equilibrium): Linearize the cart–pendulum equations around \( x=0 \), \( \theta=0 \), \( \dot{x}=0 \), \( \dot{\theta}=0 \).

Solution: Use \( \sin\theta \approx \theta \), \( \cos\theta \approx 1 \), and neglect the second-order term \( \theta\dot{\theta}^2 \). The nonlinear equations reduce to

\[ (M+m)\ddot{x} + ml\ddot{\theta} + c\dot{x} + kx = u, \]

\[ ml\ddot{x} + (J+ml^2)\ddot{\theta} + b\dot{\theta} + mgl\,\theta = 0. \]

This is a linear coupled second-order model suitable for state-space conversion in later analysis.

Problem 5 (Energy Dissipation Check): For a mechanical system with \( L=T-V \) and Rayleigh dissipation \( \mathcal{D}=\frac{1}{2}\dot{q}^TR\dot{q} \), show that in the absence of external input the rate of total energy satisfies \( \frac{dE}{dt}=-\dot{q}^TR\dot{q} \) for autonomous systems.

Solution: Multiply each Euler–Lagrange equation by \( \dot{q}_i \) and sum over \( i \):

\[ \sum_i \dot{q}_i \left( \frac{d}{dt}\frac{\partial L}{\partial \dot{q}_i} - \frac{\partial L}{\partial q_i} + \frac{\partial \mathcal{D} }{\partial \dot{q}_i} \right)=0. \]

Rearranging and using standard identities for \( E=\sum_i \dot{q}_i\frac{\partial L}{\partial \dot{q}_i}-L \) gives

\[ \frac{dE}{dt} = -\dot{q}^T\frac{\partial \mathcal{D} }{\partial \dot{q} } = -\dot{q}^TR\dot{q}. \]

Since \( R \) is positive semidefinite, the energy is non-increasing.

We introduced generalized coordinates, virtual work, and the Euler–Lagrange formulation as a rigorous energy-based method for deriving system dynamics. The lesson derived the equations from a variational argument, showed the standard matrix form used in control and simulation, and developed a complete nonlinear cart–pendulum example with damping and actuation. The accompanying Python, C++, Java, MATLAB, and Mathematica implementations demonstrate how the same Lagrangian model can be simulated and validated across platforms.