$\def\rmHS{{\rm HS}}$
$\def\rmFe{{\rm Fe}}$
$\def\rmm{{\rm m}}$
$\def\vr{\vec{r}}$
$\def\vH{\vec{H}}$
$\def\vB{\vec{B}}$
$\def\vA{\vec{A}}$
$\def\nR{{\mathbb{R}}}$
$\def\ph{{\varphi}}$
$\def\grad{\mathrm{grad}}$
$\def\div{\mathrm{div}}$
A simple example demonstrates how it works with concentrated parameters.
We discuss a horseshoe magnet standing on an iron bar (see the following Figure).
We assume that the field within the horseshoe magnet is almost
homogeneous in longitudinal direction and denote the longitudinal
field strength with
$H_\rmHS$
and the flux density with $B_\rmHS$
(that works quite good if the thickness of the arc segment is small
enough).
The field in the magnet can be described by
$$
B_\rmHS = \mu_0 \left(H_\rmHS + M_\rmHS\right)
$$
with the magnetization $M_\rmHS$ of the horseshoe magnet.
Thereby, we assume that the changes in magnetization displacement on
the iron bar are so small that $M_\rmHS$ can be approximated as a
linear-affine function of the field strength:
$$
M_\rmHS = M_{\rmHS0} + \chi H_\rmHS
$$
We define the flux $\Phi:=A_\rmHS B_\rmHS$ leaving the horseshoe
at $x_2$ and the magnetic voltage drop as the integral of
$H_\rmHS$ from $x_1$ to $x_2$ over the mean path inside the
horseshoe $V_\rmm:= l_\rmHS H_\rmHS$. Substituting this into the last two formulae gives
$$
\Phi = A_\rmHS B_\rmHS = A_\rmHS\mu_0\left((1+\chi)\frac{V_\rmm}{l_\rmHS} + M_{\rmHS0}\right) = G_\rmHS V_\rmm + \Phi_{\rmHS0}
$$
The magnet acts at its poles like a parallel connection of a permeance $G_\rmHS=\frac{(1+\chi)\mu_0A}{l_\rmHS}$ with a flux source $\Phi_0=\mu_0AM_{\rmHS0}$.
On the other side the flux goes into the iron bar (we neglect the
stray flux here). The permeance $G_\rmFe(x)$ of the iron bar depends on
the position of the magnet on the bar.
With the magnetic voltage drop $-V_\rmm$ from $x_2$ to $x_1$ on the iron permeance we get the equation
$$
G_\rmHS V_\rmm + \Phi_0 = \Phi = G_\rmFe(x)\cdot(-V_\rmm)
$$
for the flux balance or equivalently
$$
V_\rmm = \frac{\Phi_0}{G_\rmHS+G_\rmFe(x)}.
$$
The stored free magnetic energy is
$$
W_\rmm = \frac12\int_{\nR^n} \vB\cdot\vH d V.
$$
We neglect the stray flux and its energy. Therefore, the magnetic
energy can be decomposed into the part $W_\rmHS$ stored in the
horseshoe magnet and the part $W_\rmFe$ stored in the iron.
In the magnet we get
$$
W_\rmHS = \frac12 A_\rmHS l_\rmHS B_\rmHS H_\rmHS = \frac12 V_\rmm \Phi
$$
For the iron material we get similarly
$$
W_\rmFe = \frac12 V_\rmm\Phi
$$
and the overall energy can be calculated in the following way:
$W_\rmm = \frac12 V_\rmm\Phi + \frac12 V_\rmm\Phi$
$\phantom{W_\rmm}= \frac12 V_\rmm(G_\rmHS V_\rmm + \Phi_0) + \frac12 V_\rmm G_\rmFe(x) V_\rmm$
$\phantom{W_\rmm}= \frac12\left(G_\rmHS+G_\rmFe(x)\right)V_\rmm^2 + \frac12 V_\rmm\Phi_0$
and with $V_\rmm = \frac{\Phi_0}{G_\rmHS+G_\rmFe(x)}$:
$$
W_\rmm = \frac{\Phi_0^2}{G_\rmHS+G_\rmFe(x)}.
$$
And with this formula the magnetic energy has a strong local minimum at a strong local maximum of $G_\rmFe(x)$.
There remains to discuss the magnetic energy within the iron block.
The field within the iron block can be calculated with the help of the
magneto-static potential $\ph_\rmm$ with $H_\rmFe = -\grad\ph_\rmm$
defined through the equation
$$
0=\div(\vB)=-\div(\mu\grad\ph_\rmm)
$$
and the boundary conditions:
$\ph_\rmm(\vr_1)=0$ for points $\vr_1$ on the pole 1,
$\ph_\rmm(\vr_2)=V_\rmm$ for points $\vr_2$ on pole 2 and
for points $\vr_2$ on pole 2 and
$\grad\ph(\vr) \cdot d\vA=0$ at all the other boundary points $\vr$ whereby $d\vA$ is the outwards directed area element.
The stored energy is
$$
W_\rmFe = \int_{V_\rmFe} \vH\cdot\vB dV
= \int_{V_\rmFe} (\grad\ph_\rmm)\cdot(\mu\grad\ph_\rmm) dV
$$
With the differentiation rule $\div(\ph_\rmm\mu\grad\ph_\rmm) = (\grad\ph_\rmm)\cdot(\mu\grad\ph_\rmm) + \ph_\rmm \div(\mu\grad\ph_\rmm)$ and Gauß's theorem for the divergence operator one obtains
$$
W_\rmFe = \oint_{\partial V_\rmFe} \ph_\rmm \mu (\grad\ph_\rmm) \cdot d\vA - \int_{V_\rmFe} \ph_\rmm \underbrace{\div(\mu\grad\ph_\rmm)}_{=0}d V
$$
Thereby, the volume integral is zero because of $\div(\mu\grad\ph_\rmm)=0$.
The surface integral has only a nonzero component at pole 2 since at
pole 1 $\ph_\rmm=0$ and on the complement surface $\grad\ph_\rmm \cdot
d \vA=0$. Furthermore, on pole 2 the magnetic potential
$\ph_\rmm=V_\rmm$ is assumed to be constant. Thus, we have
$$
W_\rmFe = (-V_\rmm) \int_{A_2} \mu (-\grad\ph_\rmm)\cdot d\vA = (-V_\rmm) \int_{A_2} \vB \cdot d\vA = - V_\rmm \Phi.
$$