This chapter shows how the complete port‑Hamiltonian network is assembled into a single nonlinear system and solved at each audio sample. The assembly is fully automatic given the component–port–connection graph; the implementer only needs to provide the per‑component residual functions.
All state variables are concatenated into one vector \(\mathbf{X}\):
$$ \mathbf{X} = \bigl( \underbrace{m_i, p_{\text{mom},i}, E_{\text{tot},i}}_{\text{pipe cells}},\; \underbrace{m_{\text{vol}}, U_{\text{vol}}}_{\text{volumes}},\; \underbrace{\lambda_k}_{\text{Lagrange multipliers}} \bigr)^T. $$For an implicit midpoint step from \(t^n\) to \(t^{n+1}\), we seek \(\mathbf{X}^{n+1}\) such that
$$ \mathbf{F}(\mathbf{X}^{n+1}) = \mathbf{X}^{n+1} - \mathbf{X}^n - \Delta t\, (\mathbf{J} - \mathbf{R})\,\nabla H\!\left(\tfrac{\mathbf{X}^{n+1} + \mathbf{X}^n}{2}\right) - \Delta t\,\mathbf{B}\,\mathbf{u}_{\text{ext}}^{n+1/2} \;+\; \mathbf{C}^T \boldsymbol{\lambda} = \mathbf{0}, $$together with the algebraic constraints \(\mathbf{C}(\mathbf{X}^{n+1}) = \mathbf{0}\). Here \(\mathbf{C}\) collects the junction equations (pressure equality and net mass flow zero). The Lagrange multiplier term \(\mathbf{C}^T \boldsymbol{\lambda}\) adds the port constraint forces into the dynamic equations.
For each component we implement a method that computes its local residual and adds it to the global vector. The implicit midpoint rule requires evaluation of the internal fluxes at the midpoint state \(\mathbf{X}^{n+1/2} = \tfrac12(\mathbf{X}^n + \mathbf{X}^{n+1})\), so the local residual functions receive the midpoint state as input.
The global Jacobian \(\mathbf{J}_F = \partial\mathbf{F} / \partial\mathbf{X}\) is block‑sparse:
This structure is known at graph construction time, allowing a static sparse matrix format with pre‑allocated non‑zero entries for fast Jacobian assembly.
The nonlinear system \(\mathbf{F}(\mathbf{X}) = \mathbf{0}\) is solved with an inexact Newton method using a preconditioned Krylov solver for the linear sub‑problems. The high‑level loop is:
// Initial guess: previous solution (or extrapolated)
X = X_old;
for (int iter = 0; iter < maxNewtonIter; iter++)
{
// 1. Evaluate global residual F at current X (midpoint already inside F)
double[] F = EvaluateResidual(X);
if (Norm(F) < tolerance) break; // converged
// 2. Assemble Jacobian J = dF/dX (analytical or finite‑difference)
SparseMatrix J = AssembleJacobian(X);
// 3. Solve J ΔX = -F approximately using GMRES with ILU preconditioner
double[] dX = GMRES(J, -F, ILU(J), maxGMRESIter, gmresTol);
// 4. Line search to enforce global convergence (optional)
double alpha = 1.0;
while (Norm(EvaluateResidual(X + alpha*dX)) > (1 - 1e-4*alpha) * Norm(F))
alpha *= 0.5;
X += alpha * dX;
}
// X now holds the converged X^{n+1}
A typical setup uses maxNewtonIter = 4, tolerance = 1e‑8, and a restarted GMRES(30) with ILU(0) preconditioner.
The block‑banded structure allows a simple but effective ILU(0) preconditioner. For pipes, the diagonal \(3\times3\) blocks are always non‑singular as long as the mass is positive. Alternative: block‑Jacobi, where each cell’s \(3\times3\) block plus its immediate neighbours is solved exactly.
class PHSWorld {
List<Component> components;
List<Connection> connections;
double[] X, F, dX;
SparseMatrix J;
void BuildStateVector() {
int offset = 0;
foreach (var c in components) {
c.stateOffset = offset;
c.PackState(X, offset);
offset += c.stateSize;
}
// Lagrange multipliers appended after all component states
lagrangeOffset = offset;
}
double[] EvaluateResidual(double[] Xk) {
Array.Clear(F);
foreach (var c in components)
c.AddResidual(Xk, F, c.stateOffset); // dynamic equations
foreach (var conn in connections)
conn.AddConstraintResidual(Xk, F, lagrangeOffset); // C(X)=0
return F;
}
SparseMatrix AssembleJacobian(double[] Xk) {
J.Zero();
foreach (var c in components)
c.AddJacobianBlocks(Xk, J);
foreach (var conn in connections)
conn.AddConstraintJacobianBlocks(Xk, J);
return J;
}
void Step() {
BuildStateVector();
NewtonKrylov();
UnpackState();
}
}