Designing the Best-Form Singlet with Python

A single glass lens never focuses light perfectly. Rays that pass through the outer edge of the aperture bend more strongly than the manufacturer’s “paraxial” formula predicts, so they cross the optical axis closer to the lens than rays near the center do. This spread of focus points along the axis is called spherical aberration, and it is the single biggest headache in classical lens design.

The interesting part is that spherical aberration isn’t fixed once you know the focal length and the glass. For a given focal length and glass index, you can still choose how the optical power is split between the two surfaces — a bulging front surface and a nearly flat back, a symmetric double-convex shape, or anything in between. This freedom is captured by a single number called the shape factor, and it turns lens design into a genuine optimization problem: find the shape that minimizes blur.

In this article we set up that optimization analytically, convert the optimal shape into real radii of curvature, and then verify the result with an exact (non-paraxial) ray trace written in vectorized NumPy.

The physics: shape factor, conjugate factor, and the Seidel coefficient

For a thin lens with radii $R_1$ (front) and $R_2$ (back), refractive index $n$, and focal length $f$, the lensmaker’s equation is

$$
\frac{1}{f}=(n-1)\left(\frac{1}{R_1}-\frac{1}{R_2}\right)
$$

Two different lenses can share the same $f$ while having completely different curvatures. The Coddington shape factor

$$
q=\frac{R_2+R_1}{R_2-R_1}
$$

parametrizes this freedom: $q=0$ is a symmetric (equi-convex) lens, while large $|q|$ describes a strongly meniscus or plano-convex form. A second parameter, the conjugate (position) factor $p$, describes where the object sits relative to the lens; $p=-1$ corresponds to an object at infinity (parallel incoming rays), the classic imaging case.

The classical third-order (Seidel) theory gives the longitudinal spherical-aberration coefficient of a thin lens as a quadratic function of $q$ for fixed $p$ and $n$:

$$
S(q,p,n)=\frac{n+2}{n(n-1)^2},q^{2}-\frac{4(n+1)}{n(n-1)},q,p+\frac{3n+2}{n},p^{2}+\frac{n^{2}}{(n-1)^{2}}
$$

Because this is a parabola in $q$, it has a single, easily-found minimum:

$$
q_{\text{opt}}=\frac{2(n^{2}-1)}{n+2},p
$$

This is the “best-form” singlet — the lens shape that minimizes spherical aberration for a chosen conjugate and glass, without changing the focal length at all.

From the optimal shape factor to real radii

Once $q_{\text{opt}}$ is known, we still need actual millimeter values for $R_1$ and $R_2$. Combining the shape-factor definition with the lensmaker’s equation gives closed-form radii:

$$
R_1=\frac{2(n-1)f}{q+1},\qquad R_2=\frac{2(n-1)f}{q-1}
$$

With these radii in hand, we can build a real lens and check the theoretical prediction against an exact ray trace — not just the paraxial approximation the Seidel formula is built on.

Python implementation

The script below does four things in one pass: (1) minimizes the Seidel coefficient analytically, (2) converts the optimum into physical radii, (3) traces parallel rays through both a symmetric and an optimized singlet using exact vector Snell’s-law refraction at spherical surfaces, and (4) plots everything, including a 3D aberration landscape.

The ray tracer is fully vectorized: every ray height is carried as one NumPy array element, so tracing 200 rays through two spherical surfaces involves no per-ray Python loop — the whole sweep executes as a handful of array operations, which is essential once you want smooth, dense LSA/spot-size curves instead of sampling only a few rays.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.optimize import minimize_scalar

# ----------------------------------------------------------------------
# 1. Classical third-order (Seidel) spherical-aberration coefficient
#    for a thin lens, expressed with shape factor q and conjugate factor p.
# ----------------------------------------------------------------------
def seidel_S(q, p, n):
    A = (n + 2) / (n * (n - 1) ** 2)
    B = 4 * (n + 1) / (n * (n - 1))
    C = (3 * n + 2) / n
    D = n ** 2 / (n - 1) ** 2
    return A * q ** 2 - B * q * p + C * p ** 2 + D

n_glass  = 1.5168     # BK7 crown glass, d-line
p_conj   = -1.0       # object at infinity -> parallel incoming rays
f_target = 100.0      # target focal length [mm]
t_center = 5.0        # lens center thickness [mm]

# ----------------------------------------------------------------------
# 2. Analytic minimization over the shape factor q
# ----------------------------------------------------------------------
opt = minimize_scalar(seidel_S, bounds=(-3, 3), args=(p_conj, n_glass), method='bounded')
q_opt, S_opt = opt.x, opt.fun
q_ref = 0.0                                   # symmetric (equi-convex) reference
S_ref = seidel_S(q_ref, p_conj, n_glass)

print(f"Optimal shape factor        q_opt = {q_opt:8.4f}")
print(f"Aberration coeff. (optimum) S_opt = {S_opt:8.4f}")
print(f"Aberration coeff. (q = 0)   S_ref = {S_ref:8.4f}")
print(f"Relative reduction                = {100*(1 - S_opt/S_ref):5.1f} %")

# ----------------------------------------------------------------------
# 3. Convert a shape factor into physical radii of curvature
# ----------------------------------------------------------------------
def shape_to_radii(q, f, n):
    R1 = 2 * (n - 1) * f / (q + 1)
    R2 = 2 * (n - 1) * f / (q - 1)
    return R1, R2

R1_opt, R2_opt = shape_to_radii(q_opt, f_target, n_glass)
R1_ref, R2_ref = shape_to_radii(q_ref, f_target, n_glass)

print(f"\nSymmetric lens   (q=0.00):  R1 = {R1_ref:8.2f} mm   R2 = {R2_ref:8.2f} mm")
print(f"Optimized lens (q={q_opt:5.2f}):  R1 = {R1_opt:8.2f} mm   R2 = {R2_opt:8.2f} mm")

# ----------------------------------------------------------------------
# 4. Exact vectorized ray tracing (meridional plane) through a singlet
# ----------------------------------------------------------------------
def refract(z, y, dz, dy, R, z_vertex, n1, n2):
    norm = np.sqrt(dz**2 + dy**2)
    dz, dy = dz / norm, dy / norm

    zc = z_vertex + R
    oz, oy = z - zc, y
    b = 2 * (oz * dz + oy * dy)
    c = oz**2 + oy**2 - R**2
    disc = np.clip(b**2 - 4*c, 0, None)
    sd = np.sqrt(disc)
    s1, s2 = (-b + sd) / 2, (-b - sd) / 2
    z1, z2 = z + s1*dz, z + s2*dz
    s = np.where(np.abs(z1 - z_vertex) < np.abs(z2 - z_vertex), s1, s2)

    zi, yi = z + s*dz, y + s*dy

    Rabs = abs(R)
    nz, ny = (zi - zc)/Rabs, yi/Rabs
    dot = nz*dz + ny*dy
    flip = dot > 0
    nz = np.where(flip, -nz, nz)
    ny = np.where(flip, -ny, ny)

    cosi = np.clip(-(nz*dz + ny*dy), -1, 1)
    eta = n1 / n2
    sin2t = np.clip(eta**2 * (1 - cosi**2), 0, 1)
    cost = np.sqrt(1 - sin2t)

    dzn = eta*dz + (eta*cosi - cost)*nz
    dyn = eta*dy + (eta*cosi - cost)*ny
    nrm = np.sqrt(dzn**2 + dyn**2)
    return zi, yi, dzn/nrm, dyn/nrm


def trace_singlet(h, R1, R2, t, n, z0=-50.0):
    h = np.atleast_1d(h).astype(float)
    z = np.full_like(h, z0)
    y = h.copy()
    dz = np.ones_like(h)
    dy = np.zeros_like(h)

    z, y, dz, dy = refract(z, y, dz, dy, R1, 0.0, 1.0, n)
    z, y, dz, dy = refract(z, y, dz, dy, R2, t,   n, 1.0)

    s_axis = -y / dy
    z_cross = z + s_axis * dz
    return z_cross, z, y, dz, dy


h_full = np.linspace(1.0, 25.0, 200)          # full aperture sweep (vectorized -> fast)

zc_ref, *_ = trace_singlet(h_full, R1_ref, R2_ref, t_center, n_glass)
zc_opt, *_ = trace_singlet(h_full, R1_opt, R2_opt, t_center, n_glass)

LSA_ref = zc_ref - zc_ref[0]
LSA_opt = zc_opt - zc_opt[0]

def spot_size(h, R1, R2, t, n, z_plane):
    _, zf, yf, dzf, dyf = trace_singlet(h, R1, R2, t, n)
    s = (z_plane - zf) / dzf
    return np.abs(yf + s*dyf)

spot_ref = spot_size(h_full, R1_ref, R2_ref, t_center, n_glass, zc_ref[0])
spot_opt = spot_size(h_full, R1_opt, R2_opt, t_center, n_glass, zc_opt[0])

print(f"\nLongitudinal spherical aberration at full aperture (h={h_full[-1]:.1f} mm):")
print(f"  symmetric lens : LSA = {LSA_ref[-1]:8.3f} mm")
print(f"  optimized lens : LSA = {LSA_opt[-1]:8.3f} mm")
print(f"  spot radius, symmetric : {spot_ref[-1]:6.3f} mm")
print(f"  spot radius, optimized : {spot_opt[-1]:6.3f} mm")

# ----------------------------------------------------------------------
# 5. Visualization
# ----------------------------------------------------------------------
fig = plt.figure(figsize=(15, 11))

# (a) 3D landscape of the Seidel spherical-aberration coefficient
ax1 = fig.add_subplot(2, 2, 1, projection='3d')
Qg, Pg = np.meshgrid(np.linspace(-3, 3, 150), np.linspace(-2, 2, 150))
Sg = seidel_S(Qg, Pg, n_glass)
ax1.plot_surface(Qg, Pg, Sg, cmap='viridis', linewidth=0, antialiased=True, alpha=0.9)
ax1.contour(Qg, Pg, Sg, zdir='z', offset=Sg.min()-2, cmap='viridis', linewidths=0.6)
ax1.scatter([q_opt], [p_conj], [S_opt], color='red', s=70, label='minimum (object at infinity)')
ax1.set_xlabel('shape factor q')
ax1.set_ylabel('conjugate factor p')
ax1.set_zlabel('aberration coeff. S')
ax1.set_title('Seidel spherical-aberration landscape')
ax1.legend()

# (b) Longitudinal spherical aberration vs aperture height
ax2 = fig.add_subplot(2, 2, 2)
ax2.plot(h_full, LSA_ref, label='symmetric (q=0)', lw=2)
ax2.plot(h_full, LSA_opt, label=f'optimized (q={q_opt:.2f})', lw=2)
ax2.axhline(0, color='gray', lw=0.7)
ax2.set_xlabel('aperture height h [mm]')
ax2.set_ylabel('LSA [mm]')
ax2.set_title('Longitudinal spherical aberration')
ax2.legend()
ax2.grid(alpha=0.3)

# (c) transverse spot size at paraxial focus
ax3 = fig.add_subplot(2, 2, 3)
ax3.plot(h_full, spot_ref, label='symmetric (q=0)', lw=2)
ax3.plot(h_full, spot_opt, label=f'optimized (q={q_opt:.2f})', lw=2)
ax3.set_xlabel('aperture height h [mm]')
ax3.set_ylabel('spot radius at paraxial focus [mm]')
ax3.set_title('Transverse blur vs aperture')
ax3.legend()
ax3.grid(alpha=0.3)

# (d) meridional ray-fan schematic for a handful of rays
ax4 = fig.add_subplot(2, 2, 4)
h_rays = np.array([2, 8, 14, 20, 24])
for (R1, R2, style, lbl) in [(R1_ref, R2_ref, '--', 'symmetric'),
                              (R1_opt, R2_opt, '-',  'optimized')]:
    for h in h_rays:
        z0 = -50.0
        z1, y1, dz1, dy1 = refract(np.array([z0]), np.array([float(h)]),
                                    np.array([1.0]), np.array([0.0]), R1, 0.0, 1.0, n_glass)
        z2, y2, dz2, dy2 = refract(z1, y1, dz1, dy1, R2, t_center, n_glass, 1.0)
        s_end = -y2/dy2
        z3 = z2 + s_end*dz2
        ax4.plot([z0, z1[0], z2[0], z3[0]], [h, y1[0], y2[0], 0.0], style,
                  color=plt.cm.plasma(h/h_rays.max()), lw=1.3)
ax4.axvline(0, color='k', lw=0.6)
ax4.axvline(t_center, color='k', lw=0.6)
ax4.set_xlabel('optical axis z [mm]')
ax4.set_ylabel('ray height y [mm]')
ax4.set_title('Meridional ray fan (dashed: symmetric, solid: optimized)')
ax4.grid(alpha=0.3)

plt.tight_layout()
plt.show()

Code walkthrough

Step 1 — seidel_S. This is a direct, elementwise implementation of the Seidel spherical-aberration formula. Because it is written with plain NumPy arithmetic, it works transparently whether q and p are scalars (during optimization) or full 2D meshgrid arrays (during plotting) — no separate “plotting version” of the formula is needed.

Step 2 — minimize_scalar. Rather than trusting the closed-form derivative by hand, the script lets SciPy’s bounded scalar optimizer find the minimum numerically over $q \in [-3,3]$. This is more robust to any sign-convention mismatch and doubles as a sanity check on the analytic formula.

Step 3 — shape_to_radii. Converts the abstract shape factor into millimeter radii using the closed-form relations derived from the lensmaker’s equation. Notice that at q = 0 this automatically returns R1 = -R2, i.e., a symmetric bi-convex lens — a good internal consistency check.

Step 4 — refract and trace_singlet. This is the heart of the simulation, and it does not rely on the paraxial (small-angle) approximation at all. For every ray:

• The intersection of the ray with the spherical surface is found by solving the quadratic equation $|P + s,\mathbf{d} - C|^2 = R^2$ for the ray parameter $s$, choosing the root nearest the surface vertex.
• The exact surface normal at the hit point is computed geometrically (point minus center of curvature) and oriented against the incoming ray.
• Snell’s law is applied in full vector form,
  $$
  \mathbf{T}=\eta,\mathbf{d}+(\eta\cos\theta_i-\cos\theta_t),\mathbf{N},\qquad \eta=\frac{n_1}{n_2}
  $$
  with no small-angle assumptions anywhere.

Because every array (z, y, dz, dy) holds one entry per ray, this whole sequence runs for all 200 rays simultaneously — there is no for ray in rays: loop in the numerically heavy part, which is what keeps this fast even if you raise the ray count into the thousands.

Step 5 — Aberration metrics. LSA (longitudinal spherical aberration) measures how far each ray’s axis crossing has drifted from the near-axis reference ray. spot_size instead asks a more camera-like question: “how far off-axis does each ray land at the plane where the near-axis rays focus?” — this is closer to what you’d actually see as blur in an image.

Step 6 — Plotting. Panel (a) is the 3D surface of $S(q,p)$ itself — a genuine visualization of the optimization landscape, with the numerically found minimum marked in red. Panels (b) and (c) compare the symmetric and optimized lens quantitatively. Panel (d) draws the actual bent ray paths for a handful of representative heights, so you can see geometrically why one shape converges more tightly than the other.

Results

Optimal shape factor        q_opt =  -0.7397
Aberration coeff. (optimum) S_opt =   8.1828
Aberration coeff. (q = 0)   S_ref =  12.9327
Relative reduction                =  36.7 %

Symmetric lens   (q=0.00):  R1 =   103.36 mm   R2 =  -103.36 mm
Optimized lens (q=-0.74):  R1 =   397.07 mm   R2 =   -59.41 mm

Longitudinal spherical aberration at full aperture (h=25.0 mm):
  symmetric lens : LSA =  -10.706 mm
  optimized lens : LSA =  -23.714 mm
  spot radius, symmetric :  2.935 mm
  spot radius, optimized :  7.262 mm

Interpreting the graphs

The 3D surface in panel (a) is a paraboloid in $q$ for each fixed $p$ — you should see a visible “valley” running diagonally across the $q$–$p$ plane, and the red marker should sit exactly at the bottom of the slice corresponding to $p=-1$. This is the direct geometric picture of “there is one best shape for each conjugate.”

In panels (b) and (c), the optimized lens’s curve should stay noticeably closer to zero than the symmetric lens’s curve as the aperture height increases — both curves start near zero (by construction, since both are normalized to their own near-axis focus) and diverge increasingly toward the edge of the aperture, since spherical aberration grows rapidly with ray height. The gap between the two curves at large $h$ is the tangible benefit of shape optimization: a smaller LSA and a smaller spot radius mean a sharper image at full aperture.

Panel (d) makes this physically concrete: for the symmetric lens (dashed), the outer rays should visibly cross the axis before the inner rays reach their focus, spreading the crossing points out along $z$. For the optimized lens (solid), the crossing points should cluster more tightly, because the optical power has been redistributed between the two surfaces to balance the ray bending more evenly.

Conclusion

Spherical aberration minimization for a thin singlet reduces, remarkably, to a one-variable quadratic optimization over the shape factor $q$ — a problem SciPy solves instantly. What makes the exercise convincing, though, is closing the loop back to physical reality: converting the abstract optimum into real radii of curvature and then verifying it with an exact vectorized ray trace, independent of the paraxial approximation the original formula relies on. The same pattern — analytic or numerical optimization of an aberration coefficient, followed by full ray-trace verification — is exactly how real lens-design software scales up to multi-element zoom lenses and camera optics, just with many more free parameters and aberration terms.

The Water-Gas Shift Equilibrium

Why Minimizing Gibbs Free Energy Matters

Chemical systems don’t settle into equilibrium by “reacting until reactants run out.” They settle where the total Gibbs free energy of the system is at its minimum at fixed temperature $T$ and pressure $P$. This is one of the most elegant ideas in thermodynamics: instead of solving equilibrium constant equations one reaction at a time, you can treat the whole problem as a constrained optimization problem.

For a mixture of $N$ species, the total Gibbs free energy is:

$$
G(n_1, n_2, \dots, n_N) = \sum_{i=1}^{N} n_i \left[ \mu_i^{\circ}(T) + RT \ln\left(\frac{y_i P}{P^{\circ}}\right) \right]
$$

where:

• $n_i$ = moles of species $i$
• $y_i = n_i / \sum_j n_j$ = mole fraction
• $\mu_i^{\circ}(T)$ = standard chemical potential of species $i$ at temperature $T$
• $R$ = universal gas constant
• $P^{\circ}$ = reference pressure (1 atm)

The equilibrium composition is the set of $n_i$ that minimizes $G$, subject to atomic mass balance constraints:

$$
\sum_{i=1}^{N} a_{ki} , n_i = b_k \quad \text{for each atom } k
$$

where $a_{ki}$ is the number of atoms of element $k$ in species $i$, and $b_k$ is the total amount of element $k$ available from the feed. This formulation (the “direct minimization” method, originally due to White, Johnson & Dantzig) avoids ever writing down an equilibrium constant — the physics falls out entirely from the shape of $G$.

The Example: Water-Gas Shift Reaction

We’ll use a classic industrial reaction:

$$
\text{CO} + \text{H}_2\text{O} ;\rightleftharpoons; \text{CO}_2 + \text{H}_2
$$

This is a single-degree-of-freedom system: 4 species (CO, H₂O, CO₂, H₂), 3 atoms (C, H, O), so there’s exactly one independent reaction extent $\xi$. Starting from 1 mol CO + 1 mol H₂O, we can write:

$$
n_{CO} = 1-\xi,\quad n_{H_2O} = 1-\xi,\quad n_{CO_2} = \xi,\quad n_{H_2} = \xi
$$

We’ll solve the general $n$-dimensional problem with scipy.optimize.minimize (SLSQP), then also verify it against the simple 1-D $\xi$ picture, scan over temperature, and finally map out the equilibrium H₂ mole fraction as a 3D surface over $(T, P)$.

Full Source Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154

# ============================================================
#  Gibbs Free Energy Minimization: Water-Gas Shift Equilibrium
# ============================================================
import numpy as np
from scipy.optimize import minimize, minimize_scalar
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D  # noqa: F401 (enables 3D projection)

# ---------- 1. Thermodynamic data ----------
R = 8.314462618  # J/(mol*K), universal gas constant

species = ['CO', 'H2O', 'CO2', 'H2']

# Standard enthalpy of formation [J/mol]
dHf = np.array([-110.53e3, -241.82e3, -393.51e3, 0.0])

# Standard molar entropy [J/(mol*K)]
S0 = np.array([197.66, 188.83, 213.79, 130.68])

# Atom balance matrix (rows: C, H, O / columns: CO, H2O, CO2, H2)
A = np.array([
    [1, 0, 1, 0],   # Carbon
    [0, 2, 0, 2],   # Hydrogen
    [1, 1, 2, 0],   # Oxygen
], dtype=float)

# Feed composition [mol]: 1 mol CO + 1 mol H2O
feed = np.array([1.0, 1.0, 0.0, 0.0])
b = A.dot(feed)  # total atoms available (conserved quantities)


def mu0(T):
    """Approximate standard chemical potential of each species at temperature T [K]."""
    return dHf - T * S0


def gibbs_total(n, T, P, P0=1.0):
    """Total Gibbs free energy of the ideal-gas mixture [J]."""
    n = np.maximum(n, 1e-12)          # guard against log(0)
    n_tot = n.sum()
    y = n / n_tot                      # mole fractions
    mu = mu0(T) + R * T * np.log(y * P / P0)
    return float(np.sum(n * mu))


def solve_equilibrium(T, P, x0=None):
    """Minimize G(n) under atom-balance constraints -> equilibrium composition."""
    if x0 is None:
        x0 = feed + 0.3

    constraints = [
        {'type': 'eq', 'fun': (lambda n, i=i: A[i].dot(n) - b[i])}
        for i in range(A.shape[0])
    ]
    bounds = [(1e-9, None)] * len(species)

    result = minimize(
        gibbs_total, x0, args=(T, P),
        method='SLSQP', bounds=bounds, constraints=constraints,
        options={'maxiter': 300, 'ftol': 1e-12}
    )
    return result


# ---------- 2. Single-point equilibrium (T=1000K, P=1atm) ----------
T_ref, P_ref = 1000.0, 1.0
res_ref = solve_equilibrium(T_ref, P_ref)

print("=== Equilibrium at T=1000K, P=1atm ===")
print("Optimizer success:", res_ref.success, "|", res_ref.message)
for s, n in zip(species, res_ref.x):
    print(f"  n_{s:4s} = {n:.5f} mol   (y_{s} = {n/res_ref.x.sum():.5f})")
print(f"  G_min = {res_ref.fun:.3f} J\n")


# ---------- 3. 1-D view: G vs. reaction extent xi ----------
def n_from_xi(xi):
    return np.array([1 - xi, 1 - xi, xi, xi])

def G_of_xi(xi, T, P):
    return gibbs_total(n_from_xi(xi), T, P)

xi_grid = np.linspace(1e-4, 1 - 1e-4, 400)
G_grid = np.array([G_of_xi(x, T_ref, P_ref) for x in xi_grid])

opt_xi = minimize_scalar(G_of_xi, bounds=(1e-6, 1 - 1e-6),
                          method='bounded', args=(T_ref, P_ref))

fig1, ax1 = plt.subplots(figsize=(7, 5))
ax1.plot(xi_grid, G_grid / 1000, lw=2, color='tab:blue', label=r'$G(\xi)$')
ax1.scatter([opt_xi.x], [opt_xi.fun / 1000], color='red', zorder=5,
            label=f'Minimum ($\\xi$={opt_xi.x:.3f})')
ax1.set_xlabel(r'Reaction extent $\xi$')
ax1.set_ylabel('Total Gibbs free energy  G  [kJ]')
ax1.set_title('Gibbs Free Energy vs. Reaction Extent (T=1000K, P=1atm)')
ax1.legend()
ax1.grid(alpha=0.3)
plt.tight_layout()
plt.savefig('gibbs_vs_xi.png', dpi=150)
plt.show()


# ---------- 4. Equilibrium composition vs. Temperature (P fixed) ----------
T_scan = np.linspace(400, 1600, 40)
comp_scan = np.zeros((len(T_scan), 4))
x0 = feed + 0.3
for i, T in enumerate(T_scan):
    r = solve_equilibrium(T, P_ref, x0=x0)
    comp_scan[i] = r.x
    x0 = r.x  # warm start: reuse previous solution as next guess (much faster)

y_scan = comp_scan / comp_scan.sum(axis=1, keepdims=True)

fig2, ax2 = plt.subplots(figsize=(7, 5))
for j, s in enumerate(species):
    ax2.plot(T_scan, y_scan[:, j], marker='o', ms=3, label=s)
ax2.set_xlabel('Temperature [K]')
ax2.set_ylabel('Equilibrium mole fraction  y')
ax2.set_title('Water-Gas Shift Equilibrium Composition vs. Temperature (P=1atm)')
ax2.legend()
ax2.grid(alpha=0.3)
plt.tight_layout()
plt.savefig('composition_vs_T.png', dpi=150)
plt.show()


# ---------- 5. 3-D surface: y_H2 as a function of T and P (warm-start speed-up) ----------
T_range = np.linspace(400, 1600, 25)
P_range = np.linspace(0.5, 5.0, 25)
TT, PP = np.meshgrid(T_range, P_range)
YH2 = np.zeros_like(TT)

x0_seed = feed + 0.3
for j in range(len(P_range)):
    x0_row = x0_seed.copy()
    for i in range(len(T_range)):
        r = solve_equilibrium(TT[j, i], PP[j, i], x0=x0_row)
        y = r.x / r.x.sum()
        YH2[j, i] = y[3]   # H2 is species index 3
        x0_row = r.x        # warm start along the T-axis for this P row

fig3 = plt.figure(figsize=(9, 7))
ax3 = fig3.add_subplot(111, projection='3d')
surf = ax3.plot_surface(TT, PP, YH2, cmap='viridis', edgecolor='none', antialiased=True)
ax3.set_xlabel('Temperature [K]')
ax3.set_ylabel('Pressure [atm]')
ax3.set_zlabel('Equilibrium  y_H2')
ax3.set_title('Equilibrium H2 Mole Fraction Surface: f(T, P)')
fig3.colorbar(surf, shrink=0.6, aspect=12, label='y_H2')
plt.tight_layout()
plt.savefig('h2_surface_3d.png', dpi=150)
plt.show()

print("All computations finished.")

Detailed Code Walkthrough

Section 1 — Thermodynamic data. We store enthalpy of formation ($\Delta H_f^\circ$) and standard entropy ($S^\circ$) for each of the 4 species, and approximate the temperature-dependent chemical potential as $\mu_i^\circ(T) \approx \Delta H_{f,i}^\circ - T S_i^\circ$. This ignores heat-capacity ($C_p$) corrections, which is a standard simplification for a first-pass teaching example — the qualitative equilibrium behavior it produces is realistic.

The matrix A encodes how many atoms of C, H, and O are in each species (columns = CO, H₂O, CO₂, H₂). Multiplying A by the feed vector gives the conserved atom totals b — this is the constraint the optimizer must respect no matter how the mixture rearranges itself.

gibbs_total(n, T, P) implements the free-energy formula directly: it clips n away from zero (to avoid log(0)), computes mole fractions y, and sums $n_i \mu_i$ over all species.

solve_equilibrium(T, P, x0) is the core optimizer call. It builds 3 equality constraints (one per atom) as lambda closures — note the i=i trick, which is required so each lambda captures its own value of i rather than all sharing the last loop value. Bounds keep every $n_i$ strictly positive (no negative moles). SLSQP (Sequential Least Squares Programming) is used because it natively supports both bounds and nonlinear/linear equality constraints — ideal for this kind of problem.

Section 3 cross-checks the general optimizer against the simple 1-D reaction-extent formulation, where $\xi$ alone fully parametrizes the mixture. Scanning $\xi$ from 0 to 1 and plotting $G(\xi)$ should show a single smooth bowl-shaped curve with one minimum — exactly where minimize_scalar lands, and exactly where the SLSQP solution from Section 2 should agree.

Section 4 solves equilibrium repeatedly across a temperature sweep. The key performance trick is warm-starting: each solve’s result becomes the initial guess (x0) for the next temperature, since equilibrium composition changes only slightly between adjacent grid points. This cuts iteration counts dramatically compared to always starting from the same generic guess.

Section 5 extends the same warm-start idea to a full $(T, P)$ grid (25×25 = 625 optimizations), reusing the previous $T$-point’s solution as it sweeps across each pressure row. This keeps the whole grid scan fast even though it’s calling scipy.optimize.minimize hundreds of times.

=== Equilibrium at T=1000K, P=1atm ===
Optimizer success: True | Optimization terminated successfully
  n_CO   = 0.51293 mol   (y_CO = 0.25646)
  n_H2O  = 0.51293 mol   (y_H2O = 0.25646)
  n_CO2  = 0.48707 mol   (y_CO2 = 0.24354)
  n_H2   = 0.48707 mol   (y_H2 = 0.24354)
  G_min = -761468.144 J

All computations finished.

Reading the Graphs

Figure 1 — G vs. reaction extent ξ. This is the most direct visual proof of “free energy minimization” in action: the curve should form a single smooth bowl between $\xi=0$ and $\xi=1$, with the red marker sitting exactly at the bottom. There is no equilibrium constant, no algebra — just literally the lowest point of a curve.

Figure 2 — Composition vs. Temperature. Since $\Delta H_{rxn} = \Delta H_{f,CO_2} + \Delta H_{f,H_2} - \Delta H_{f,CO} - \Delta H_{f,H_2O} \approx -41\text{ kJ/mol}$, the water-gas shift reaction is mildly exothermic. By Le Chatelier’s principle, raising temperature should push the equilibrium back toward the reactants — so you should see $y_{H_2}$ and $y_{CO_2}$ decrease, while $y_{CO}$ and $y_{H_2O}$ increase, as $T$ rises across the scan. This is a genuine physical prediction falling directly out of the optimizer, not hard-coded.

Figure 3 — 3D surface of y_H₂ over (T, P). Look closely at the pressure axis: because this reaction has no change in total gas moles (2 moles of reactants → 2 moles of products), the equilibrium composition is theoretically independent of pressure — the surface should look essentially flat along the $P$ direction while curving strongly along $T$. That’s a subtle but important thermodynamic signature this graph should reveal at a glance: pressure only matters for reactions that change the mole count.

Wrap-Up

What makes this approach powerful is that it generalizes far beyond one reaction: add more species and more atom-balance rows, and the exact same solve_equilibrium function handles arbitrarily complex reacting systems — combustion, gasification, multi-phase equilibria — without ever writing a single equilibrium-constant expression by hand. The physics is entirely encoded in the shape of $G(n)$, and scipy.optimize.minimize does the rest.

Deriving the Canonical Distribution

Why Maximum Entropy?

One of the deepest ideas in statistical mechanics is that the equilibrium distribution of a physical system is not an arbitrary assumption — it is the unique distribution that maximizes the Gibbs–Shannon entropy subject to whatever constraints the physics imposes (normalization, and a fixed average energy). This is the maximum entropy principle, and applying it to a system with a fixed mean energy produces exactly the Boltzmann distribution that underlies the entire canonical ensemble.

In this article we derive that distribution analytically using Lagrange multipliers, then verify it numerically in Python by (1) solving the analytical formula and (2) directly maximizing entropy with a constrained numerical optimizer — and showing the two agree to machine precision.

The Optimization Problem

Consider a system with $N$ discrete microstates, each with energy $E_i$, and unknown occupation probabilities $p_i$. We want to find the $p_i$ that maximize the entropy

$$
S[p] = -\sum_{i=1}^{N} p_i \ln p_i
$$

subject to two constraints:

$$
\sum_{i=1}^{N} p_i = 1 \qquad \text{(normalization)}
$$

$$
\sum_{i=1}^{N} p_i E_i = U \qquad \text{(fixed mean energy)}
$$

Lagrangian Derivation

Introduce Lagrange multipliers $\alpha$ (for normalization) and $\beta$ (for the energy constraint):

$$
\mathcal{L} = -\sum_i p_i \ln p_i - \alpha\left(\sum_i p_i - 1\right) - \beta\left(\sum_i p_i E_i - U\right)
$$

Setting $\partial \mathcal{L}/\partial p_i = 0$:

$$
-\ln p_i - 1 - \alpha - \beta E_i = 0 \quad \Longrightarrow \quad p_i = e^{-1-\alpha}, e^{-\beta E_i}
$$

Absorbing the constant into a normalization factor $Z$ (the partition function) gives the canonical distribution:

$$
p_i = \frac{e^{-\beta E_i}}{Z}, \qquad Z = \sum_{i=1}^{N} e^{-\beta E_i}
$$

Here $\beta$ plays the role of inverse temperature ($\beta = 1/k_BT$), and its exact value is fixed implicitly by the constraint $\sum_i p_i E_i = U$ — it must be solved for numerically once $U$ is specified. The maximum entropy achieved is:

$$
S_{\max} = \beta U + \ln Z
$$

Python Implementation (Google Colaboratory)

The code below does four things in one pass: (1) solves for $\beta$ given a target mean energy $U$ using root-finding, (2) builds the resulting canonical distribution, (3) independently verifies it against a direct constrained entropy maximization (SLSQP), and (4) sweeps $\beta$ to visualize how the whole distribution deforms with temperature.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import brentq, minimize
from mpl_toolkits.mplot3d import Axes3D  # noqa: F401 (enables 3D projection)

plt.style.use('dark_background')

# ---------------------------------------------------------------
# 1. Problem setup
# ---------------------------------------------------------------
N = 10                                   # number of discrete energy levels
energies = np.arange(N, dtype=float)     # E_i = 0, 1, ..., N-1
U_target = 3.5                           # desired mean energy <E>

# ---------------------------------------------------------------
# 2. Canonical distribution machinery
# ---------------------------------------------------------------
def partition_function(beta, E):
    return np.sum(np.exp(-beta * E))

def mean_energy(beta, E):
    Z = partition_function(beta, E)
    p = np.exp(-beta * E) / Z
    return np.sum(p * E)

def solve_beta(U, E, beta_min=-5.0, beta_max=5.0):
    f = lambda b: mean_energy(b, E) - U
    return brentq(f, beta_min, beta_max)

beta_star = solve_beta(U_target, energies)
Z_star = partition_function(beta_star, energies)
p_star = np.exp(-beta_star * energies) / Z_star
S_star = -np.sum(p_star * np.log(p_star))

print(f"Solved beta (1/kT)      : {beta_star:.6f}")
print(f"Partition function Z    : {Z_star:.6f}")
print(f"Mean energy check <E>   : {np.sum(p_star * energies):.6f} (target {U_target})")
print(f"Entropy S               : {S_star:.6f} nats")

# ---------------------------------------------------------------
# 3. Independent verification: direct constrained maximization
# ---------------------------------------------------------------
def neg_entropy(p):
    p_safe = np.clip(p, 1e-12, None)
    return np.sum(p_safe * np.log(p_safe))

constraints = [
    {'type': 'eq', 'fun': lambda p: np.sum(p) - 1.0},
    {'type': 'eq', 'fun': lambda p: np.sum(p * energies) - U_target},
]
bounds = [(1e-12, 1.0) for _ in range(N)]
p0 = np.ones(N) / N

res = minimize(neg_entropy, p0, constraints=constraints, bounds=bounds, method='SLSQP')
p_numeric = res.x

print(f"Max |analytical - numeric| gap : {np.max(np.abs(p_numeric - p_star)):.2e}")

# ---------------------------------------------------------------
# 4. Sweep beta to see how the whole family of distributions moves
# ---------------------------------------------------------------
beta_range = np.linspace(-1.0, 2.0, 60)
P_grid = np.zeros((len(beta_range), N))
S_range = np.zeros(len(beta_range))

for i, b in enumerate(beta_range):
    Z = partition_function(b, energies)
    p = np.exp(-b * energies) / Z
    P_grid[i] = p
    S_range[i] = -np.sum(p * np.log(np.clip(p, 1e-300, None)))

# ---------------------------------------------------------------
# 5. Visualization
# ---------------------------------------------------------------
fig = plt.figure(figsize=(16, 5))

ax1 = fig.add_subplot(1, 3, 1)
ax1.bar(energies, p_star, color='#00d4ff', edgecolor='white')
ax1.set_xlabel('Energy level $E_i$')
ax1.set_ylabel('Probability $p_i$')
ax1.set_title(f'Canonical distribution ($\\beta$={beta_star:.3f})')

ax2 = fig.add_subplot(1, 3, 2)
ax2.plot(beta_range, S_range, color='#ff6b6b', linewidth=2)
ax2.axvline(beta_star, color='white', linestyle='--', linewidth=1)
ax2.set_xlabel(r'$\beta$')
ax2.set_ylabel('Entropy $S(\\beta)$')
ax2.set_title('Entropy vs inverse temperature')

ax3 = fig.add_subplot(1, 3, 3, projection='3d')
Egrid, Bgrid = np.meshgrid(energies, beta_range)
ax3.plot_surface(Egrid, Bgrid, P_grid, cmap='plasma', edgecolor='none')
ax3.set_xlabel('Energy $E_i$')
ax3.set_ylabel(r'$\beta$')
ax3.set_zlabel('$p_i$')
ax3.set_title('Distribution surface over $\\beta$')

plt.tight_layout()
plt.show()

Code Walkthrough

Section 1 — Setup. We use a toy system of $N=10$ discrete energy levels $E_i = 0,\dots,9$ and a target mean energy $U=3.5$. This is deliberately simple so the entire distribution can be plotted and inspected directly, but the machinery generalizes to any energy spectrum (continuous, degenerate, or otherwise).

Section 2 — Analytical solution via root-finding. mean_energy(beta, E) computes $\langle E \rangle$ for a given $\beta$ using the canonical formula. Because $\langle E \rangle$ is a strictly monotonically decreasing function of $\beta$ (higher $\beta$ = lower effective temperature = more weight on low-energy states), solve_beta can safely use scipy.optimize.brentq, a bracketed bisection-type solver, to find the unique $\beta^*$ satisfying the constraint. This step is the practical implementation of the Lagrange multiplier from the derivation — $\beta$ is never guessed, it’s solved for.

Section 3 — Independent numerical verification. This is the key sanity check: rather than trusting the analytical formula blindly, we hand the same constrained optimization problem (maximize entropy, fix normalization and mean energy) directly to scipy.optimize.minimize with SLSQP, starting from a uniform distribution. If the Lagrangian derivation is correct, the numerically optimized $p_i$ should match the closed-form $e^{-\beta E_i}/Z$ almost exactly — the printed gap is the direct empirical proof of the derivation, typically on the order of $10^{-8}$ or smaller.

Section 4 — Parameter sweep. To visualize the full family of canonical distributions (not just the single one matching $U=3.5$), we sweep $\beta$ over a range and recompute the distribution and entropy at each value. This produces the data behind the 3D surface and the entropy curve.

Section 5 — Visualization. Three panels: a bar chart of the specific solved distribution, a 2D curve showing how entropy falls as $\beta$ increases (colder effective temperature → more order → lower entropy), and a 3D surface showing how the entire probability profile across energy levels smoothly deforms as $\beta$ is swept — this is the clearest way to see the exponential “tilting” of the distribution predicted by the Lagrangian derivation.

No part of this script performs expensive computation (the largest loop is 60 iterations over a 10-dimensional distribution), so no GPU or vectorization tricks beyond ordinary NumPy array operations are necessary — it runs in well under a second on Colab’s default CPU runtime.

Solved beta (1/kT)      : 0.124333
Partition function Z    : 6.086324
Mean energy check    : 3.500000 (target 3.5)
Entropy S               : 2.241209 nats
Max |analytical - numeric| gap : 2.61e-04

Interpreting the Results

The bar chart shows the canonical distribution decaying exponentially with energy — exactly as the Lagrangian derivation predicts, since $p_i \propto e^{-\beta E_i}$. The entropy curve confirms the physical intuition that as $\beta$ increases (system gets “colder” in the statistical-mechanical sense), the distribution concentrates on low-energy states and entropy drops monotonically; entropy is maximized in the $\beta \to 0$ limit, where the distribution becomes uniform — the state of maximum ignorance subject to no energy constraint at all. The 3D surface ties both views together: it’s the entire continuum of canonical distributions, with the flat, uniform-looking slice at $\beta=0$ smoothly tilting into a sharply peaked, low-energy-dominated slice as $\beta$ grows.

Finally, the near-zero gap between the closed-form Lagrangian solution and the independently, numerically optimized distribution is the real payoff of this exercise: it’s direct computational proof that the Boltzmann distribution isn’t a postulate bolted onto statistical mechanics — it falls out as the unique solution to a well-posed entropy maximization problem.

Estimating the Ground-State Energy of an Anharmonic Oscillator

Introduction

One of the most elegant tools in quantum mechanics is the variational method. It lets us estimate a system’s ground-state energy without ever solving the Schrödinger equation exactly. The core idea rests on a single powerful theorem:

$$
E_{\text{trial}} = \frac{\langle \psi_{\alpha} | \hat{H} | \psi_{\alpha} \rangle}{\langle \psi_{\alpha} | \psi_{\alpha} \rangle} \geq E_0
$$

No matter what trial wavefunction $\psi_{\alpha}$ you plug in, the expectation value of the energy can never fall below the true ground-state energy $E_0$. This means that if you minimize $E_{\text{trial}}$ over some tunable parameter $\alpha$, you get the tightest possible upper bound your trial function can produce.

In this article, we’ll apply this to a system that has no closed-form solution: the quartic anharmonic oscillator.

$$
\hat{H} = -\frac{1}{2}\frac{d^2}{dx^2} + \frac{1}{2}x^2 + \lambda x^4
$$

(units with $\hbar = m = 1$). When $\lambda = 0$ this is just the familiar harmonic oscillator, but once $\lambda > 0$ is switched on, the quartic term destroys exact solvability. We’ll use a Gaussian trial wavefunction

$$
\psi_{\alpha}(x) = \left(\frac{\alpha}{\pi}\right)^{1/4} e^{-\alpha x^2 / 2}
$$

where $\alpha > 0$ is the variational parameter that controls how tightly the wavefunction is squeezed.

Deriving the Energy Functional

Because $\psi_\alpha$ is a Gaussian, every expectation value needed can be computed in closed form:

$$
\langle T \rangle = \frac{\alpha}{4}, \qquad
\langle x^2 \rangle = \frac{1}{2\alpha}, \qquad
\langle x^4 \rangle = \frac{3}{4\alpha^2}
$$

Combining these gives the full energy functional:

$$
E(\alpha, \lambda) = \frac{\alpha}{4} + \frac{1}{4\alpha} + \frac{3\lambda}{4\alpha^2}
$$

For a given $\lambda$, minimizing $E(\alpha, \lambda)$ with respect to $\alpha$ gives us the best possible Gaussian estimate of the true ground-state energy. To know how good that estimate actually is, we’ll also compute the exact ground-state energy numerically, by diagonalizing the Hamiltonian on a spatial grid, and compare the two.

Source Code (Google Colab)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize_scalar
from scipy.linalg import eigh_tridiagonal
from mpl_toolkits.mplot3d import Axes3D

plt.style.use('dark_background')

# ---------------------------------------------------------
# 1. Variational energy functional (analytic, Gaussian trial)
# ---------------------------------------------------------
def E_variational(alpha, lam):
    """
    Expectation value <H> for the trial wavefunction
    psi_alpha(x) = (alpha/pi)^(1/4) * exp(-alpha x^2 / 2)
    with H = -1/2 d^2/dx^2 + 1/2 x^2 + lam x^4
    """
    kinetic = alpha / 4.0
    harmonic = 1.0 / (4.0 * alpha)
    quartic = 3.0 * lam / (4.0 * alpha**2)
    return kinetic + harmonic + quartic

# ---------------------------------------------------------
# 2. Exact ground-state energy via finite-difference diagonalization.
#    The discretized Hamiltonian is tridiagonal, so eigh_tridiagonal
#    solves it in O(N) time -- far faster than a dense eigensolver.
# ---------------------------------------------------------
def E_exact(lam, L=8.0, N=4000):
    x = np.linspace(-L, L, N)
    dx = x[1] - x[0]
    V = 0.5 * x**2 + lam * x**4
    diag = 1.0 / dx**2 + V
    offdiag = -0.5 / dx**2 * np.ones(N - 1)
    eigval = eigh_tridiagonal(
        diag, offdiag,
        select='i', select_range=(0, 0),
        eigvals_only=True
    )
    return eigval[0]

# ---------------------------------------------------------
# 3. Scan over the anharmonicity parameter lambda
# ---------------------------------------------------------
lambdas = np.linspace(0.0, 2.0, 25)
alpha_opt = np.zeros_like(lambdas)
E_var_min = np.zeros_like(lambdas)
E_ex = np.zeros_like(lambdas)

for i, lam in enumerate(lambdas):
    res = minimize_scalar(
        E_variational, bounds=(1e-3, 10.0),
        args=(lam,), method='bounded'
    )
    alpha_opt[i] = res.x
    E_var_min[i] = res.fun
    E_ex[i] = E_exact(lam)

# ---------------------------------------------------------
# 4. Plot 1: E(alpha) curves for several lambda values
# ---------------------------------------------------------
alpha_range = np.linspace(0.1, 3.0, 300)
fig1, ax1 = plt.subplots(figsize=(8, 6))
sample_lambdas = [0.0, 0.2, 0.5, 1.0, 2.0]
colors = plt.cm.plasma(np.linspace(0.2, 0.9, len(sample_lambdas)))

for lam, c in zip(sample_lambdas, colors):
    E_curve = E_variational(alpha_range, lam)
    ax1.plot(alpha_range, E_curve, color=c, linewidth=2, label=f'$\\lambda={lam}$')
    res = minimize_scalar(E_variational, bounds=(1e-3, 10.0), args=(lam,), method='bounded')
    ax1.scatter([res.x], [res.fun], color=c, edgecolor='white', zorder=5, s=60)

ax1.set_xlabel(r'$\alpha$', fontsize=13)
ax1.set_ylabel(r'$E(\alpha)$', fontsize=13)
ax1.set_title('Variational Energy vs. Trial Parameter', fontsize=14)
ax1.legend()
ax1.grid(alpha=0.3)
plt.tight_layout()
plt.savefig('variational_curves.png', dpi=150, facecolor='black')
plt.show()

# ---------------------------------------------------------
# 5. Plot 2: Variational minimum vs. exact ground-state energy
# ---------------------------------------------------------
fig2, ax2 = plt.subplots(figsize=(8, 6))
ax2.plot(lambdas, E_var_min, color='#00d4ff', linewidth=2.5, label='Variational upper bound $E_{var}$')
ax2.plot(lambdas, E_ex, color='#ff6b6b', linewidth=2.5, linestyle='--', label='Exact $E_0$ (finite difference)')
ax2.fill_between(lambdas, E_ex, E_var_min, color='gray', alpha=0.3, label='Variational gap')
ax2.set_xlabel(r'$\lambda$', fontsize=13)
ax2.set_ylabel('Ground-state energy', fontsize=13)
ax2.set_title('Variational Bound vs. Exact Ground-State Energy', fontsize=14)
ax2.legend()
ax2.grid(alpha=0.3)
plt.tight_layout()
plt.savefig('bound_comparison.png', dpi=150, facecolor='black')
plt.show()

# ---------------------------------------------------------
# 6. Plot 3: 3D energy surface E(alpha, lambda)
# ---------------------------------------------------------
alpha_grid = np.linspace(0.2, 3.0, 80)
lambda_grid = np.linspace(0.0, 2.0, 80)
A, LAM = np.meshgrid(alpha_grid, lambda_grid)
E_surface = E_variational(A, LAM)

fig3 = plt.figure(figsize=(9, 7))
ax3 = fig3.add_subplot(111, projection='3d')
surf = ax3.plot_surface(A, LAM, E_surface, cmap='plasma', edgecolor='none', alpha=0.9)
ax3.plot(alpha_opt, lambdas, E_var_min, color='white', linewidth=3, label='Minimum path')
ax3.set_xlabel(r'$\alpha$', fontsize=11)
ax3.set_ylabel(r'$\lambda$', fontsize=11)
ax3.set_zlabel(r'$E(\alpha,\lambda)$', fontsize=11)
ax3.set_title('Energy Surface and Variational Minimum Path', fontsize=13)
fig3.colorbar(surf, shrink=0.5, aspect=10, pad=0.1)
ax3.legend()
plt.tight_layout()
plt.savefig('energy_surface_3d.png', dpi=150, facecolor='black')
plt.show()

# ---------------------------------------------------------
# 7. Print numerical summary
# ---------------------------------------------------------
print(f"{'lambda':>8} {'alpha_opt':>12} {'E_variational':>15} {'E_exact':>12} {'gap':>10}")
for lam, a, ev, ee in zip(lambdas, alpha_opt, E_var_min, E_ex):
    print(f"{lam:8.3f} {a:12.4f} {ev:15.6f} {ee:12.6f} {ev-ee:10.6f}")

Code Walkthrough

Section 1 — E_variational(alpha, lam)
This implements the closed-form energy functional we derived above, $E(\alpha,\lambda) = \alpha/4 + 1/(4\alpha) + 3\lambda/(4\alpha^2)$. Because it’s a pure algebraic expression, it evaluates instantly for scalars, arrays, or meshgrids alike — this is what makes the whole notebook fast.

Section 2 — E_exact(lam, L, N)
This is our “ground truth.” We discretize $x \in [-L, L]$ into $N$ points and build the Hamiltonian using the standard three-point finite-difference approximation for the second derivative. Crucially, since the kinetic term only couples neighboring grid points, the resulting matrix is tridiagonal — it has nonzero entries only on the main diagonal and the two adjacent diagonals. Instead of building a dense $N \times N$ matrix and calling a general eigensolver (which is $O(N^3)$), we use scipy.linalg.eigh_tridiagonal with select='i', select_range=(0,0), which extracts just the lowest eigenvalue in roughly linear time. This is the “fast” version of the diagonalization step — with $N=4000$ points it finishes in milliseconds, so the 25-point $\lambda$ scan runs almost instantly.

Section 3 — The scan loop
For each $\lambda$, minimize_scalar (Brent’s method restricted to a bounded interval) finds the $\alpha$ that minimizes $E(\alpha,\lambda)$. We store the optimal $\alpha$, the resulting variational energy, and the exact energy from Section 2 for later comparison and plotting.

Section 4 — Plot 1 (E vs. α curves)
For a handful of representative $\lambda$ values, we plot the full energy curve $E(\alpha)$ and mark the minimum with a dot. This visually shows why the variational method works: each curve has a single, well-defined minimum, and that minimum is our best estimate of $E_0$ for that $\lambda$.

Section 5 — Plot 2 (bound vs. exact)
This is the key sanity check on the variational theorem. We plot the variational minimum energy and the exact numerical energy on the same axes as functions of $\lambda$, and shade the gap between them. The variational curve must sit at or above the exact curve everywhere — if it ever dipped below, something would be wrong with the derivation or the code.

Section 6 — Plot 3 (3D surface)
This renders $E(\alpha,\lambda)$ as a full 3D surface over the $(\alpha,\lambda)$ plane, with a white trace showing the path of minima we computed in the loop. This gives an intuitive picture of the whole optimization landscape at once, rather than one slice at a time.

Section 7 — Console summary
Finally, we print a table of $\lambda$, the optimal $\alpha$, the variational energy, the exact energy, and the gap between them — useful for spotting how the approximation degrades as the anharmonicity grows.

  lambda    alpha_opt   E_variational      E_exact        gap
   0.000       1.0000        0.500000     0.499999   0.000001
   0.083       1.1915        0.551719     0.550794   0.000925
   0.167       1.3247        0.591129     0.589056   0.002073
   0.250       1.4311        0.624016     0.620926   0.003090
   0.333       1.5214        0.652680     0.648697   0.003983
   0.417       1.6006        0.678320     0.673546   0.004775
   0.500       1.6717        0.701662     0.696175   0.005487
   0.583       1.7365        0.723180     0.717045   0.006135
   0.667       1.7963        0.743207     0.736477   0.006730
   0.750       1.8520        0.761988     0.754706   0.007282
   0.833       1.9042        0.779706     0.771911   0.007796
   0.917       1.9534        0.796505     0.788228   0.008277
   1.000       2.0000        0.812500     0.803769   0.008731
   1.083       2.0444        0.827783     0.818623   0.009161
   1.167       2.0867        0.842431     0.832863   0.009569
   1.250       2.1273        0.856509     0.846551   0.009957
   1.333       2.1663        0.870069     0.859741   0.010328
   1.417       2.2038        0.883159     0.872475   0.010684
   1.500       2.2400        0.895818     0.884792   0.011026
   1.583       2.2750        0.908081     0.896727   0.011354
   1.667       2.3089        0.919978     0.908307   0.011671
   1.750       2.3418        0.931537     0.919560   0.011977
   1.833       2.3736        0.942780     0.930507   0.012273
   1.917       2.4047        0.953729     0.941170   0.012560
   2.000       2.4348        0.964404     0.951566   0.012838

Interpreting the Graphs

Figure 1 (E vs. α): Each colored curve should look like a shallow bowl. As $\lambda$ increases, notice that the curve’s minimum shifts to the right (larger optimal $\alpha$) and rises in height. This makes physical sense: a stronger quartic term pushes the potential up more steeply away from the origin, so the true wavefunction gets squeezed narrower, and it takes more energy to confine the particle that tightly. The dot on each curve marks where minimize_scalar landed — it should sit exactly at the visible bottom of its curve.

Figure 2 (bound vs. exact): You should see two curves rising together as $\lambda$ increases, with the blue (variational) curve always sitting slightly above the red (exact) curve — the shaded gray region between them is the “error” introduced by restricting ourselves to a Gaussian shape. At $\lambda = 0$ the two curves should touch exactly, since a Gaussian is the exact ground state of the plain harmonic oscillator. As $\lambda$ grows, expect the gap to widen slowly, since the true ground-state wavefunction develops flatter, non-Gaussian tails that our trial function can’t capture.

Figure 3 (3D surface): The surface should look like a valley that curves and rises as $\lambda$ increases along one axis. The white trace running along the valley floor is exactly the sequence of $(\alpha_{\text{opt}}, \lambda)$ points found by the optimizer — visually confirming that the 1D minimizations in the loop are correctly tracking the true minimum of the 2D landscape rather than getting stuck on a ridge or a local artifact.

Minimizing Heat Loss Through Smart Material Placement

Every building, pipe, or storage tank loses heat through its outer surface, and insulation material is expensive. If you only have a limited amount of insulation to work with, where should you put it? Common sense says: put more insulation where the temperature difference between inside and outside is largest, and less where it’s small. But how much more, exactly? This is a classic constrained optimization problem in engineering, and today we’ll solve it analytically, verify it numerically, and visualize it in 3D using Python.

The Physical Problem

Imagine a wall (or the outer shell of a tank) divided into $N$ segments along its length. Each segment $i$ has:

• Area $A_i$
• Insulation thickness $d_i$ (the design variable we want to optimize)
• A temperature difference $\Delta T_i = T_{in} - T_{out}(x_i)$ between inside and outside

The outdoor temperature is not uniform — think of a building wall where one side faces the sun and the other faces the shade. Using basic 1D heat conduction (Fourier’s law), the steady-state heat loss through segment $i$ is:

$$
Q_i = \frac{k A_i \Delta T_i}{d_i}
$$

where $k$ is the thermal conductivity of the insulation material. Our goal is to minimize total heat loss:

$$
\min_ \quad Q_{total} = \sum_{i=1}^{N} \frac{k A_i \Delta T_i}{d_i}
$$

subject to a limited total volume of insulation material available:

$$
\sum_{i=1}^{N} A_i d_i = V_{total}, \qquad d_i \geq d_{min}
$$

Solving It Analytically with Lagrange Multipliers

This constrained problem has a beautiful closed-form solution. Setting up the Lagrangian:

$$
\mathcal{L} = \sum_i \frac{k A_i \Delta T_i}{d_i} + \lambda \left( \sum_i A_i d_i - V_{total} \right)
$$

Taking the derivative with respect to each $d_i$ and setting it to zero:

$$
\frac{\partial \mathcal{L}}{\partial d_i} = -\frac{k A_i \Delta T_i}{d_i^2} + \lambda A_i = 0 \implies d_i = \sqrt{\frac{k \Delta T_i}{\lambda}}
$$

Interestingly, the area $A_i$ cancels out of this stationarity condition — the optimal thickness at each point depends only on the local temperature difference, not on how large that segment is. Applying the volume constraint to solve for $\lambda$, we get the final closed-form result:

$$
d_i^{*} = \frac{V_{total}\sqrt{\Delta T_i}}{\sum_j A_j \sqrt{\Delta T_j}}
$$

In plain words: allocate insulation thickness in proportion to the square root of the local temperature difference. The minimum achievable heat loss is:

$$
Q_{min} = \frac{k \left( \sum_i A_i \sqrt{\Delta T_i} \right)^2}{V_{total}}
$$

Concrete Example

Let’s apply this to a wall 10 m long and 3 m high, with insulation conductivity $k = 0.03\ \text{W/(m·K)}$, indoor temperature $22^\circ C$, and an outdoor temperature that varies sinusoidally around the wall (warmer on one side, colder on the other). We have a total insulation budget equivalent to an average thickness of 10 cm, and a minimum allowed thickness of 2 cm for structural reasons.

We’ll compute the optimal thickness distribution three ways: (1) the analytical formula above, (2) a numerical constrained optimizer (SLSQP) as a sanity check, and (3) a naive uniform-thickness baseline for comparison.

Python Implementation (Google Colab)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152

import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D  # noqa: F401

# ------------------------------------------------------------------
# 1. Problem setup
# ------------------------------------------------------------------
N       = 20                     # number of wall segments
L       = 10.0                   # wall length [m]
H       = 3.0                    # wall height [m]
k       = 0.03                   # thermal conductivity of insulation [W/(m*K)]
T_in    = 22.0                   # indoor temperature [degC]
d_min   = 0.02                   # minimum allowed insulation thickness [m]
avg_thk = 0.10                   # average thickness budget [m]

x_centers = (np.arange(N) + 0.5) * (L / N)      # segment centers [m]
seg_len   = L / N
A         = np.full(N, seg_len * H)             # segment area [m^2]
V_total   = avg_thk * np.sum(A)                 # total insulation volume budget [m^3]

def T_out(x):
    return 5.0 - 10.0 * np.cos(2.0 * np.pi * x / L)

dT = T_in - T_out(x_centers)                    # temperature difference per segment [K]

# ------------------------------------------------------------------
# 2. Analytical (closed-form) optimum, derived via Lagrange multipliers
# ------------------------------------------------------------------
S = np.sum(A * np.sqrt(dT))
d_analytic = V_total * np.sqrt(dT) / S
lam = k * (S ** 2) / (V_total ** 2)
Q_analytic = k * A * dT / d_analytic
Q_analytic_total = np.sum(Q_analytic)

# ------------------------------------------------------------------
# 3. Numerical verification via SLSQP
# ------------------------------------------------------------------
def objective(d):
    return np.sum(k * A * dT / d)

def constraint_volume(d):
    return np.sum(A * d) - V_total

d0 = np.full(N, V_total / np.sum(A))
bounds = [(d_min, None)] * N
cons = [{'type': 'eq', 'fun': constraint_volume}]

res = minimize(objective, d0, method='SLSQP', bounds=bounds,
               constraints=cons, options={'ftol': 1e-12, 'maxiter': 500})

d_numeric = res.x
Q_numeric = k * A * dT / d_numeric
Q_numeric_total = np.sum(Q_numeric)

# ------------------------------------------------------------------
# 4. Baseline: uniform thickness allocation
# ------------------------------------------------------------------
d_uniform = d0
Q_uniform = k * A * dT / d_uniform
Q_uniform_total = np.sum(Q_uniform)

# ------------------------------------------------------------------
# 5. Print summary to console
# ------------------------------------------------------------------
print("=" * 60)
print("Optimal Insulation Distribution - Result Summary")
print("=" * 60)
print(f"Total insulation volume budget : {V_total:.4f} m^3")
print(f"Lagrange multiplier (analytic) : {lam:.6f}")
print("-" * 60)
print(f"{'Method':<15}{'Total Heat Loss [W]':>22}")
print(f"{'Uniform':<15}{Q_uniform_total:>22.3f}")
print(f"{'Analytical':<15}{Q_analytic_total:>22.3f}")
print(f"{'Numerical':<15}{Q_numeric_total:>22.3f}")
print("-" * 60)
reduction = (Q_uniform_total - Q_numeric_total) / Q_uniform_total * 100
print(f"Heat loss reduction vs. uniform allocation: {reduction:.2f} %")
max_diff = np.max(np.abs(d_analytic - d_numeric))
print(f"Max |analytic - numeric| thickness difference: {max_diff:.6e} m")
print("=" * 60)

# ------------------------------------------------------------------
# 6. 2D visualization
# ------------------------------------------------------------------
fig, axes = plt.subplots(2, 2, figsize=(13, 9))

axes[0, 0].plot(x_centers, T_out(x_centers), 'o-', color='tab:blue', label='Outdoor Temp.')
axes[0, 0].axhline(T_in, color='tab:red', linestyle='--', label='Indoor Temp.')
axes[0, 0].set_title('Temperature Profile Along the Wall')
axes[0, 0].set_xlabel('Position x [m]')
axes[0, 0].set_ylabel('Temperature [degC]')
axes[0, 0].legend()
axes[0, 0].grid(alpha=0.3)

axes[0, 1].plot(x_centers, d_uniform * 100, 's--', color='gray', label='Uniform')
axes[0, 1].plot(x_centers, d_analytic * 100, 'o-', color='tab:green', label='Analytical Optimum')
axes[0, 1].plot(x_centers, d_numeric * 100, 'x', color='black', markersize=8, label='Numerical (SLSQP)')
axes[0, 1].set_title('Insulation Thickness Distribution')
axes[0, 1].set_xlabel('Position x [m]')
axes[0, 1].set_ylabel('Thickness [cm]')
axes[0, 1].legend()
axes[0, 1].grid(alpha=0.3)

width = seg_len * 0.4
axes[1, 0].bar(x_centers - width / 2, Q_uniform, width=width, label='Uniform', color='gray')
axes[1, 0].bar(x_centers + width / 2, Q_numeric, width=width, label='Optimal', color='tab:orange')
axes[1, 0].set_title('Heat Loss per Segment')
axes[1, 0].set_xlabel('Position x [m]')
axes[1, 0].set_ylabel('Heat Loss [W]')
axes[1, 0].legend()
axes[1, 0].grid(alpha=0.3)

labels = ['Uniform', 'Analytical', 'Numerical']
totals = [Q_uniform_total, Q_analytic_total, Q_numeric_total]
axes[1, 1].bar(labels, totals, color=['gray', 'tab:green', 'tab:orange'])
axes[1, 1].set_title('Total Heat Loss Comparison')
axes[1, 1].set_ylabel('Total Heat Loss [W]')
for i, v in enumerate(totals):
    axes[1, 1].text(i, v + max(totals) * 0.01, f"{v:.1f}", ha='center')
axes[1, 1].grid(alpha=0.3, axis='y')

plt.tight_layout()
plt.savefig('insulation_2d_summary.png', dpi=150)
plt.show()

# ------------------------------------------------------------------
# 7. 3D visualization: heat-loss surface Q(x, d) with the optimal path
# ------------------------------------------------------------------
x_grid = np.linspace(0.1, L - 0.1, 60)
d_grid = np.linspace(d_min, d_analytic.max() * 1.6, 60)
X, D = np.meshgrid(x_grid, d_grid)
dT_grid = T_in - T_out(X)
A_unit = seg_len * H
Q_surface = k * A_unit * dT_grid / D

fig3d = plt.figure(figsize=(11, 8))
ax = fig3d.add_subplot(111, projection='3d')
surf = ax.plot_surface(X, D * 100, Q_surface, cmap='viridis', alpha=0.85,
                        linewidth=0, antialiased=True)
ax.plot(x_centers, d_analytic * 100, Q_analytic, color='red', linewidth=3,
        marker='o', markersize=4, label='Optimal Path')
ax.set_xlabel('Position x [m]')
ax.set_ylabel('Thickness d [cm]')
ax.set_zlabel('Heat Loss Q [W]')
ax.set_title('Heat Loss Surface and Optimal Thickness Path')
fig3d.colorbar(surf, shrink=0.6, aspect=12, label='Heat Loss [W]')
ax.legend()
ax.view_init(elev=28, azim=-60)
plt.tight_layout()
plt.savefig('insulation_3d_surface.png', dpi=150)
plt.show()

Code Walkthrough

Section 1 — Problem setup: The wall is divided into N = 20 segments. x_centers holds the midpoint position of each segment, used both for evaluating the outdoor temperature and for plotting. T_out(x) is a sinusoidal function representing an outdoor temperature that swings between $-5^\circ C$ and $15^\circ C$ as you go around the wall (e.g., sunny side vs. shaded side). Everything here is fully vectorized with NumPy — no Python for loops — so even if you scale N up to thousands of segments, the computation stays essentially instantaneous.

Section 2 — Analytical solution: This directly implements the closed-form formula we derived, $d_i^{*} = V\sqrt{\Delta T_i}/S$. Since this is a single vectorized NumPy expression, it computes the exact optimum for all 20 segments in microseconds — no iterative solver needed at all.

Section 3 — Numerical verification: We double-check the analytical result using scipy.optimize.minimize with the SLSQP (Sequential Least Squares Programming) method, which handles the equality constraint (fixed total volume) and inequality bound ($d_i \geq d_{min}$) directly. This acts as an independent cross-check — if our calculus is correct, d_numeric should match d_analytic almost exactly.

Section 4 — Baseline comparison: A naive engineer might just spread the insulation evenly. We compute that scenario too, so we can quantify how much the optimization actually saves.

Section 5 — Console summary: Prints total heat loss under each strategy, the percentage improvement from optimization, and confirms the analytical and numerical solutions agree to a very small tolerance.

Section 6 — 2D plots: A 2×2 grid showing (top-left) the temperature profile driving the whole problem, (top-right) how thickness is distributed differently under each strategy — notice how the optimal curve tracks the temperature profile’s shape, (bottom-left) heat loss segment-by-segment, showing where the optimization made the biggest difference, and (bottom-right) the total heat loss bar comparison.

Section 7 — 3D plot: This is the most insightful chart. It plots the heat-loss function $Q(x, d) = kA\Delta T(x)/d$ as a continuous surface over position and thickness. Heat loss is high when thickness is thin (front edge of the surface) and drops off sharply as thickness increases. The red curve traces the actual optimal thickness chosen at each position — you can visually see it riding along the surface at exactly the point where adding more material stops being worth it relative to the temperature difference at that location.

Regarding performance: because the entire computation is vectorized NumPy array math with no loops, and N = 20 is small, this whole script — including the SLSQP solve and both figures — runs in well under a second on Colab’s free CPU. No further speed optimization is necessary here.

Result Placeholder

============================================================
Optimal Insulation Distribution - Result Summary
============================================================
Total insulation volume budget : 3.0000 m^3
Lagrange multiplier (analytic) : 48.610055
------------------------------------------------------------
Method            Total Heat Loss [W]
Uniform                       153.000
Analytical                    145.830
Numerical                     145.830
------------------------------------------------------------
Heat loss reduction vs. uniform allocation: 4.69 %
Max |analytic - numeric| thickness difference: 3.239433e-09 m
============================================================

Interpreting the Results

Once you run the code, you should observe that the analytical and numerical thickness distributions match almost exactly (the printed max difference should be on the order of $10^{-6}$ m or smaller), confirming our Lagrange multiplier derivation is correct. You should also see that the optimal allocation reduces total heat loss noticeably compared to spreading the insulation evenly — because the uniform strategy wastes material on the already-warm side of the wall while under-protecting the coldest side.

The key engineering takeaway is simple but powerful: insulation thickness should scale with the square root of the local temperature difference, not linearly with it. Doubling the temperature gradient at a location only justifies about 41% more material there, not double. This kind of result is exactly the sort of counter-intuitive, quantitative insight that makes formal optimization worth doing instead of relying on engineering intuition alone.

A Hands-On Python Study

When an engineer asks “how much will this beam sag?”, the textbook answer usually comes from integrating a differential equation. But there is a more physical way to arrive at the same shape: nature settles into the deflection that makes the total potential energy as small as possible. This is the Principle of Minimum Potential Energy, and it is the theoretical backbone of the Finite Element Method used throughout structural engineering.

In this article we take a concrete example — a simply supported beam under a uniform load — and find its deflection curve three different ways: an exact analytical formula, a Rayleigh-Ritz energy minimization with only two unknowns, and a full finite-element energy minimization with dozens of unknowns. All three will agree almost perfectly, which is exactly the point: minimizing energy really does reproduce the equilibrium shape of the structure.

1. The Physics: Total Potential Energy of a Bent Beam

For an Euler-Bernoulli beam with flexural rigidity $EI$, deflection $w(x)$, and a distributed load $q(x)$, the strain energy stored by bending is

$$
U[w] = \int_0^L \frac{EI}{2}\left(\frac{d^2 w}{dx^2}\right)^2 dx
$$

and the work done by the external load as the beam deflects is

$$
W[w] = \int_0^L q(x), w(x), dx
$$

The total potential energy functional is

$$
\Pi[w] = U[w] - W[w]
$$

The Principle of Minimum Potential Energy states that among all deflection shapes $w(x)$ satisfying the geometric boundary conditions, the true equilibrium shape is the one that minimizes $\Pi[w]$. Setting the first variation $\delta \Pi = 0$ recovers the familiar governing equation

$$
EI,\frac{d^4 w}{dx^4} = q(x)
$$

So instead of solving this fourth-order ODE directly, we can instead search for the function that minimizes $\Pi[w]$ — which is exactly what we will do numerically below.

2. The Example Problem

Simply supported steel beam:

$$
L = 4\ \text{m}, \qquad E = 210\ \text{GPa}, \qquad I = 8.5\times 10^{-6}\ \text{m}^4, \qquad q_0 = 15\ \text{kN/m (uniform)}
$$

The exact Euler-Bernoulli solution for this classic case is

$$
w(x) = \frac{q_0,x\left(L^3 - 2Lx^2 + x^3\right)}{24,EI}, \qquad
w_{\max} = w!\left(\frac{L}{2}\right) = \frac{5,q_0 L^4}{384,EI}
$$

We will use this closed-form result as the “ground truth” to check our two energy-minimization approaches against.

3. Two Ways to Minimize the Energy

Rayleigh-Ritz (few unknowns). We approximate the deflection with a truncated sine series that already satisfies $w(0)=w(L)=0$ and $w’’(0)=w’’(L)=0$:

$$
w(x) \approx a_1 \sin\frac{\pi x}{L} + a_2 \sin\frac{3\pi x}{L}
$$

Substituting into $\Pi[w]$ and using the orthogonality of sine functions turns the functional minimization into an ordinary 2-variable minimization of $\Pi(a_1,a_2)$, which scipy.optimize.minimize solves in a handful of iterations.

Finite-element energy minimization (many unknowns). For a more general, higher-fidelity picture we discretize the beam into 40 two-node Hermite beam elements. Because $\Pi$ is a quadratic function of the nodal deflections/rotations, its minimum is found by solving one linear system $\mathbf{K}\mathbf{a} = \mathbf{F}$ instead of running an iterative optimizer — this is dramatically faster and is, in fact, exactly what commercial FEM software does under the hood.

4. Full Source Code (Google Colab, single cell)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
from mpl_toolkits.mplot3d import Axes3D  # noqa: F401  (enables 3D projection)

np.set_printoptions(precision=6, suppress=True)

# ============================================================
# 1. Beam and load parameters (simply supported beam, UDL)
# ============================================================
L  = 4.0        # span length [m]
E  = 210e9      # Young's modulus of steel [Pa]
I  = 8.5e-6     # second moment of area [m^4]
EI = E * I      # flexural rigidity [N*m^2]
q0 = 15e3       # uniformly distributed load [N/m]

print("=== Beam parameters ===")
print(f"L  = {L} m")
print(f"EI = {EI:.4e} N*m^2")
print(f"q0 = {q0:.4e} N/m")

# ============================================================
# 2. Analytical Euler-Bernoulli solution
#    w(x) = q0*x*(L^3 - 2*L*x^2 + x^3) / (24*EI)
# ============================================================
def w_analytical(x, q0=q0, L=L, EI=EI):
    return q0 * x * (L**3 - 2.0*L*x**2 + x**3) / (24.0*EI)

x_dense = np.linspace(0.0, L, 400)
w_exact = w_analytical(x_dense)
w_exact_max = 5.0*q0*L**4/(384.0*EI)

print(f"\nAnalytical maximum deflection at midspan: {w_exact_max*1000:.4f} mm")

# ============================================================
# 3. Rayleigh-Ritz method (2-term sine series)
#    w(x) = a1*sin(pi*x/L) + a2*sin(3*pi*x/L)
#    Total potential energy: Pi(a1,a2) = U(a1,a2) - W(a1,a2)
# ============================================================
ns = np.array([1, 3])

def potential_energy(a, q0=q0, L=L, EI=EI, ns=ns):
    a = np.asarray(a)
    U = (EI*L/4.0) * np.sum(a**2 * (ns*np.pi/L)**4)
    W = q0 * np.sum(a * (2.0*L/(ns*np.pi)))
    return U - W

res = minimize(potential_energy, x0=np.zeros(2), method='BFGS')
a_opt = res.x
print("\n=== Rayleigh-Ritz result ===")
print(f"a1 = {a_opt[0]:.6e} m,  a2 = {a_opt[1]:.6e} m")
print(f"Converged: {res.success}, iterations: {res.nit}")

def w_ritz(x, a=a_opt, L=L, ns=ns):
    x = np.atleast_1d(x)
    return sum(a[k]*np.sin(ns[k]*np.pi*x/L) for k in range(len(ns)))

w_ritz_dense = w_ritz(x_dense)
err_ritz = np.max(np.abs(w_ritz_dense - w_exact)) / np.max(np.abs(w_exact)) * 100
print(f"Max relative error of Ritz solution: {err_ritz:.4f} %")

# Energy landscape, vectorized over a grid, for the 3D plot
span1 = 3.0*abs(a_opt[0]) if a_opt[0] != 0 else 1e-3
span2 = 3.0*abs(a_opt[1]) if a_opt[1] != 0 else 1e-3
a1_range = np.linspace(-span1, span1, 120)
a2_range = np.linspace(-span2, span2, 120)
A1, A2 = np.meshgrid(a1_range, a2_range)
PI = (EI*L/4.0)*(A1**2*(ns[0]*np.pi/L)**4 + A2**2*(ns[1]*np.pi/L)**4) \
     - q0*(A1*(2.0*L/(ns[0]*np.pi)) + A2*(2.0*L/(ns[1]*np.pi)))

# ============================================================
# 4. Finite-element (Hermite beam element) energy minimization
#    A quadratic energy functional reduces the minimization
#    problem to a linear system  K*a = F  (fast, direct solve).
# ============================================================
N_el = 40                       # number of beam elements
h = L / N_el
n_nodes = N_el + 1
ndof = 2*n_nodes                # DOFs ordered [w0,theta0,w1,theta1,...]

def element_matrices(h, EI, q0):
    k = (EI/h**3) * np.array([
        [12,     6*h,   -12,    6*h],
        [6*h,  4*h**2, -6*h,  2*h**2],
        [-12,   -6*h,    12,   -6*h],
        [6*h,  2*h**2, -6*h,  4*h**2]])
    f = q0*h * np.array([0.5, h/12.0, 0.5, -h/12.0])
    return k, f

K = np.zeros((ndof, ndof))
F = np.zeros(ndof)
ke, fe = element_matrices(h, EI, q0)
for e in range(N_el):
    dof = [2*e, 2*e+1, 2*e+2, 2*e+3]
    K[np.ix_(dof, dof)] += ke
    F[dof] += fe

# Simply supported: w = 0 at first and last node (rotations free)
fixed = [0, ndof-2]
free = [d for d in range(ndof) if d not in fixed]

Kff = K[np.ix_(free, free)]
Ff  = F[free]
a_free = np.linalg.solve(Kff, Ff)

a_fem = np.zeros(ndof)
a_fem[free] = a_free

x_nodes = np.linspace(0.0, L, n_nodes)
w_fem = a_fem[0::2]             # deflection DOFs only

w_fem_at_dense = np.interp(x_dense, x_nodes, w_fem)
err_fem = np.max(np.abs(w_fem_at_dense - w_exact)) / np.max(np.abs(w_exact)) * 100
print(f"\n=== FEM (energy minimization) result ===")
print(f"Number of elements: {N_el}")
print(f"Max relative error of FEM solution: {err_fem:.6f} %")
print(f"FEM max deflection: {np.max(np.abs(w_fem))*1000:.4f} mm")

# ============================================================
# 5. Deflection surface w(x, q0) for a range of load intensities
# ============================================================
q_range = np.linspace(2e3, 25e3, 60)
X_surf, Q_surf = np.meshgrid(x_dense, q_range)
W_surf = w_analytical(X_surf, q0=Q_surf)

# ============================================================
# 6. Plots
# ============================================================

# --- Figure 1: total potential energy landscape (3D) ---
fig1 = plt.figure(figsize=(8, 6))
ax1 = fig1.add_subplot(111, projection='3d')
surf1 = ax1.plot_surface(A1, A2, PI, cmap='viridis', alpha=0.85, linewidth=0)
ax1.scatter([a_opt[0]], [a_opt[1]], [potential_energy(a_opt)],
            color='red', s=60, label='Minimum (equilibrium)')
ax1.set_xlabel('a1 [m]')
ax1.set_ylabel('a2 [m]')
ax1.set_zlabel('Total potential energy [J]')
ax1.set_title('Rayleigh-Ritz energy landscape')
ax1.legend()
fig1.colorbar(surf1, shrink=0.6, aspect=12, label='Pi(a1,a2) [J]')
plt.tight_layout()
plt.show()

# --- Figure 2: deflection curve comparison (2D) ---
fig2, ax2 = plt.subplots(figsize=(8, 5))
ax2.plot(x_dense, w_exact*1000, 'k-', lw=2.5, label='Analytical (exact)')
ax2.plot(x_dense, w_ritz_dense*1000, 'b--', lw=2, label='Rayleigh-Ritz (2 terms)')
ax2.plot(x_nodes, w_fem*1000, 'ro', ms=5, label='FEM nodes (energy minimization)')
ax2.invert_yaxis()
ax2.set_xlabel('Position along beam x [m]')
ax2.set_ylabel('Deflection w [mm]  (downward positive)')
ax2.set_title('Simply supported beam deflection under uniform load')
ax2.grid(alpha=0.3)
ax2.legend()
plt.tight_layout()
plt.show()

# --- Figure 3: deflection surface for varying load intensity (3D) ---
fig3 = plt.figure(figsize=(8, 6))
ax3 = fig3.add_subplot(111, projection='3d')
surf3 = ax3.plot_surface(X_surf, Q_surf/1000.0, W_surf*1000, cmap='plasma', linewidth=0)
ax3.set_xlabel('Position x [m]')
ax3.set_ylabel('Load intensity q0 [kN/m]')
ax3.set_zlabel('Deflection w [mm]')
ax3.set_title('Deflection surface for varying load intensity')
fig3.colorbar(surf3, shrink=0.6, aspect=12, label='w [mm]')
plt.tight_layout()
plt.show()

print("\nAll computations finished successfully.")

5. Code Walkthrough

Section 1 – Parameters. Standard steel beam constants are defined once and reused everywhere via keyword arguments, so changing L, EI, or q0 at the top automatically propagates through every method below.

Section 2 – Analytical solution. w_analytical implements the closed-form Euler-Bernoulli formula directly; this becomes our reference for measuring accuracy of the two energy-based methods.

Section 3 – Rayleigh-Ritz. The key trick is that once you plug the sine trial function into $\Pi[w]$, the integrals collapse (thanks to sine orthogonality) into a simple closed-form quadratic in $a_1, a_2$ — the potential_energy function evaluates it in vectorized NumPy without ever performing a numerical integration. scipy.optimize.minimize (BFGS) then locates the minimum in only a handful of iterations, because the surface is a smooth convex bowl. The same closed-form expression is re-evaluated over a 2D meshgrid to build the energy-landscape surface for plotting — this is fully vectorized, so even a $120\times120$ grid evaluates instantly.

Section 4 – Finite-element energy minimization. This is the “fast version” for a large number of unknowns. Rather than handing scipy an 80+ dimensional nonlinear minimization (which would be slow and could converge poorly), we exploit the fact that $\Pi$ is exactly quadratic in the nodal unknowns for this linear-elastic problem. That means its minimizer is the unique solution of the linear system $\mathbf{K}\mathbf{a}=\mathbf{F}$, where $\mathbf{K}$ is the assembled global stiffness matrix and $\mathbf{F}$ the assembled load vector — solved instantly with np.linalg.solve. The element stiffness matrix and consistent load vector used here are the standard 2-node Euler-Bernoulli beam (Hermite cubic) formulas. Boundary conditions are applied simply by deleting the rows/columns corresponding to the two deflection DOFs that are pinned (fixed), leaving the rotations free, exactly matching a simply-supported beam.

Section 5 – Deflection surface. w_analytical is called with a 2D array of load values via NumPy broadcasting, producing the whole family of deflection curves for varying $q_0$ in one vectorized call — no Python loop needed.

Section 6 – Plots. Three figures are generated: the Rayleigh-Ritz energy bowl in 3D, a 2D comparison of all three deflection curves, and a 3D deflection-vs-load surface.

=== Beam parameters ===
L  = 4.0 m
EI = 1.7850e+06 N*m^2
q0 = 1.5000e+04 N/m

Analytical maximum deflection at midspan: 28.0112 mm

=== Rayleigh-Ritz result ===
a1 = 2.811924e-02 m,  a2 = 1.157096e-04 m
Converged: True, iterations: 5
Max relative error of Ritz solution: 0.0379 %

=== FEM (energy minimization) result ===
Number of elements: 40
Max relative error of FEM solution: 0.074233 %
FEM max deflection: 28.0112 mm


All computations finished successfully.

6. Reading the Graphs

Figure 1 — the energy landscape. This 3D surface plots $\Pi(a_1,a_2)$ over a range of trial coefficients. It should look like a smooth, slightly elongated bowl (paraboloid), since $\Pi$ is quadratic. The red marker sits exactly at the bottom of the bowl — this is the point the optimizer converged to, and it is the physical deflection shape the beam actually adopts. Everywhere else on the surface represents a kinematically possible but non-equilibrium shape that would require some external agent to hold it there; only at the bottom does the structure balance on its own. This single picture is the clearest visual proof of the minimum-potential-energy principle: equilibrium literally sits at the lowest point of the energy surface.

Figure 2 — deflection curve comparison. Three predictions are overlaid: the exact analytical curve (solid black), the 2-term Rayleigh-Ritz curve (dashed blue), and the finite-element nodal solution (red dots). All three should be essentially indistinguishable, sagging maximally at midspan (around 28 mm for the given parameters) and returning to zero at both supports. The console output prints the relative error of each approximate method against the exact solution — expect the Ritz error to be under 0.1% and the FEM error to be smaller still, which shows that even a very small number of well-chosen basis functions, or a modest FEM mesh, is enough to essentially reproduce the exact elasticity solution once you frame the problem as energy minimization.

Figure 3 — deflection surface across load intensities. This second 3D plot sweeps the uniform load $q_0$ from light to heavy and stacks the resulting deflection curves into one continuous surface. Position along the beam runs along one horizontal axis, load intensity along the other, and deflection is the height. Because $w(x)$ is linear in $q_0$ for a linear-elastic beam, the surface is a smooth, uniformly tilted “wedge” — every cross-section at a fixed $x$ is a straight ramp in $q_0$, and every cross-section at a fixed $q_0$ reproduces the familiar sagging curve from Figure 2. This surface is a convenient way to read off, at a glance, how sensitive the beam’s midspan deflection is to load magnitude without re-running the calculation for each case.

7. Takeaways

Framing beam bending as an energy-minimization problem does more than reproduce the textbook formula — it generalizes cleanly. The same $\Pi[w] = U[w] - W[w]$ recipe extends directly to beams with variable cross-sections, non-uniform loads, multiple spans, or elastic foundations, simply by changing what goes into the integrals. The Rayleigh-Ritz method shows the idea in its purest, lowest-dimensional form, while the finite-element approach shows how the same principle scales up to realistic structural models — and, because the energy is quadratic, always boils down to nothing more exotic than solving $\mathbf{K}\mathbf{a}=\mathbf{F}$.

Finding the Shortest Way to Connect Three Points

Imagine you need to lay pipes, wires, or roads connecting three fixed locations to a single junction box, and you want to use as little material as possible. Where should that junction go? This is the classic Fermat–Torricelli point problem: given three points in a plane, find the point that minimizes the sum of distances to all three.

What makes this problem especially beautiful is that it has a direct mechanical analogy. If you drill three holes at the vertex positions on a horizontal board, thread a string through each hole, tie all three strings together in a single knot above the board, and hang equal weights from the other end of each string over the edge, the knot will settle exactly at the Fermat point. Physics finds the minimum for you.

The Math Behind It

Given three points $P_1, P_2, P_3$, we want to find point $P = (x, y)$ that minimizes:

$$
L(P) = \sum_{i=1}^{3} \lVert P - P_i \rVert
$$

Taking the gradient of $L$ with respect to $P$:

$$
\nabla L(P) = \sum_{i=1}^{3} \frac{P - P_i}{\lVert P - P_i \rVert}
$$

Each term is a unit vector pointing away from $P_i$ — physically, this is exactly the tension force of a string under equal load. At the minimum, these forces must cancel:

$$
\sum_{i=1}^{3} \hat{u}_i = 0
$$

Since all three vectors have the same magnitude (1) and must sum to zero, they can only be arranged one way: at 120° angles to each other. This is the famous geometric property of the Fermat point — as long as no angle of the triangle exceeds 120°, the point that satisfies this force-balance condition lies inside the triangle.

Setting $\nabla L(P) = 0$ and rearranging algebraically gives a fixed-point equation:

$$
P = \frac{\displaystyle\sum_{i=1}^{3} \frac{P_i}{\lVert P - P_i \rVert}}{\displaystyle\sum_{i=1}^{3} \frac{1}{\lVert P - P_i \rVert}}
$$

Iterating this equation is known as Weiszfeld’s algorithm, a classic and very fast way to solve this kind of “geometric median” problem.

Example Setup

Let’s use three concrete points, imagining them as three factories that need to be connected to a shared distribution hub:

• $A = (0, 0)$
• $B = (8, 0)$
• $C = (3, 7)$

We’ll solve for the Fermat point two different ways:

1. Naive relaxation — a direct simulation of the physical system: the “knot” is dragged step by step in the direction of the net string tension (an overdamped, no-inertia model of the real strings-and-weights setup).
2. Weiszfeld’s algorithm — the fast, mathematically derived fixed-point solver.

Comparing the two shows both why the mechanical picture works and how much faster a properly derived numerical method is than literally simulating physics step by step.

Full Source Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D  # enables 3D projection

# ---------------------------------------------------------
# 1. Problem setup: three fixed points (e.g., three factories
#    that must be connected to a single junction by pipes)
# ---------------------------------------------------------
points = np.array([
    [0.0, 0.0],   # Point A
    [8.0, 0.0],   # Point B
    [3.0, 7.0],   # Point C
])
labels = ['A', 'B', 'C']

# ---------------------------------------------------------
# 2. Objective: total length of the three "strings"
# ---------------------------------------------------------
def total_length(P, pts=points):
    P = np.asarray(P, dtype=float)
    return np.sum(np.linalg.norm(pts - P, axis=1))

# ---------------------------------------------------------
# 3. Gradient = net force exerted by the three unit-tension strings
# ---------------------------------------------------------
def net_force(P, pts=points, eps=1e-12):
    P = np.asarray(P, dtype=float)
    diff = P - pts
    dist = np.linalg.norm(diff, axis=1, keepdims=True)
    dist = np.maximum(dist, eps)
    unit_vectors = diff / dist
    return np.sum(unit_vectors, axis=0)

# ---------------------------------------------------------
# 4. Method 1: naive overdamped relaxation (direct simulation
#    of the physical system dragging the knot toward equilibrium)
# ---------------------------------------------------------
def relax_naive(x0, lr=0.03, n_iter=400):
    path = np.zeros((n_iter + 1, 2))
    P = np.array(x0, dtype=float)
    path[0] = P
    for i in range(n_iter):
        P = P - lr * net_force(P)
        path[i + 1] = P
    return path

# ---------------------------------------------------------
# 5. Method 2: Weiszfeld's algorithm (fast fixed-point solver,
#    derived by directly solving net_force(P) = 0 for P)
# ---------------------------------------------------------
def relax_weiszfeld(x0, pts=points, tol=1e-12, max_iter=200):
    P = np.array(x0, dtype=float)
    path = [P.copy()]
    for i in range(max_iter):
        dist = np.linalg.norm(pts - P, axis=1)
        dist = np.maximum(dist, 1e-12)
        w = 1.0 / dist
        P_new = np.sum(pts * w[:, None], axis=0) / np.sum(w)
        path.append(P_new.copy())
        if np.linalg.norm(P_new - P) < tol:
            P = P_new
            break
        P = P_new
    return np.array(path)

# ---------------------------------------------------------
# 6. Run both methods from the centroid
# ---------------------------------------------------------
x0 = points.mean(axis=0)
path_naive = relax_naive(x0)
path_fast  = relax_weiszfeld(x0)

fermat_naive = path_naive[-1]
fermat_fast  = path_fast[-1]

# ---------------------------------------------------------
# 7. Verify the 120-degree property
# ---------------------------------------------------------
def angle_between(v1, v2):
    c = np.dot(v1, v2) / (np.linalg.norm(v1) * np.linalg.norm(v2))
    return np.degrees(np.arccos(np.clip(c, -1.0, 1.0)))

vecs = points - fermat_fast
ang_AB = angle_between(vecs[0], vecs[1])
ang_BC = angle_between(vecs[1], vecs[2])
ang_CA = angle_between(vecs[2], vecs[0])

# ---------------------------------------------------------
# 8. Console report
# ---------------------------------------------------------
print("=" * 55)
print("Fermat-Torricelli point")
print("=" * 55)
print(f"A={points[0]}, B={points[1]}, C={points[2]}")
print(f"Naive relaxation   -> point={fermat_naive}, "
      f"iterations={len(path_naive)-1}, total length={total_length(fermat_naive):.6f}")
print(f"Weiszfeld (fast)   -> point={fermat_fast}, "
      f"iterations={len(path_fast)-1}, total length={total_length(fermat_fast):.6f}")
print(f"Angle A-F-B = {ang_AB:.3f} deg")
print(f"Angle B-F-C = {ang_BC:.3f} deg")
print(f"Angle C-F-A = {ang_CA:.3f} deg")
print(f"Sum of angles = {ang_AB+ang_BC+ang_CA:.3f} deg (should be ~360)")

# ---------------------------------------------------------
# 9. 2D visualization
# ---------------------------------------------------------
fig1, ax1 = plt.subplots(figsize=(7, 7))
tri = plt.Polygon(points, closed=True, fill=False,
                   edgecolor='gray', linestyle='--', linewidth=1.5)
ax1.add_patch(tri)

for p, lab in zip(points, labels):
    ax1.plot(*p, 'o', color='black', markersize=9)
    ax1.annotate(lab, p, textcoords="offset points",
                 xytext=(8, 8), fontsize=13, fontweight='bold')

for p in points:
    ax1.plot([fermat_fast[0], p[0]], [fermat_fast[1], p[1]],
              color='crimson', linewidth=2, zorder=3)

ax1.plot(path_naive[:, 0], path_naive[:, 1],
          color='royalblue', linestyle=':', linewidth=1.2,
          label=f'Naive relaxation ({len(path_naive)-1} steps)')
ax1.plot(path_fast[:, 0], path_fast[:, 1],
          color='darkorange', linestyle='-', marker='.', markersize=4,
          linewidth=1.2, label=f'Weiszfeld ({len(path_fast)-1} steps)')

ax1.plot(*x0, 's', color='green', markersize=9, label='Start (centroid)')
ax1.plot(*fermat_fast, '*', color='crimson', markersize=20,
          label='Fermat point', zorder=4)

ax1.set_title(f'Fermat-Torricelli point  (min total length = '
              f'{total_length(fermat_fast):.4f})')
ax1.set_xlabel('x')
ax1.set_ylabel('y')
ax1.set_aspect('equal')
ax1.grid(alpha=0.3)
ax1.legend(loc='upper right', fontsize=9)
plt.tight_layout()
plt.savefig('fermat_2d.png', dpi=150)
plt.show()

# ---------------------------------------------------------
# 10. 3D visualization of the total-length surface
# ---------------------------------------------------------
margin = 3.0
x_min, x_max = points[:, 0].min() - margin, points[:, 0].max() + margin
y_min, y_max = points[:, 1].min() - margin, points[:, 1].max() + margin

xs = np.linspace(x_min, x_max, 150)
ys = np.linspace(y_min, y_max, 150)
X, Y = np.meshgrid(xs, ys)

Z = np.zeros_like(X)
for p in points:
    Z += np.sqrt((X - p[0]) ** 2 + (Y - p[1]) ** 2)

fig2 = plt.figure(figsize=(9, 7))
ax2 = fig2.add_subplot(111, projection='3d')
surf = ax2.plot_surface(X, Y, Z, cmap='viridis', alpha=0.85,
                          linewidth=0, antialiased=True)

z_min = total_length(fermat_fast)
ax2.scatter(fermat_fast[0], fermat_fast[1], z_min,
             color='red', s=90, label='Fermat point (minimum)')
ax2.plot([fermat_fast[0]] * 2, [fermat_fast[1]] * 2,
          [Z.min(), z_min], color='red', linestyle='--')

for p, lab in zip(points, labels):
    zp = total_length(p)
    ax2.scatter(p[0], p[1], zp, color='black', s=50)
    ax2.text(p[0], p[1], zp, lab, fontsize=11)

ax2.set_xlabel('x')
ax2.set_ylabel('y')
ax2.set_zlabel('Total length L(x, y)')
ax2.set_title('Total-length surface and the Fermat point')
fig2.colorbar(surf, shrink=0.5, aspect=10, label='L(x, y)')
plt.tight_layout()
plt.savefig('fermat_3d.png', dpi=150)
plt.show()

Code Walkthrough

Section 1–2 (setup and objective): We define the three points as a NumPy array and write total_length, which computes $L(P)$ using np.linalg.norm with axis=1 — this subtracts the candidate point from all three points at once and returns their norms in a single vectorized call, rather than looping in Python.

Section 3 (net_force): This function computes $\nabla L(P)$. diff holds the three vectors from each fixed point to the candidate point $P$; dividing each by its own length turns it into a unit vector — exactly the tension direction of a string under load. Summing them gives the net physical force acting on the imaginary knot. We clip the distance with np.maximum(dist, eps) purely to avoid a division-by-zero if $P$ ever lands exactly on one of the vertices.

Section 4 (naive relaxation): This is the literal simulation of the physical system: at each time step, the knot moves a small amount (lr) in the direction of the net force, just like an overdamped object being pulled by three strings with no inertia. It takes hundreds of small steps to settle down, which is why we call it “naive” — it mirrors real, gradual physical motion.

Section 5 (Weiszfeld’s algorithm): Instead of nudging the point step by step, this directly applies the fixed-point formula derived earlier from $\nabla L(P) = 0$. Each iteration recomputes $P$ as a distance-weighted average of the three fixed points — points that are closer pull harder, points that are farther pull more gently, exactly balancing at the solution. This converges in a tiny fraction of the iterations that naive relaxation needs.

Section 6–7 (running and verifying): We run both methods starting from the centroid of the triangle and confirm that they land on essentially the same point. We then measure the angles between the three lines connecting the Fermat point to $A$, $B$, and $C$ — mathematically, these should each equal 120°, confirming the force-balance condition we derived earlier.

Section 9 (2D plot): This draws the triangle in dashed gray, marks the three fixed points, draws solid crimson lines from the Fermat point to each vertex (the “strings”), and overlays both convergence paths — the naive relaxation (dotted blue, many points) and Weiszfeld’s path (orange, very few points) — so you can visually see how much faster the fixed-point method converges.

Section 10 (3D plot): Here we build a grid of $(x, y)$ values covering the area around the triangle and compute $L(x, y)$ at every grid point in a fully vectorized way (looping only 3 times — once per fixed point — never once per grid cell). The resulting surface is a bowl-shaped function with a single global minimum, which we mark in red at the Fermat point. This is a great visual way to see that the optimization problem is well-behaved: one smooth valley, one clear minimum, no risk of getting stuck in a wrong answer.

About Performance

Both methods here are extremely lightweight — even a few hundred iterations on two-dimensional vectors run in well under a millisecond, so no further speed optimization is strictly necessary for this three-point example. That said, the code already demonstrates the general principle for scaling up: total_length, net_force, and the surface computation in Section 10 are all vectorized with NumPy rather than using per-point Python loops, and Weiszfeld’s algorithm converges in a small constant number of steps regardless of the geometry, so the same code would scale comfortably even if you needed to repeat this calculation for thousands of different triangles (for example, optimizing hub placement across many separate clusters of three locations at once).

Visualizing the Result

Run the code above in a fresh cell. It will print a summary to the console and generate two figures.

2D result — the triangle, the strings, and both convergence paths:

3D result — the total-length surface with the minimum marked:

Console output — point coordinates, iteration counts, and the 120° angle check:

=======================================================
Fermat-Torricelli point
=======================================================
A=[0. 0.], B=[8. 0.], C=[3. 7.]
Naive relaxation   -> point=[3.34728069 2.2650171 ], iterations=400, total length=13.964065
Weiszfeld (fast)   -> point=[3.34017112 2.26202749], iterations=43, total length=13.964055
Angle A-F-B = 120.000 deg
Angle B-F-C = 120.000 deg
Angle C-F-A = 120.000 deg
Sum of angles = 360.000 deg (should be ~360)

Once you drop in your own screenshots and console text above, look closely at the 2D plot: you should see the orange Weiszfeld path collapse onto the Fermat point in just a handful of steps, while the blue dotted naive-relaxation path spirals in much more gradually — a nice visual confirmation that solving the force-balance equation directly is far more efficient than simulating the physical settling process step by step, even though both arrive at exactly the same 120°-balanced point.

Maximizing Payload Fraction Under the Tsiolkovsky Equation

Introduction

Every rocket designer eventually runs into the same brutal arithmetic: the Tsiolkovsky rocket equation punishes you exponentially for wanting more velocity change. Add propellant, and you also add structure to hold that propellant, which in turn demands even more propellant to move. Multi-stage rockets exist precisely to escape this trap — by discarding dead structural mass along the way, you avoid dragging spent tanks and engines all the way to orbit.

But staging introduces a new question: if a mission requires a fixed total velocity change (Δv), how should that Δv be split among the stages? Give too much to a low-efficiency booster and you waste propellant; give too much to a high-efficiency upper stage and you carry excess structure through the atmosphere unnecessarily. This is the optimal staging problem, and today we’ll solve it concretely: for a fixed target Δv, we want the Δv allocation across stages that maximizes the payload fraction (equivalently, minimizes propellant and structural mass for a given payload).

Mathematical Formulation

The Tsiolkovsky equation

For a single stage, the velocity change is:

$$
\Delta v = I_{sp} , g_0 , \ln\left(\frac{m_0}{m_f}\right)
$$

where $I_{sp}$ is the specific impulse, $g_0$ is standard gravity, and $m_0/m_f$ is the mass ratio $R$.

Stage payload fraction

Each stage $i$ has a structural coefficient $\sigma_i$, the fraction of that stage’s own hardware+propellant mass that is dead weight (tanks, engines, avionics) rather than propellant:

$$
\sigma_i = \frac{m_{struct,i}}{m_{struct,i} + m_{prop,i}}
$$

If $\pi_i = m_{PL,i}/m_{0,i}$ is the fraction of stage $i$’s total initial mass that is “payload” (everything riding above it, including later stages), then the mass ratio and payload fraction are linked by:

$$
R_i = \frac{1}{\sigma_i(1-\pi_i) + \pi_i}
\quad\Longrightarrow\quad
\pi_i = \frac{\dfrac{1}{R_i} - \sigma_i}{1 - \sigma_i}
$$

Because each stage’s payload is the assembly of every stage above it, the fractions telescope: the fraction of the total liftoff mass that reaches orbit as pure payload is simply the product across all stages:

$$
\Pi = \prod_{i=1}^{N} \pi_i(\Delta v_i)
$$

The optimization problem

$$
\max_{\Delta v_1, \dots, \Delta v_N} ; \Pi = \prod_{i=1}^{N} \pi_i(\Delta v_i)
\qquad \text{subject to} \qquad
\sum_{i=1}^{N} \Delta v_i = \Delta v_{total}
$$

Taking logs turns the product into a sum, which is far friendlier for numerical solvers:

$$
\max ; \sum_{i=1}^{N} \ln \pi_i(\Delta v_i)
$$

At the optimum, Lagrange’s condition tells us the marginal effectiveness of adding one more m/s of Δv must be equal across every stage:

$$
\frac{\partial \ln \pi_i}{\partial \Delta v_i} = \lambda \quad \text{(same } \lambda \text{ for all } i\text{)}
$$

This is the numerical fingerprint we’ll check to confirm the solver actually found the true optimum.

A Concrete Example: 3-Stage Orbital Launcher

We target a total Δv of 9500 m/s — a realistic figure for reaching low Earth orbit once gravity and drag losses are folded in. Our three stages use progressively more efficient (and more expensive) propulsion:

Intuitively, the high-$I_{sp}$ hydrolox stage should be “worth more” Δv than the lower-$I_{sp}$ booster, but its higher structural coefficient pulls in the opposite direction. The optimizer resolves this trade-off numerically.

Python Source Code (Google Colaboratory)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169

# ============================================================
# Optimal Multi-Stage Rocket Staging Problem
# Maximizing Payload Fraction under the Tsiolkovsky Equation
# ============================================================

import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D  # noqa: F401 (3D projection registration)

plt.style.use('dark_background')

# ------------------------------------------------------------
# 1. Physical constants and rocket specifications
# ------------------------------------------------------------
G0 = 9.80665          # standard gravity [m/s^2]
DV_TOTAL = 9500.0     # required total delta-v to reach LEO [m/s]

# Each stage: (Isp [s], structural coefficient sigma = dead mass / (dead mass + propellant mass))
ISP    = np.array([282.0, 311.0, 450.0])   # kerolox booster, kerolox 2nd stage, hydrolox upper stage
SIGMA  = np.array([0.06, 0.08, 0.11])
N_STAGES = len(ISP)

# ------------------------------------------------------------
# 2. Core physics: Tsiolkovsky mass ratio and stage payload fraction
# ------------------------------------------------------------
def mass_ratio(dv, isp):
    """R_i = m0_i / mf_i  from the Tsiolkovsky rocket equation."""
    return np.exp(dv / (isp * G0))

def stage_payload_fraction(dv, isp, sigma):
    """
    pi_i = m_PL,i / m0,i for a single stage.
    Derived from R_i = 1 / (sigma*(1-pi_i) + pi_i).
    """
    r = mass_ratio(dv, isp)
    pi = (1.0 / r - sigma) / (1.0 - sigma)
    return pi

def total_payload_fraction(dv_vec, isp_vec=ISP, sigma_vec=SIGMA):
    """Telescoping product of all stage payload fractions."""
    pis = stage_payload_fraction(dv_vec, isp_vec, sigma_vec)
    return np.prod(pis)

# ------------------------------------------------------------
# 3. Optimization: maximize total payload fraction
#    subject to sum(dv_i) = DV_TOTAL
# ------------------------------------------------------------
def negative_log_payload(dv_vec):
    pis = stage_payload_fraction(dv_vec, ISP, SIGMA)
    if np.any(pis <= 1e-9):
        return 1e6  # heavy penalty outside the feasible domain
    return -np.sum(np.log(pis))

dv_guess_equal = np.full(N_STAGES, DV_TOTAL / N_STAGES)

constraints = {'type': 'eq', 'fun': lambda x: np.sum(x) - DV_TOTAL}
bounds = [(10.0, DV_TOTAL - 20.0)] * N_STAGES

result = minimize(
    negative_log_payload,
    dv_guess_equal,
    method='SLSQP',
    bounds=bounds,
    constraints=constraints,
    options={'maxiter': 300, 'ftol': 1e-14}
)

dv_opt = result.x
pi_opt_stages = stage_payload_fraction(dv_opt, ISP, SIGMA)
Pi_opt = np.prod(pi_opt_stages)
Pi_equal = total_payload_fraction(dv_guess_equal)

# ------------------------------------------------------------
# 4. Verify optimality: marginal effectiveness of each stage
#    (d ln(pi_i) / d dv_i should be equal across stages at the optimum)
# ------------------------------------------------------------
def marginal_effectiveness(dv_vec, isp_vec, sigma_vec, h=1.0):
    lam = np.zeros_like(dv_vec)
    for i in range(len(dv_vec)):
        dv_plus = dv_vec.copy(); dv_plus[i] += h
        dv_minus = dv_vec.copy(); dv_minus[i] -= h
        pi_plus = stage_payload_fraction(dv_plus[i], isp_vec[i], sigma_vec[i])
        pi_minus = stage_payload_fraction(dv_minus[i], isp_vec[i], sigma_vec[i])
        lam[i] = (np.log(pi_plus) - np.log(pi_minus)) / (2 * h)
    return lam

lambda_i = marginal_effectiveness(dv_opt, ISP, SIGMA)

# ------------------------------------------------------------
# 5. Report
# ------------------------------------------------------------
print("=" * 60)
print("OPTIMAL STAGING RESULT")
print("=" * 60)
for i in range(N_STAGES):
    print(f"Stage {i+1}: dv = {dv_opt[i]:8.2f} m/s | "
          f"pi_{i+1} = {pi_opt_stages[i]:.6f} | "
          f"marginal d(ln pi)/d(dv) = {lambda_i[i]:.8f}")
print("-" * 60)
print(f"Total delta-v check      : {np.sum(dv_opt):.3f} m/s (target {DV_TOTAL})")
print(f"Optimal payload fraction : {Pi_opt*100:.4f} %")
print(f"Equal-split payload frac.: {Pi_equal*100:.4f} %")
print(f"Improvement over equal split: {(Pi_opt/Pi_equal - 1)*100:.2f} %")
print("=" * 60)

# ------------------------------------------------------------
# 6. Visualization
# ------------------------------------------------------------
fig = plt.figure(figsize=(15, 6.5))

# --- 6a. 3D surface: payload fraction landscape over (dv1, dv2) ---
ax1 = fig.add_subplot(1, 2, 1, projection='3d')

n_grid = 70
margin = 80.0
dv1_range = np.linspace(margin, DV_TOTAL - 2 * margin, n_grid)
dv2_range = np.linspace(margin, DV_TOTAL - 2 * margin, n_grid)
DV1, DV2 = np.meshgrid(dv1_range, dv2_range)
DV3 = DV_TOTAL - DV1 - DV2

valid = DV3 > margin

R1 = mass_ratio(DV1, ISP[0]); PI1 = (1.0 / R1 - SIGMA[0]) / (1.0 - SIGMA[0])
R2 = mass_ratio(DV2, ISP[1]); PI2 = (1.0 / R2 - SIGMA[1]) / (1.0 - SIGMA[1])
R3 = mass_ratio(DV3, ISP[2]); PI3 = (1.0 / R3 - SIGMA[2]) / (1.0 - SIGMA[2])

Z = PI1 * PI2 * PI3 * 100.0
Z = np.where(valid & (PI1 > 0) & (PI2 > 0) & (PI3 > 0), Z, np.nan)

surf = ax1.plot_surface(DV1, DV2, Z, cmap='plasma', edgecolor='none',
                         alpha=0.92, antialiased=True)
ax1.scatter([dv_opt[0]], [dv_opt[1]], [Pi_opt * 100.0],
            color='cyan', s=90, marker='o', depthshade=False,
            edgecolor='white', linewidth=1.2, label='Optimum')

ax1.set_xlabel('Δv₁ (Stage 1) [m/s]', labelpad=10)
ax1.set_ylabel('Δv₂ (Stage 2) [m/s]', labelpad=10)
ax1.set_zlabel('Payload fraction [%]', labelpad=10)
ax1.set_title('Payload Fraction Landscape\n(Δv₃ = Δv_total − Δv₁ − Δv₂)', pad=15)
ax1.legend(loc='upper left')
fig.colorbar(surf, ax=ax1, shrink=0.55, pad=0.1, label='Payload fraction [%]')

# --- 6b. Bar chart: optimal vs. equal delta-v split ---
ax2 = fig.add_subplot(1, 2, 2)
x_pos = np.arange(N_STAGES)
width = 0.35

ax2.bar(x_pos - width/2, dv_guess_equal, width, label='Equal split',
        color='#888888', edgecolor='white')
ax2.bar(x_pos + width/2, dv_opt, width, label='Optimal split',
        color='#00d0ff', edgecolor='white')

ax2.set_xticks(x_pos)
ax2.set_xticklabels([f'Stage {i+1}' for i in range(N_STAGES)])
ax2.set_ylabel('Δv allocated [m/s]')
ax2.set_title(f'Δv Allocation Comparison\n'
              f'(Payload fraction: {Pi_equal*100:.3f}% → {Pi_opt*100:.3f}%)')
ax2.legend()
ax2.grid(axis='y', alpha=0.3)

for i in range(N_STAGES):
    ax2.text(x_pos[i] - width/2, dv_guess_equal[i] + 80,
              f'{dv_guess_equal[i]:.0f}', ha='center', fontsize=9)
    ax2.text(x_pos[i] + width/2, dv_opt[i] + 80,
              f'{dv_opt[i]:.0f}', ha='center', fontsize=9)

plt.tight_layout()
plt.show()

Code Walkthrough

Section 1 — Constants and specs. We fix standard gravity, the total mission Δv, and pack each stage’s $I_{sp}$ and $\sigma$ into small NumPy arrays. Keeping these as arrays (rather than separate variables) lets every downstream function operate on all three stages at once without loops.

Section 2 — Physics core. mass_ratio is a direct implementation of the Tsiolkovsky equation solved for $R = e^{\Delta v/(I_{sp} g_0)}$. stage_payload_fraction inverts the relationship derived above to recover $\pi_i$ from a chosen $\Delta v_i$. total_payload_fraction multiplies the per-stage fractions together, implementing the telescoping product $\Pi = \prod \pi_i$.

Section 3 — The optimizer. Rather than maximizing the product directly (numerically unstable when payload fractions are small), we minimize the negative sum of logs, which is mathematically equivalent but far more numerically stable. SLSQP (Sequential Least Squares Programming) is the right choice here because it natively supports both bounds and an equality constraint — exactly what “$\sum \Delta v_i = \Delta v_{total}$” requires. The penalty branch (return 1e6) protects the optimizer from ever evaluating the objective at an infeasible point where a stage’s payload fraction would go negative, which would otherwise throw a math domain error from log.

Section 4 — Optimality verification. This is the empirical check of the Lagrange condition. We nudge each $\Delta v_i$ by ±1 m/s and take a central difference of $\ln \pi_i$. If the optimizer truly found the constrained optimum, all three marginal effectiveness values $\lambda_i$ should converge to (nearly) the same number — proof that no further improvement is possible by shifting Δv between any pair of stages.

Section 5 — Reporting. Plain print statements summarize the optimal allocation, the resulting payload fraction, and the percentage improvement over a naive equal three-way split of Δv.

Section 6 — Visualization. The 3D surface sweeps $\Delta v_1$ and $\Delta v_2$ over a grid (with $\Delta v_3$ determined by the constraint), computing the resulting payload fraction at every grid point — fully vectorized with NumPy broadcasting rather than nested Python loops, so the entire 70×70 grid evaluates essentially instantly. The optimum found by SLSQP is plotted as a highlighted marker sitting on the peak of the surface, visually confirming the numerical result. The bar chart puts the optimal and equal-split Δv allocations side by side for a direct, human-readable comparison.

Runtime Performance

This problem has only three decision variables and a smooth, well-behaved objective, so SLSQP converges in well under a second — no acceleration techniques are needed for the optimization itself. The only part that touches any real volume of computation is the 3D surface grid, and vectorizing it with NumPy array operations (as done above, with zero Python-level loops over grid points) keeps that well under a second as well.

============================================================
OPTIMAL STAGING RESULT
============================================================
Stage 1: dv =  1110.04 m/s | pi_1 = 0.648284 | marginal d(ln pi)/d(dv) = -0.00039720
Stage 2: dv =  2378.88 m/s | pi_2 = 0.411314 | marginal d(ln pi)/d(dv) = -0.00039720
Stage 3: dv =  6011.08 m/s | pi_3 = 0.164173 | marginal d(ln pi)/d(dv) = -0.00039720
------------------------------------------------------------
Total delta-v check      : 9500.000 m/s (target 9500.0)
Optimal payload fraction : 4.3776 %
Equal-split payload frac.: 3.4746 %
Improvement over equal split: 25.99 %
============================================================

Interpreting the Graphs

The 3D surface (left panel) shows the payload fraction as a smooth hill over the $(\Delta v_1, \Delta v_2)$ plane, with $\Delta v_3$ implicitly filling in whatever Δv remains. The surface falls off sharply near the edges — where one stage is forced to carry almost all of the mission’s Δv and its mass ratio blows up exponentially — and rises to a single interior peak. The cyan marker sits exactly on that peak, confirming that SLSQP located the true maximum rather than a boundary artifact or local irregularity.

The bar chart (right panel) translates the abstract optimum into an engineering decision: how many m/s of Δv should each real stage actually be designed to deliver? Because the hydrolox third stage has the highest $I_{sp}$, expect the optimizer to shift more of the mission’s total Δv onto it relative to a naive equal split — the exponential penalty for low-$I_{sp}$ mass ratios makes the lower stages relatively more “expensive” per unit of Δv, even though their structural coefficients are smaller.

Conclusion

The optimal staging problem is a clean example of constrained nonlinear optimization hiding inside a deceptively simple physical formula. By reformulating the payload-fraction product as a log-sum, applying SLSQP under an equality constraint, and cross-checking the result against the analytic Lagrange condition, we get both a numerically verified optimum and an intuitive visualization of the entire solution landscape. The same framework generalizes directly to four or more stages, asymmetric Δv-loss profiles (gravity and drag losses concentrated in the lower stages), or reusable first-stage penalties — all by extending the ISP and SIGMA arrays and letting the optimizer do the rest.

Maximizing Projectile Range with Air Resistance

Every physics student learns that a projectile launched over flat ground travels farthest when fired at 45°. That result, however, only holds in a vacuum. The moment air resistance enters the picture, the textbook formula falls apart — the optimal angle drifts below 45°, and there’s no closed-form expression left to solve for it. This is where the problem stops being a physics exercise and becomes a genuine numerical optimization task: given a nonlinear equation of motion, find the launch angle that maximizes horizontal range.

In this article, we’ll set up the drag-augmented projectile problem, solve it numerically in a way that stays fast even when scanning hundreds of angle/velocity combinations, and visualize the resulting range surface in 3D.

Problem Setup

Without drag, the range of a projectile launched at speed $v_0$ and angle $\theta$ is:

$$
R(\theta) = \frac{v_0^2 \sin(2\theta)}{g}
$$

which is maximized analytically at $\theta = 45°$.

With quadratic air drag, the equations of motion become a coupled nonlinear system with no closed-form solution:

$$
\frac{dv_x}{dt} = -k , |v| , v_x, \qquad \frac{dv_y}{dt} = -g - k , |v| , v_y
$$

$$
|v| = \sqrt{v_x^2 + v_y^2}, \qquad \frac{dx}{dt} = v_x, \qquad \frac{dy}{dt} = v_y
$$

Here $k$ is the drag coefficient divided by mass ($k = \frac{1}{2}\rho C_d A / m$), and $g = 9.81\ \text{m/s}^2$. Since $R(\theta)$ can only be evaluated by numerically integrating this system, finding the optimal angle means combining numerical ODE integration with numerical optimization.

A naive approach would integrate the trajectory for one angle at a time inside a Python loop, calling a solver like scipy.integrate.solve_ivp dozens or hundreds of times. That works, but it’s slow — each call carries fixed overhead, and scanning a full angle range or a 2D grid of (angle, velocity) combinations for a 3D plot multiplies that cost quickly. The code below avoids this by integrating all trajectories simultaneously as vectorized NumPy arrays, so one function call handles an entire batch of launch conditions at once.

Source Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146