In earlier chapters, we modeled systems using Newton/Kirchhoff laws and state-space equations. In Lesson 2 of this chapter, we introduced the Lagrangian \( L(q,\dot{q},t)=T-V \). The Hamiltonian perspective re-parameterizes the dynamics using generalized coordinates \( q \) and conjugate momenta \( p \), and introduces a scalar energy function \( H(q,p,t) \).

For many physical systems, especially conservative ones, \( H \) equals total stored energy. For open systems (with inputs and dissipation), the port-Hamiltonian (pH) form preserves the same energy-centered interpretation while explicitly exposing exchanged power through input/output ports.

\[ \text{Hamiltonian (conservative):}\quad \dot{x} = \mathbf{J}\nabla H(x),\qquad \mathbf{J}^\top = -\mathbf{J} \]

\[ \text{Port-Hamiltonian (open system):}\quad \dot{x} = \big(\mathbf{J}(x)-\mathbf{R}(x)\big)\nabla H(x)+\mathbf{g}(x)u,\qquad y=\mathbf{g}(x)^\top \nabla H(x) \]

Here \( x \) is the energy state, \( \mathbf{J}(x) \) encodes interconnection (power-conserving exchange), \( \mathbf{R}(x) \) encodes dissipation, and \( (u,y) \) are power-port variables.

Consider an \( n \)-DOF system with generalized coordinates \( q=[q_1,\dots,q_n]^\top \). If the Lagrangian is regular, meaning

\[ \det\!\left(\frac{\partial^2 L}{\partial \dot{q}^2}\right)\neq 0, \]

then the momentum mapping is locally invertible: \( p_i=\frac{\partial L}{\partial \dot{q}_i} \), and we can solve for \( \dot{q}=\dot{q}(q,p,t) \).

The Hamiltonian is defined as the Legendre transform of \( L \) with respect to \( \dot{q} \):

\[ H(q,p,t)=\sum_{i=1}^{n} p_i \dot{q}_i - L(q,\dot{q},t), \quad \text{with } \dot{q}=\dot{q}(q,p,t). \]

Derivation of Hamilton's equations. Take the differential of \( H \):

\[ dH = \sum_{i=1}^{n}\dot{q}_i\,dp_i + \sum_{i=1}^{n}p_i\,d\dot{q}_i - \sum_{i=1}^{n}\frac{\partial L}{\partial q_i}dq_i - \sum_{i=1}^{n}\frac{\partial L}{\partial \dot{q}_i}d\dot{q}_i - \frac{\partial L}{\partial t}dt. \]

Using \( p_i=\frac{\partial L}{\partial \dot{q}_i} \), the \( d\dot{q}_i \) terms cancel:

\[ dH = \sum_{i=1}^{n}\dot{q}_i\,dp_i - \sum_{i=1}^{n}\frac{\partial L}{\partial q_i}dq_i - \frac{\partial L}{\partial t}dt. \]

By Euler–Lagrange equations from the previous lesson, \( \frac{d}{dt}\frac{\partial L}{\partial \dot{q}_i}-\frac{\partial L}{\partial q_i}=0 \), hence \( \dot{p}_i=\frac{\partial L}{\partial q_i} \). Comparing with the standard differential form

\[ dH = \sum_{i=1}^{n}\frac{\partial H}{\partial q_i}dq_i +\sum_{i=1}^{n}\frac{\partial H}{\partial p_i}dp_i +\frac{\partial H}{\partial t}dt, \]

we obtain the canonical equations:

\[ \boxed{\dot{q}_i=\frac{\partial H}{\partial p_i}},\qquad \boxed{\dot{p}_i=-\frac{\partial H}{\partial q_i}},\qquad \boxed{\frac{\partial H}{\partial t}=-\frac{\partial L}{\partial t}}. \]

flowchart TD
  A["Physical system and generalized coordinates q"] --> B["Lagrangian L(q,qdot,t) = T - V"]
  B --> C["Momentum definition p = dL/dqdot"]
  C --> D["Invert map qdot = qdot(q,p,t)"]
  D --> E["Hamiltonian H = p.qdot - L"]
  E --> F["Canonical equations for q and p"]
  F --> G["Energy interpretation and conservation test"]
        

Define the state vector \( x=\begin{bmatrix}q \\ p\end{bmatrix}\in\mathbb{R}^{2n} \) and the canonical skew-symmetric matrix

\[ \mathbf{J}= \begin{bmatrix} \mathbf{0} & \mathbf{I} \\ -\mathbf{I} & \mathbf{0} \end{bmatrix}, \qquad \mathbf{J}^\top = -\mathbf{J}. \]

Then Hamilton’s equations can be written compactly as

\[ \dot{x} = \mathbf{J}\nabla H(x,t). \]

For any smooth scalar function \( f(q,p,t) \), the time derivative along trajectories is

\[ \frac{df}{dt} = \nabla f^\top \dot{x} + \frac{\partial f}{\partial t} = \nabla f^\top \mathbf{J}\nabla H + \frac{\partial f}{\partial t}. \]

This motivates the (canonical) Poisson bracket:

\[ \{f,H\} := \nabla f^\top \mathbf{J}\nabla H, \qquad \Rightarrow \qquad \frac{df}{dt} = \{f,H\} + \frac{\partial f}{\partial t}. \]

Energy conservation (autonomous conservative systems). Let \( f=H \) and assume \( \frac{\partial H}{\partial t}=0 \). Then

\[ \dot{H} = \nabla H^\top \mathbf{J}\nabla H. \]

Since \( \mathbf{J} \) is skew-symmetric, the quadratic form is identically zero:

\[ \big(\nabla H^\top \mathbf{J}\nabla H\big)^\top = \nabla H^\top \mathbf{J}^\top \nabla H = -\nabla H^\top \mathbf{J}\nabla H \;\Rightarrow\; \nabla H^\top \mathbf{J}\nabla H=0. \]

Therefore,

\[ \boxed{\dot{H}=0} \]

for autonomous conservative Hamiltonian systems. This is the key reason the Hamiltonian formulation is ideal for energy-based analysis.

Real engineering systems are often open: they interact with sources, loads, actuators, and sensors. The port-Hamiltonian form extends the Hamiltonian structure by adding dissipation and external ports:

\[ \dot{x} = \big(\mathbf{J}(x)-\mathbf{R}(x)\big)\nabla H(x) + \mathbf{g}(x)u, \qquad y = \mathbf{g}(x)^\top \nabla H(x), \]

with structural conditions \( \mathbf{J}(x)^\top = -\mathbf{J}(x) \) and \( \mathbf{R}(x)=\mathbf{R}(x)^\top \succeq 0 \).

Passivity / energy-balance proof. Differentiate the Hamiltonian:

\[ \dot{H} = \nabla H^\top \dot{x} = \nabla H^\top(\mathbf{J}-\mathbf{R})\nabla H + \nabla H^\top \mathbf{g}u. \]

The skew term vanishes: \( \nabla H^\top \mathbf{J}\nabla H = 0 \). Thus,

\[ \dot{H} = -\nabla H^\top \mathbf{R}\nabla H + \underbrace{\nabla H^\top \mathbf{g}}_{y^\top}u = -\nabla H^\top \mathbf{R}\nabla H + y^\top u. \]

Since \( \mathbf{R}\succeq 0 \), the dissipation term is nonnegative, so

\[ \boxed{\dot{H}\le y^\top u.} \]

This is the standard passivity inequality: the rate of increase of stored energy is bounded by supplied power.

flowchart TD
  U["Input port u (actuator/source)"] --> SYS["Interconnection + storage + dissipation"]
  SYS --> Y["Output port y (power conjugate)"]
  SYS --> H["Stored energy H(x)"]
  H --> B["Energy rate = supply - dissipation"]
        

One major advantage of the port-Hamiltonian view is that it uses the same mathematical skeleton across domains. The main changes are the physical interpretation of the state variables and the Hamiltonian.

5.1 Mechanical Mass–Spring–Damper with Force Input

Let \( q \) be displacement and \( p=m\dot{q} \) momentum. Define

\[ H(q,p)=\frac{1}{2}kq^2+\frac{p^2}{2m}. \]

Choose

\[ x=\begin{bmatrix}q\\p\end{bmatrix},\quad \mathbf{J}=\begin{bmatrix}0 & 1\\ -1 & 0\end{bmatrix},\quad \mathbf{R}=\begin{bmatrix}0 & 0\\ 0 & d\end{bmatrix},\quad \mathbf{g}=\begin{bmatrix}0\\1\end{bmatrix}, \]

where \( d \) is the viscous damping coefficient and \( u \) is applied force. Since \( \nabla H = \begin{bmatrix}kq \\ p/m\end{bmatrix} \), the pH equations become

\[ \dot{q}=\frac{p}{m}, \qquad \dot{p}=-kq-d\frac{p}{m}+u, \qquad y=\frac{p}{m}. \]

The power port is force–velocity: \( y^\top u = \dot{q}\,u \).

5.2 Series RLC Circuit with Voltage Input

Let the energy variables be inductor flux \( \phi \) and capacitor charge \( q_c \). Then

\[ H(\phi,q_c)=\frac{\phi^2}{2L}+\frac{q_c^2}{2C}. \]

Use the state \( x=\begin{bmatrix}\phi \\ q_c\end{bmatrix} \), with \( \nabla H = \begin{bmatrix}i \\ v_C\end{bmatrix} = \begin{bmatrix}\phi/L \\ q_c/C\end{bmatrix} \). A series resistor \( R_e \) and source voltage \( u \) give

\[ \mathbf{J}=\begin{bmatrix}0 & -1\\ 1 & 0\end{bmatrix},\quad \mathbf{R}=\begin{bmatrix}R_e & 0\\ 0 & 0\end{bmatrix},\quad \mathbf{g}=\begin{bmatrix}1\\0\end{bmatrix}. \]

Hence

\[ \dot{x}=(\mathbf{J}-\mathbf{R})\nabla H + \mathbf{g}u \]

expands to

\[ \dot{\phi}=u-R_e i-v_C,\qquad \dot{q}_c=i, \qquad y=i. \]

The power port is voltage–current: \( y^\top u = iu \), exactly the electrical source power.

5.3 Why this matters for later lessons

The same pH structure appears in mechanics, circuits, fluid systems, and electromechanical couplings. In later lessons (bond graphs and multi-domain energy-based models), the interconnection matrices and ports will become systematic tools for building large models by composition rather than re-deriving every ODE from scratch.

Standard explicit integrators (Euler, RK4) are accurate for many engineering tasks, but in conservative Hamiltonian systems they may exhibit small long-term energy drift. This is a numerical artifact, not physical dissipation.

For the canonical system \( \dot{q}=\frac{\partial H}{\partial p},\; \dot{p}=-\frac{\partial H}{\partial q} \), a simple symplectic Euler update is

\[ p_{k+1}=p_k-h\,\frac{\partial H}{\partial q}(q_k,p_{k+1}), \qquad q_{k+1}=q_k+h\,\frac{\partial H}{\partial p}(q_k,p_{k+1}). \]

In this introductory lesson, we use RK4 in the code examples because it is familiar from Chapter 15 and easy to implement in all languages. However, the scripts also verify the pH power balance numerically:

\[ H(T)-H(0)\approx \int_0^T \left(-\nabla H^\top \mathbf{R}\nabla H + y^\top u\right)\,dt. \]

All implementations below simulate the same mass–spring–damper port-Hamiltonian system, compute the Hamiltonian \( H(t) \), and verify the energy balance numerically. This reinforces the core modeling concept independently of programming language.

Chapter18_Lesson3.py

# Chapter18_Lesson3.py
# Hamiltonian and Port-Hamiltonian simulation (mass-spring-damper with input port)
# x = [q, p], H = 0.5*k*q^2 + p^2/(2*m)

import numpy as np
import matplotlib.pyplot as plt

m = 1.5      # mass
k = 12.0     # spring constant
d = 0.8      # viscous damping coefficient
omega = 1.4  # input frequency

J = np.array([[0.0, 1.0],
              [-1.0, 0.0]])
R = np.array([[0.0, 0.0],
              [0.0, d]])
g = np.array([[0.0],
              [1.0]])

def grad_H(x):
    q, p = x
    return np.array([k * q, p / m])

def H(x):
    q, p = x
    return 0.5 * k * q * q + 0.5 * p * p / m

def u(t):
    return np.sin(omega * t)

def y_output(x):
    # y = g^T * gradH = p/m
    return x[1] / m

def dynamics(t, x):
    gh = grad_H(x)
    dx = (J - R) @ gh + (g[:, 0] * u(t))
    return dx

def rk4_step(t, x, h):
    k1 = dynamics(t, x)
    k2 = dynamics(t + 0.5 * h, x + 0.5 * h * k1)
    k3 = dynamics(t + 0.5 * h, x + 0.5 * h * k2)
    k4 = dynamics(t + h, x + h * k3)
    return x + (h / 6.0) * (k1 + 2*k2 + 2*k3 + k4)

T = 20.0
h = 0.002
N = int(T / h) + 1
t_grid = np.linspace(0.0, T, N)

x = np.array([0.15, 0.0])  # initial [q, p]
X = np.zeros((N, 2))
Hvals = np.zeros(N)
yvals = np.zeros(N)
uvals = np.zeros(N)
diss = np.zeros(N)  # gradH^T R gradH

for i, t in enumerate(t_grid):
    X[i, :] = x
    gh = grad_H(x)
    Hvals[i] = H(x)
    yvals[i] = y_output(x)
    uvals[i] = u(t)
    diss[i] = gh @ (R @ gh)
    if i < N - 1:
        x = rk4_step(t, x, h)

# Numerical check of pH power balance:
# dH/dt = -diss + y*u
supply = yvals * uvals
rhs = -diss + supply
energy_change = Hvals[-1] - Hvals[0]
rhs_integral = np.trapz(rhs, t_grid)
residual = energy_change - rhs_integral

print("Final state [q, p] =", X[-1])
print("Initial energy =", Hvals[0])
print("Final energy   =", Hvals[-1])
print("Energy balance residual (should be small) =", residual)

# Save a CSV for cross-language comparison
data = np.column_stack([t_grid, X[:, 0], X[:, 1], Hvals, uvals, yvals, diss, supply])
header = "t,q,p,H,u,y,dissipation,supply"
np.savetxt("Chapter18_Lesson3_python_output.csv", data, delimiter=",", header=header, comments="")

# Plots
plt.figure(figsize=(10, 4))
plt.plot(t_grid, X[:, 0], label="q(t)")
plt.plot(t_grid, X[:, 1], label="p(t)")
plt.xlabel("Time [s]")
plt.ylabel("States")
plt.title("Port-Hamiltonian states")
plt.legend()
plt.tight_layout()

plt.figure(figsize=(10, 4))
plt.plot(t_grid, Hvals, label="H(t)")
plt.plot(t_grid, np.cumsum(rhs) * h + Hvals[0], label="H(0)+int(-diss+y*u)dt", linestyle="--")
plt.xlabel("Time [s]")
plt.ylabel("Energy")
plt.title("Energy balance verification")
plt.legend()
plt.tight_layout()

plt.show()
      

Chapter18_Lesson3.cpp

// Chapter18_Lesson3.cpp
// Port-Hamiltonian mass-spring-damper simulation using RK4
// Compile: g++ -O2 -std=c++17 Chapter18_Lesson3.cpp -o Chapter18_Lesson3

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

struct State {
    double q;
    double p;
};

const double m = 1.5;
const double k = 12.0;
const double d = 0.8;
const double omega_u = 1.4;

double input_u(double t) {
    return std::sin(omega_u * t);
}

double H(const State& x) {
    return 0.5 * k * x.q * x.q + 0.5 * x.p * x.p / m;
}

double y_output(const State& x) {
    return x.p / m;
}

double dissipation(const State& x) {
    // gradH^T R gradH = d * (p/m)^2
    double v = x.p / m;
    return d * v * v;
}

State dynamics(double t, const State& x) {
    State dx;
    // qdot = p/m
    dx.q = x.p / m;
    // pdot = -k q - d*(p/m) + u
    dx.p = -k * x.q - d * (x.p / m) + input_u(t);
    return dx;
}

State add(const State& a, const State& b_) {
    return {a.q + b_.q, a.p + b_.p};
}

State scale(double s, const State& x) {
    return {s * x.q, s * x.p};
}

State rk4_step(double t, const State& x, double h) {
    State k1 = dynamics(t, x);
    State k2 = dynamics(t + 0.5 * h, add(x, scale(0.5 * h, k1)));
    State k3 = dynamics(t + 0.5 * h, add(x, scale(0.5 * h, k2)));
    State k4 = dynamics(t + h, add(x, scale(h, k3)));
    State sum = add(add(k1, scale(2.0, k2)), add(scale(2.0, k3), k4));
    return add(x, scale(h / 6.0, sum));
}

int main() {
    const double T = 20.0;
    const double h = 0.002;
    const int N = static_cast<int>(T / h) + 1;

    State x{0.15, 0.0};
    double H0 = H(x);
    double rhs_int = 0.0;

    std::ofstream csv("Chapter18_Lesson3_cpp_output.csv");
    csv << "t,q,p,H,u,y,dissipation,supply\\n";
    csv << std::setprecision(12);

    double prev_rhs = 0.0;
    bool first = true;

    for (int i = 0; i < N; ++i) {
        double t = i * h;
        double u = input_u(t);
        double y = y_output(x);
        double diss = dissipation(x);
        double supply = y * u;
        double rhs = -diss + supply;

        if (first) {
            first = false;
        } else {
            rhs_int += 0.5 * h * (prev_rhs + rhs);
        }
        prev_rhs = rhs;

        csv << t << "," << x.q << "," << x.p << "," << H(x) << ","
            << u << "," << y << "," << diss << "," << supply << "\\n";

        if (i < N - 1) {
            x = rk4_step(t, x, h);
        }
    }

    csv.close();

    double Hf = H(x);
    double residual = (Hf - H0) - rhs_int;

    std::cout << "Final state [q, p] = [" << x.q << ", " << x.p << "]\\n";
    std::cout << "Initial energy     = " << H0 << "\\n";
    std::cout << "Final energy       = " << Hf << "\\n";
    std::cout << "Energy balance residual = " << residual << "\\n";
    std::cout << "CSV written to Chapter18_Lesson3_cpp_output.csv\\n";
    return 0;
}
      

Chapter18_Lesson3.java

// Chapter18_Lesson3.java
// Port-Hamiltonian mass-spring-damper simulation using RK4
// Compile: javac Chapter18_Lesson3.java
// Run:     java Chapter18_Lesson3

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

public class Chapter18_Lesson3 {
    static class State {
        double q, p;
        State(double q, double p) { this.q = q; this.p = p; }
    }

    static final double m = 1.5;
    static final double k = 12.0;
    static final double d = 0.8;
    static final double omegaU = 1.4;

    static double u(double t) {
        return Math.sin(omegaU * t);
    }

    static double H(State x) {
        return 0.5 * k * x.q * x.q + 0.5 * x.p * x.p / m;
    }

    static double y(State x) {
        return x.p / m;
    }

    static double dissipation(State x) {
        double v = x.p / m;
        return d * v * v;
    }

    static State dynamics(double t, State x) {
        double qdot = x.p / m;
        double pdot = -k * x.q - d * (x.p / m) + u(t);
        return new State(qdot, pdot);
    }

    static State add(State a, State b) {
        return new State(a.q + b.q, a.p + b.p);
    }

    static State scale(double s, State x) {
        return new State(s * x.q, s * x.p);
    }

    static State rk4Step(double t, State x, double h) {
        State k1 = dynamics(t, x);
        State k2 = dynamics(t + 0.5 * h, add(x, scale(0.5 * h, k1)));
        State k3 = dynamics(t + 0.5 * h, add(x, scale(0.5 * h, k2)));
        State k4 = dynamics(t + h, add(x, scale(h, k3)));
        State sum = add(add(k1, scale(2.0, k2)), add(scale(2.0, k3), k4));
        return add(x, scale(h / 6.0, sum));
    }

    public static void main(String[] args) throws IOException {
        double T = 20.0;
        double h = 0.002;
        int N = (int)(T / h) + 1;

        State x = new State(0.15, 0.0);
        double H0 = H(x);

        double rhsIntegral = 0.0;
        double prevRhs = 0.0;
        boolean first = true;

        try (PrintWriter out = new PrintWriter(new FileWriter("Chapter18_Lesson3_java_output.csv"))) {
            out.println("t,q,p,H,u,y,dissipation,supply");

            for (int i = 0; i < N; i++) {
                double t = i * h;
                double ui = u(t);
                double yi = y(x);
                double diss = dissipation(x);
                double supply = yi * ui;
                double rhs = -diss + supply;

                if (first) {
                    first = false;
                } else {
                    rhsIntegral += 0.5 * h * (prevRhs + rhs);
                }
                prevRhs = rhs;

                out.printf("%.10f,%.10f,%.10f,%.10f,%.10f,%.10f,%.10f,%.10f%n",
                        t, x.q, x.p, H(x), ui, yi, diss, supply);

                if (i < N - 1) {
                    x = rk4Step(t, x, h);
                }
            }
        }

        double Hf = H(x);
        double residual = (Hf - H0) - rhsIntegral;

        System.out.printf("Final state [q, p] = [%.8f, %.8f]%n", x.q, x.p);
        System.out.printf("Initial energy     = %.10f%n", H0);
        System.out.printf("Final energy       = %.10f%n", Hf);
        System.out.printf("Energy balance residual = %.6e%n", residual);
        System.out.println("CSV written to Chapter18_Lesson3_java_output.csv");
    }
}
      

Chapter18_Lesson3.m

This MATLAB script is directly runnable and also maps cleanly to a Simulink implementation: two Integrator blocks for \( q \) and \( p \), a MATLAB Function block for the pH right-hand side, and Scope blocks for states and energy.

% Chapter18_Lesson3.m
% Port-Hamiltonian mass-spring-damper simulation (RK4) + energy-balance check
% Can be run in MATLAB or GNU Octave.

clear; clc;

m = 1.5;
k = 12.0;
d = 0.8;
omega = 1.4;

u = @(t) sin(omega*t);
H = @(x) 0.5*k*x(1)^2 + 0.5*(x(2)^2)/m;
yout = @(x) x(2)/m;
diss = @(x) d*(x(2)/m)^2;
f = @(t,x) [x(2)/m;
           -k*x(1) - d*(x(2)/m) + u(t)];

T = 20.0;
h = 0.002;
t = 0:h:T;
N = numel(t);

x = zeros(2,N);
x(:,1) = [0.15; 0.0];

Hvals = zeros(1,N);
yvals = zeros(1,N);
uvals = zeros(1,N);
dvals = zeros(1,N);

for i = 1:N
    xi = x(:,i);
    ti = t(i);
    Hvals(i) = H(xi);
    yvals(i) = yout(xi);
    uvals(i) = u(ti);
    dvals(i) = diss(xi);

    if i < N
        k1 = f(ti, xi);
        k2 = f(ti + 0.5*h, xi + 0.5*h*k1);
        k3 = f(ti + 0.5*h, xi + 0.5*h*k2);
        k4 = f(ti + h, xi + h*k3);
        x(:,i+1) = xi + (h/6)*(k1 + 2*k2 + 2*k3 + k4);
    end
end

supplyVals = yvals .* uvals;
rhs = -dvals + supplyVals;
rhsIntegral = trapz(t, rhs);
residual = (Hvals(end) - Hvals(1)) - rhsIntegral;

fprintf('Final state [q, p] = [%g, %g]\n', x(1,end), x(2,end));
fprintf('Initial energy     = %.8f\n', Hvals(1));
fprintf('Final energy       = %.8f\n', Hvals(end));
fprintf('Energy balance residual = %.4e\n', residual);

% Optional CSV export
M = [t(:), x(1,:).', x(2,:).', Hvals(:), uvals(:), yvals(:), dvals(:), supplyVals(:)];
writematrix(M, 'Chapter18_Lesson3_matlab_output.csv');

figure;
plot(t, x(1,:), t, x(2,:));
grid on;
xlabel('Time [s]');
ylabel('States');
legend('q(t)', 'p(t)');
title('Port-Hamiltonian states');

figure;
plot(t, Hvals, t, Hvals(1) + cumtrapz(t, rhs), '--');
grid on;
xlabel('Time [s]');
ylabel('Energy');
legend('H(t)', 'H(0)+int(-diss+y*u)dt');
title('Energy balance verification');

% Simulink note:
% A Simulink model can be built with two Integrator blocks (q, p),
% a MATLAB Function block implementing pH dynamics, and scopes for q, p, H.
      

Chapter18_Lesson3.nb

(* Chapter18_Lesson3.nb *)
m = 1.5; k = 12.0; d = 0.8; omega = 1.4;
u[t_] := Sin[omega t];
H[q_, p_] := 1/2 k q^2 + p^2/(2 m);
y[q_, p_] := p/m;
diss[q_, p_] := d (p/m)^2;

eqns = {
   q'[t] == p[t]/m,
   p'[t] == -k q[t] - d (p[t]/m) + u[t],
   q[0] == 0.15,
   p[0] == 0.0
};

tmax = 20;
sol = NDSolveValue[eqns, {q, p}, {t, 0, tmax}];
qf[t_] := sol[[1]][t];
pf[t_] := sol[[2]][t];

energy[t_] := H[qf[t], pf[t]];
rhs[t_] := -diss[qf[t], pf[t]] + y[qf[t], pf[t]] u[t];

energyChange = energy[tmax] - energy[0];
rhsIntegral = NIntegrate[rhs[t], {t, 0, tmax}];
residual = N[energyChange - rhsIntegral];

Print["Energy balance residual = ", residual];

Plot[{qf[t], pf[t]}, {t, 0, tmax},
 PlotLegends -> {"q(t)", "p(t)"},
 AxesLabel -> {"t", "states"},
 PlotLabel -> "Port-Hamiltonian states"];

Plot[{energy[t], energy[0] + NIntegrate[rhs[s], {s, 0, t}]}, {t, 0, tmax},
 PlotLegends -> {"H(t)", "H(0)+int(-diss+y*u)dt"},
 AxesLabel -> {"t", "energy"},
 PlotLabel -> "Energy balance verification"];
      

Problem 1 (Legendre transform for a 1-DOF oscillator): Consider the undamped mass–spring system with \( L(q,\dot{q})=\frac{1}{2}m\dot{q}^2-\frac{1}{2}kq^2 \). Derive the Hamiltonian and Hamilton’s equations.

Solution:

The conjugate momentum is \( p=\frac{\partial L}{\partial \dot{q}}=m\dot{q} \), so \( \dot{q}=p/m \). The Hamiltonian is

\[ H(q,p)=p\dot{q}-L = p\left(\frac{p}{m}\right) - \left(\frac{1}{2}m\left(\frac{p}{m}\right)^2 - \frac{1}{2}kq^2\right) = \frac{p^2}{2m}+\frac{1}{2}kq^2. \]

Hamilton’s equations are

\[ \dot{q}=\frac{\partial H}{\partial p}=\frac{p}{m},\qquad \dot{p}=-\frac{\partial H}{\partial q}=-kq. \]

Eliminating \( p \) gives \( m\ddot{q}+kq=0 \), which matches the original dynamics.

Problem 2 (Conservation of Hamiltonian): Show that for an autonomous Hamiltonian system \( \dot{x}=\mathbf{J}\nabla H(x) \) with \( \mathbf{J}^\top=-\mathbf{J} \), the Hamiltonian is constant along trajectories.

Solution:

\[ \dot{H} = \nabla H^\top \dot{x} = \nabla H^\top \mathbf{J}\nabla H. \]

Let \( s=\nabla H^\top \mathbf{J}\nabla H \). Since it is a scalar, \( s=s^\top \). But

\[ s^\top = \nabla H^\top \mathbf{J}^\top \nabla H = \nabla H^\top (-\mathbf{J}) \nabla H = -s. \]

Hence \( s=-s \Rightarrow s=0 \), and therefore \( \dot{H}=0 \).

Problem 3 (Passivity of the pH model): For the pH system \( \dot{x}=(\mathbf{J}-\mathbf{R})\nabla H+\mathbf{g}u,\; y=\mathbf{g}^\top \nabla H \), prove that \( \dot{H}\le y^\top u \) if \( \mathbf{J}^\top=-\mathbf{J} \) and \( \mathbf{R}=\mathbf{R}^\top \succeq 0 \).

Solution:

\[ \dot{H} = \nabla H^\top \dot{x} = \nabla H^\top(\mathbf{J}-\mathbf{R})\nabla H + \nabla H^\top \mathbf{g}u. \]

The skew term is zero, and the port term is \( y^\top u \). Thus

\[ \dot{H} = -\nabla H^\top \mathbf{R}\nabla H + y^\top u. \]

Since \( \mathbf{R}\succeq 0 \), \( \nabla H^\top \mathbf{R}\nabla H \ge 0 \), so

\[ \dot{H}\le y^\top u. \]

This is the pH passivity inequality.

Problem 4 (RLC pH realization): For a series RLC circuit with source voltage \( u \), use states \( \phi \) (inductor flux) and \( q_c \) (capacitor charge) to write a pH model.

Solution:

Define the Hamiltonian \( H=\frac{\phi^2}{2L}+\frac{q_c^2}{2C} \). Then \( \nabla H=[\phi/L,\; q_c/C]^\top=[i,\; v_C]^\top \). Choose

\[ \mathbf{J}=\begin{bmatrix}0 & -1\\ 1 & 0\end{bmatrix},\quad \mathbf{R}=\begin{bmatrix}R_e & 0\\ 0 & 0\end{bmatrix},\quad \mathbf{g}=\begin{bmatrix}1\\0\end{bmatrix}. \]

Then

\[ \dot{x}=(\mathbf{J}-\mathbf{R})\nabla H+\mathbf{g}u \]

yields

\[ \dot{\phi}=u-R_e i-v_C,\qquad \dot{q}_c=i, \qquad y=i. \]

This exactly matches the series RLC dynamics and makes the source power explicit as \( ui \).

Problem 5 (Mechanical pH energy balance with zero input): For the mass–spring–damper pH model in Section 5.1 with \( u=0 \), derive \( \dot{H} \) and interpret the result.

Solution:

With \( H=\frac{1}{2}kq^2+\frac{p^2}{2m} \), \( \nabla H=[kq,\; p/m]^\top \), and \( \mathbf{R}=\mathrm{diag}(0,d) \),

\[ \dot{H} = -\nabla H^\top \mathbf{R}\nabla H = -d\left(\frac{p}{m}\right)^2. \]

Because \( d\ge 0 \), we have \( \dot{H}\le 0 \). The stored energy decreases monotonically and the loss term is exactly the mechanical dissipation power in the damper.

We transformed the energy-based Lagrangian viewpoint into the Hamiltonian formulation by introducing conjugate momentum and the Legendre transform, obtained Hamilton’s canonical equations, and proved energy conservation for autonomous conservative systems. We then extended the framework to port-Hamiltonian systems, where interconnection, dissipation, and external power ports are encoded explicitly. This provides a unified, physically transparent modeling template for mechanical, electrical, and other engineering domains, and it prepares the ground for systematic interconnection methods in the next lesson (bond graphs) and later multi-domain models.