Solving Steady-State Heat Conduction with Energy Minimization

August 11, 2026

Many physical systems reach a steady (equilibrium) state not because something is being “solved” in the traditional sense, but because nature is minimizing a quantity — usually an energy functional. Steady-state heat conduction is a textbook example: instead of solving the heat equation directly, we can find the temperature distribution that minimizes a scalar energy functional. This is the Dirichlet principle, and it turns out the exact same mathematical machinery describes irrotational, incompressible fluid flow (potential flow), since both are governed by the Laplace/Poisson equation.

In this article, we’ll take a concrete example — a 2D metal plate heated from the left and containing an internal heat source — and find its steady-state temperature field by minimizing energy directly, rather than solving a linear PDE the “usual” way.

The Physics: From Energy to Equation

For steady-state heat conduction with conductivity $k$ and a heat source density $q(x,y)$, the temperature field $T(x,y)$ that the system settles into is the one that minimizes the functional:

$$
E[T] = \int_\Omega \left[ \frac{k}{2} , |\nabla T|^2 ;-; q(x,y), T \right] , dx, dy
$$

The first term is analogous to elastic/kinetic energy stored in the temperature gradient; the second term is the work done by the heat source. Taking the variational derivative and setting it to zero (Euler–Lagrange equation) gives:

$$
\frac{\delta E}{\delta T} = -k \nabla^2 T - q = 0 \quad \Longrightarrow \quad k,\nabla^2 T + q = 0
$$

This is exactly the steady-state heat equation. So minimizing $E[T]$ and solving $k\nabla^2 T + q = 0$ are mathematically equivalent — this is the essence of the Dirichlet principle.

The same functional form, with $T$ replaced by a velocity potential $\phi$ and $q=0$, is the kinetic energy of an incompressible, irrotational fluid:

$$
E[\phi] = \int_\Omega \frac{1}{2},|\nabla \phi|^2 , dx, dy, \qquad \nabla^2 \phi = 0
$$

which is why the code below applies equally well to potential-flow problems.

Discretizing the Problem

On a finite-difference grid, the continuous functional becomes a quadratic form in the vector of unknown (interior) temperatures $\mathbf{T}$:

$$
E(\mathbf{T}) = \frac{1}{2},\mathbf{T}^\top A, \mathbf{T} - \mathbf{f}^\top \mathbf{T}
$$

where $A$ is the discrete (negative) Laplacian operator (positive-definite, sparse) and $\mathbf{f}$ bundles the heat source and boundary-condition contributions. Its gradient is:

$$
\nabla E(\mathbf{T}) = A\mathbf{T} - \mathbf{f}
$$

Setting this to zero recovers the linear system $A\mathbf{T} = \mathbf{f}$ — but instead of solving it as a linear system, we’ll hand $E(\mathbf{T})$ and $\nabla E(\mathbf{T})$ to an optimizer and let it find the minimum directly, which is the whole point of this exercise.

The Example Problem

A square plate, $1\text{m} \times 1\text{m}$:

Left wall ($x=0$) held at $T = 100°C$ (a hot wall)
Top, bottom, and right walls held at $T = 0°C$
A localized internal heat source (a Gaussian “heater”) embedded inside the plate near $(0.65, 0.5)$

We’ll compute the steady-state temperature field two different ways — a general-purpose optimizer (L-BFGS-B) and a specialized, much faster energy-minimizing algorithm (Conjugate Gradient) — and cross-check both against a direct sparse linear solve.

Python Source Code

# ==========================================================
# Steady-State Heat Conduction via Energy Minimization
# ==========================================================
import numpy as np
import scipy.sparse as sp
import scipy.sparse.linalg as spla
from scipy.optimize import minimize
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import time

# ---------- 1. Grid and physical parameters ----------
Nx, Ny = 51, 51
Lx, Ly = 1.0, 1.0
dx = Lx / (Nx - 1)
dy = Ly / (Ny - 1)
x = np.linspace(0, Lx, Nx)
y = np.linspace(0, Ly, Ny)
X, Y = np.meshgrid(x, y, indexing='ij')

k_cond = 1.0      # thermal conductivity
T_hot  = 100.0    # left wall temperature
T_cold = 0.0      # other walls temperature

# internal heat source (Gaussian "heater")
q0, xs, ys, sigma = 3000.0, 0.65, 0.5, 0.06
Q = q0 * np.exp(-((X - xs)**2 + (Y - ys)**2) / (2 * sigma**2))

# ---------- 2. Build the discrete energy operator A and load vector f ----------
mx, my = Nx - 2, Ny - 2   # number of interior nodes in each direction

main_x = -2.0 * np.ones(mx)
off_x  = np.ones(mx - 1)
D2x = sp.diags([off_x, main_x, off_x], [-1, 0, 1], format='csr') / dx**2

main_y = -2.0 * np.ones(my)
off_y  = np.ones(my - 1)
D2y = sp.diags([off_y, main_y, off_y], [-1, 0, 1], format='csr') / dy**2

Ix = sp.identity(mx, format='csr')
Iy = sp.identity(my, format='csr')

L = sp.kron(D2x, Iy, format='csr') + sp.kron(Ix, D2y, format='csr')
A = (-L).tocsr()  # positive-definite discrete energy operator

Q_int = Q[1:-1, 1:-1]
boundary_term = np.zeros((mx, my))
boundary_term[0, :] += T_hot / dx**2   # contribution from the hot left wall

f = (Q_int / k_cond + boundary_term).flatten()

# ---------- 3. Discrete energy functional and its gradient ----------
def energy(t):
    return 0.5 * t @ (A @ t) - f @ t

def grad_energy(t):
    return A @ t - f

# ---------- 4. Reference solution (direct sparse solve) ----------
t0 = time.time()
T_ref = spla.spsolve(A.tocsc(), f)
t_direct = time.time() - t0

# ---------- 5. Energy minimization: general-purpose optimizer (L-BFGS-B) ----------
energy_history_bfgs = []
def callback_bfgs(t):
    energy_history_bfgs.append(energy(t))

t0 = time.time()
res_bfgs = minimize(energy, np.zeros(mx * my), jac=grad_energy,
                     method='L-BFGS-B', callback=callback_bfgs,
                     options={'maxiter': 800, 'ftol': 1e-14, 'gtol': 1e-8})
t_bfgs = time.time() - t0
T_bfgs = res_bfgs.x

# ---------- 6. Energy minimization: fast specialized version (Conjugate Gradient) ----------
energy_history_cg = []
def callback_cg(t):
    energy_history_cg.append(energy(t))

t0 = time.time()
T_cg, info = spla.cg(A, f, rtol=1e-12, maxiter=1000, callback=callback_cg)
t_cg = time.time() - t0

print(f"Direct sparse solve : {t_direct*1000:8.2f} ms")
print(f"L-BFGS-B minimizer  : {t_bfgs*1000:8.2f} ms  ({len(energy_history_bfgs)} iterations)")
print(f"Conjugate Gradient  : {t_cg*1000:8.2f} ms  ({len(energy_history_cg)} iterations)")
print(f"max|T_bfgs - T_ref| = {np.max(np.abs(T_bfgs - T_ref)):.3e}")
print(f"max|T_cg   - T_ref| = {np.max(np.abs(T_cg   - T_ref)):.3e}")

# ---------- 7. Assemble the full temperature field (with boundary values) ----------
def assemble_full(t_int):
    T_full = np.zeros((Nx, Ny))
    T_full[-1, :] = T_cold
    T_full[:, 0]  = T_cold
    T_full[:, -1] = T_cold
    T_full[0, :]  = T_hot
    T_full[1:-1, 1:-1] = t_int.reshape(mx, my)
    return T_full

T_field = assemble_full(T_cg)

# ---------- 8. Heat flux field  q = -k*grad(T)  (a fluid-like transport field) ----------
dTdx, dTdy = np.gradient(T_field, dx, dy)
qx, qy = -k_cond * dTdx, -k_cond * dTdy

# ---------- 9. Visualization ----------
fig1, ax1 = plt.subplots(figsize=(7, 6))
c = ax1.contourf(X, Y, T_field, levels=40, cmap='inferno')
skip = 4
ax1.quiver(X[::skip, ::skip], Y[::skip, ::skip],
           qx[::skip, ::skip], qy[::skip, ::skip],
           color='white', scale=4000, width=0.003)
ax1.set_xlabel('x [m]')
ax1.set_ylabel('y [m]')
ax1.set_title('Steady-State Temperature Field and Heat Flux')
fig1.colorbar(c, ax=ax1, label='Temperature [°C]')
plt.tight_layout()
plt.show()

fig2, ax2 = plt.subplots(figsize=(7, 5))
ax2.plot(energy_history_bfgs, 'o-', label='L-BFGS-B', markersize=3)
ax2.plot(energy_history_cg, 's-', label='Conjugate Gradient', markersize=3)
ax2.set_xlabel('Iteration')
ax2.set_ylabel('Energy E[T]')
ax2.set_title('Convergence of the Energy Functional')
ax2.legend()
ax2.grid(alpha=0.3)
plt.tight_layout()
plt.show()

fig3 = plt.figure(figsize=(8, 6))
ax3 = fig3.add_subplot(111, projection='3d')
surf = ax3.plot_surface(X, Y, T_field, cmap='inferno', linewidth=0, antialiased=True)
ax3.set_xlabel('x [m]')
ax3.set_ylabel('y [m]')
ax3.set_zlabel('T [°C]')
ax3.set_title('3D Temperature Distribution')
fig3.colorbar(surf, ax=ax3, shrink=0.6, label='Temperature [°C]')
plt.tight_layout()
plt.show()

Direct sparse solve :    36.16 ms
L-BFGS-B minimizer  :  1754.44 ms  (222 iterations)
Conjugate Gradient  :    86.39 ms  (169 iterations)
max|T_bfgs - T_ref| = 8.850e-05
max|T_cg   - T_ref| = 8.624e-11

Code Walkthrough

Section 1 — Grid setup. We define a $51 \times 51$ grid over a unit square. X, Y = np.meshgrid(..., indexing='ij') keeps the array’s first axis aligned with $x$ and the second with $y$, which matters later when we flatten arrays into vectors.

Section 2 — Building $A$ and $\mathbf{f}$. This is the heart of the discretization. D2x and D2y are standard tridiagonal second-derivative operators (the 1D Laplacian stencil $[1, -2, 1]/h^2$) built only over the interior grid points — boundary points are not unknowns, they’re known data. We combine the two 1D operators into a 2D Laplacian using a Kronecker sum: sp.kron(D2x, Iy) + sp.kron(Ix, D2y). This is the standard trick for turning a 2D finite-difference stencil into a single sparse matrix without writing nested loops. We negate the Laplacian (A = -L) so that $A$ is positive-definite, matching the sign convention required for an energy-minimization (rather than energy-maximization) problem.

Because the hot boundary ($T=100°C$ at $x=0$) is not one of our unknowns, its influence has to be folded into the right-hand side vector $\mathbf{f}$ manually — this is the boundary_term array, which only affects the row of interior nodes immediately adjacent to the hot wall.

Section 3 — Energy and gradient functions. energy(t) directly implements $E(\mathbf{T}) = \frac12 \mathbf{T}^\top A \mathbf{T} - \mathbf{f}^\top \mathbf{T}$ and grad_energy(t) implements its exact analytic gradient $A\mathbf{T}-\mathbf{f}$. Supplying the exact gradient (rather than letting the optimizer estimate it by finite differences) is the single biggest performance factor here — without it, scipy.optimize.minimize would need $O(n)$ extra function evaluations per iteration just to approximate the gradient.

Section 4 — Reference solution. spla.spsolve solves $A\mathbf{T}=\mathbf{f}$ directly using sparse LU factorization. This is our ground truth to validate the optimization-based approaches against.

Sections 5 & 6 — Two flavors of energy minimization.

L-BFGS-B is a general-purpose quasi-Newton optimizer. It works for any differentiable energy functional (including nonlinear ones), which is why it’s the natural first choice pedagogically.
Conjugate Gradient (scipy.sparse.linalg.cg) is not just “a solver” — it is literally an algorithm that minimizes a quadratic energy functional of exactly this form, one conjugate direction at a time. Because it’s purpose-built for symmetric positive-definite quadratic problems (like ours) rather than general nonlinear ones, it converges in dramatically fewer, cheaper iterations. This is our “fast version”: same energy-minimization idea, specialized to the structure of the problem for a large speedup, while still being philosophically the same energy-minimization approach as L-BFGS-B.

Section 7 — Reassembling the full grid. The optimizer only ever sees the vector of interior unknowns; assemble_full puts the known boundary values back around it to produce a complete $(N_x, N_y)$ temperature field for plotting.

Section 8 — Heat flux. Fourier’s law says heat flows down the temperature gradient: $\mathbf{q} = -k\nabla T$. This vector field is mathematically identical in form to a fluid velocity field in potential flow, which is why we visualize it with arrows — it makes the heat-conduction/fluid-flow analogy visually concrete.

Why Conjugate Gradient Wins on Speed

L-BFGS-B has to build up an approximation to the inverse Hessian using gradient history, which costs extra bookkeeping every iteration and typically needs many iterations to reach high accuracy on an ill-conditioned quadratic like a discretized Laplacian. Conjugate Gradient, by contrast, exploits the fact that our energy is exactly quadratic with a known, sparse, symmetric positive-definite matrix — each iteration is just one sparse matrix-vector product plus a few dot products, and CG is mathematically guaranteed to reach the exact minimum within $n$ steps (and in practice, far fewer, since the Poisson matrix’s eigenvalues cluster). The console output above should show CG finishing in a small fraction of the time L-BFGS-B needs, while both agree with the direct solve to within numerical precision.

Understanding the Results

Figure 1 — Temperature field and heat flux. The heatmap shows temperature decaying smoothly from the hot left wall (bright) toward the cooler edges, with a distinct hot spot around the embedded heat source. The white arrows show the heat flux vector field $\mathbf{q}=-k\nabla T$ — heat flows from hot to cold, “downhill” on the temperature surface, exactly like a fluid flowing from high to low potential.

Figure 2 — Energy convergence. This plot shows the value of the energy functional $E(\mathbf{T})$ at each iteration for both algorithms. Both curves should descend monotonically toward the same minimum value (since $A$ is positive-definite, the quadratic has a unique global minimum), but the Conjugate Gradient curve should reach convergence in far fewer iterations — visually demonstrating that both algorithms are doing the same conceptual thing (rolling downhill on the energy landscape) at very different speeds.

Figure 3 — 3D temperature surface. This is the same data as Figure 1, but viewed as a literal energy landscape — high near the hot wall and the internal heater, sloping down toward the cold boundaries. Seeing it as a physical surface makes the “minimization” framing intuitive: the true physical steady state is the shape a stretched elastic membrane would settle into if pinned at the given boundary heights and pushed up by the heat source, minimizing its stored elastic energy.

Closing Thoughts

Framing steady-state heat conduction as an energy-minimization problem does more than provide an alternative numerical method — it reveals the same mathematics underlying incompressible potential flow, electrostatics, and membrane mechanics. Once a problem is expressed as minimizing $\frac12\mathbf{T}^\top A\mathbf{T} - \mathbf{f}^\top\mathbf{T}$, any tool from the optimization world — gradient descent, L-BFGS, conjugate gradient, or beyond — becomes a legitimate physics solver, and the choice between them becomes purely a question of computational efficiency rather than modeling correctness.

Minimizing the Weight of a Truss Structure with Python

August 10, 2026

A Two-Bar Truss Case Study

Structural weight minimization is one of the classic entry points into engineering optimization. The idea is simple to state but rich in behavior: given a structure that has to carry a load safely, find the combination of member cross-sectional areas and member lengths (or, equivalently, the geometry that determines those lengths) that uses the least material while still satisfying stress and stiffness requirements.

In this article we work through a concrete, fully worked example — a symmetric two-bar truss — set up the mechanics equations, formulate it as a constrained nonlinear optimization problem in Python, solve it with scipy.optimize, and then visualize the result from several angles, including a 3D view of the design space.

The Structure

Picture two support points fixed on the ground, separated by a fixed half-width $b$, and a single apex node above them where a vertical load $P$ is applied. Two bars connect the apex to each support, forming a symmetric “A-frame” truss. Both bars share the same cross-sectional area $A$.

The two design variables are:

$A$ — the cross-sectional area of each bar
$h$ — the height of the apex above the base

Because the base width $b$ is fixed, the bar length is fully determined by $h$:

$$
L(h) = \sqrt{b^2 + h^2}
$$

So treating $h$ as a design variable is equivalent to treating the member length directly as a design variable — raising the apex lengthens the bars and simultaneously changes the load angle.

The Mechanics

By vertical equilibrium at the apex, each bar carries an axial force:

$$
F(h) = \frac{P , L(h)}{2h}
$$

The axial stress in each bar is:

$$
\sigma(A,h) = \frac{F(h)}{A} = \frac{P , L(h)}{2Ah}
$$

The vertical deflection of the apex node under the load (from the bar’s elastic elongation, projected back through the load angle) is:

$$
\delta(A,h) = \frac{P , L(h)^3}{2AEh^2}
$$

where $E$ is Young’s modulus of the bar material.

The total structural weight is:

$$
W(A,h) = 2 \rho A L(h) = 2 \rho A \sqrt{b^2+h^2}
$$

where $\rho$ is the material density.

The Optimization Problem

We want to minimize weight subject to an allowable stress $\sigma_{allow}$ and an allowable deflection $\delta_{allow}$:

$$
\min_{A,,h} \quad W(A,h) = 2\rho A \sqrt{b^2+h^2}
$$

$$
\text{subject to} \quad \sigma(A,h) \le \sigma_{allow}, \qquad \delta(A,h) \le \delta_{allow}
$$

$$
A_{min} \le A \le A_{max}, \qquad h_{min} \le h \le h_{max}
$$

This is a genuinely interesting problem because increasing $h$ makes the load angle steeper (reducing the axial force, and therefore the area needed to satisfy the stress limit), but it also increases $L$, which pushes the deflection up as $L^3$. There is no free lunch: some intermediate height minimizes the total weight, and the optimizer has to find it.

Python Implementation

The code below is a single, self-contained cell. It sets up the mechanics, runs a gradient-based constrained optimizer (SLSQP), cross-checks the numerical result against a fast analytical reduction of the problem, and produces four plots.

# ==========================================================
# Weight Minimization of a Two-Bar Truss
# Design Variables: Cross-sectional Area (A) and Height (h)
# ==========================================================

import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt

# ----------------------------------------------------------
# 1. Fixed Physical Parameters
# ----------------------------------------------------------
P           = 50000.0     # Applied vertical load at apex [N]
b           = 1.0         # Fixed half base width [m]
rho         = 7850.0      # Density of steel [kg/m^3]
E           = 200e9       # Young's modulus of steel [Pa]
sigma_allow = 165e6       # Allowable axial stress [Pa]
delta_allow = 5e-3        # Allowable vertical deflection at apex [m]

# ----------------------------------------------------------
# 2. Geometry / Mechanics helper functions
# ----------------------------------------------------------
def bar_length(h):
    """Length of one bar as a function of apex height h."""
    return np.sqrt(b**2 + h**2)

def axial_force(h):
    """Axial force carried by one bar (symmetric 2-bar truss)."""
    L = bar_length(h)
    return P * L / (2.0 * h)

def axial_stress(A, h):
    """Axial stress in one bar."""
    return axial_force(h) / A

def apex_deflection(A, h):
    """Vertical deflection at the apex node."""
    L = bar_length(h)
    return (P * L**3) / (2.0 * A * E * h**2)

def total_weight(x):
    """Total weight of the two bars: x = [A, h]."""
    A, h = x
    L = bar_length(h)
    return 2.0 * rho * A * L

# ----------------------------------------------------------
# 3. Constraint functions (SLSQP requires g(x) >= 0 form)
# ----------------------------------------------------------
def stress_constraint(x):
    A, h = x
    return sigma_allow - axial_stress(A, h)

def deflection_constraint(x):
    A, h = x
    return delta_allow - apex_deflection(A, h)

# ----------------------------------------------------------
# 4. Optimization set-up
# ----------------------------------------------------------
A_min, A_max = 1.0e-5, 5.0e-3     # [m^2]  (0.1 cm^2 -- 50 cm^2)
h_min, h_max = 0.5, 4.0           # [m]

bounds = [(A_min, A_max), (h_min, h_max)]
x0     = np.array([5.0e-4, 1.5])   # initial guess [A, h]

history = [x0.copy()]
def record(xk):
    history.append(xk.copy())

result = minimize(
    total_weight,
    x0,
    method='SLSQP',
    bounds=bounds,
    constraints=[
        {'type': 'ineq', 'fun': stress_constraint},
        {'type': 'ineq', 'fun': deflection_constraint},
    ],
    callback=record,
    options={'maxiter': 200, 'ftol': 1e-10}
)

A_opt, h_opt = result.x
W_opt        = result.fun
history      = np.array(history)

# ----------------------------------------------------------
# 5. Report
# ----------------------------------------------------------
print("=" * 55)
print(" OPTIMIZATION RESULT ")
print("=" * 55)
print(f"Success              : {result.success}")
print(f"Iterations           : {result.nit}")
print(f"Optimal area  A*     : {A_opt*1e4:.4f} cm^2")
print(f"Optimal height h*    : {h_opt:.4f} m")
print(f"Optimal bar length L*: {bar_length(h_opt):.4f} m")
print(f"Minimum weight W*    : {W_opt:.4f} kg")
print(f"Stress at optimum    : {axial_stress(A_opt,h_opt)/1e6:.3f} MPa "
      f"(allowable {sigma_allow/1e6:.1f} MPa)")
print(f"Deflection at optimum: {apex_deflection(A_opt,h_opt)*1e3:.4f} mm "
      f"(allowable {delta_allow*1e3:.1f} mm)")
print("=" * 55)

# ----------------------------------------------------------
# 6. Analytical cross-check: reduce to a 1-D problem in h
#    (since weight increases monotonically with A, the
#     optimal A for a given h is exactly the tighter of the
#     two constraints)
# ----------------------------------------------------------
h_fine = np.linspace(h_min, h_max, 2000)
L_fine = bar_length(h_fine)

A_stress_curve = P * L_fine / (2.0 * h_fine * sigma_allow)
A_defl_curve   = (P * L_fine**3) / (2.0 * E * h_fine**2 * delta_allow)
A_active_curve = np.maximum(A_stress_curve, A_defl_curve)
W_active_curve = 2.0 * rho * A_active_curve * L_fine

idx_min   = np.argmin(W_active_curve)
h_star_1d = h_fine[idx_min]
W_star_1d = W_active_curve[idx_min]

print(f"1-D cross-check   h* : {h_star_1d:.4f} m,  W*: {W_star_1d:.4f} kg")

# ----------------------------------------------------------
# 7. Figure 1: 3-D weight surface with optimization path
# ----------------------------------------------------------
A_grid = np.linspace(A_min, A_max, 120)
h_grid = np.linspace(h_min, h_max, 120)
AG, HG = np.meshgrid(A_grid, h_grid)
WG = 2.0 * rho * AG * bar_length(HG)

fig1 = plt.figure(figsize=(10, 7))
ax1 = fig1.add_subplot(111, projection='3d')
surf = ax1.plot_surface(AG*1e4, HG, WG, cmap='viridis',
                         alpha=0.75, linewidth=0, antialiased=True)

path_W = np.array([total_weight(pt) for pt in history])
ax1.plot(history[:, 0]*1e4, history[:, 1], path_W,
         color='red', marker='o', markersize=4, linewidth=2,
         label='Optimization path')
ax1.scatter([A_opt*1e4], [h_opt], [W_opt],
            color='gold', s=120, edgecolor='black',
            marker='*', label='Optimum', zorder=5)

ax1.set_xlabel(r'Area $A$ [cm$^2$]')
ax1.set_ylabel(r'Height $h$ [m]')
ax1.set_zlabel(r'Weight $W$ [kg]')
ax1.set_title('Weight Surface and Optimization Path')
fig1.colorbar(surf, shrink=0.6, aspect=12, label='Weight [kg]')
ax1.legend()
plt.tight_layout()
plt.show()

# ----------------------------------------------------------
# 8. Figure 2: 2-D contour with feasible region & boundaries
# ----------------------------------------------------------
feasible = ((AG >= P*bar_length(HG)/(2*HG*sigma_allow)) &
            (AG >= (P*bar_length(HG)**3)/(2*E*HG**2*delta_allow)))

fig2, ax2 = plt.subplots(figsize=(9, 7))
cont = ax2.contourf(AG*1e4, HG, WG, levels=30, cmap='viridis', alpha=0.85)
fig2.colorbar(cont, label='Weight [kg]')

ax2.contourf(AG*1e4, HG, feasible.astype(int), levels=[-0.5, 0.5],
             colors=['white'], alpha=0.6)
ax2.plot(A_stress_curve*1e4, h_fine, color='cyan', linewidth=2.5,
         label='Stress boundary')
ax2.plot(A_defl_curve*1e4, h_fine, color='magenta', linewidth=2.5,
         label='Deflection boundary')
ax2.plot(history[:, 0]*1e4, history[:, 1], color='red',
         marker='o', markersize=4, linewidth=2, label='Optimization path')
ax2.scatter([A_opt*1e4], [h_opt], color='gold', s=160,
            edgecolor='black', marker='*', zorder=5, label='Optimum')

ax2.set_xlim(A_min*1e4, A_max*1e4)
ax2.set_ylim(h_min, h_max)
ax2.set_xlabel(r'Area $A$ [cm$^2$]')
ax2.set_ylabel(r'Height $h$ [m]')
ax2.set_title('Design Space: Feasible Region and Weight Contours')
ax2.legend(loc='upper right')
plt.tight_layout()
plt.show()

# ----------------------------------------------------------
# 9. Figure 3: Convergence history
# ----------------------------------------------------------
fig3, ax3 = plt.subplots(figsize=(8, 5))
ax3.plot(range(len(path_W)), path_W, marker='o', color='steelblue')
ax3.axhline(W_opt, color='red', linestyle='--',
            label=f'Converged weight = {W_opt:.3f} kg')
ax3.set_xlabel('Iteration')
ax3.set_ylabel('Weight [kg]')
ax3.set_title('Convergence History of the Optimizer')
ax3.grid(alpha=0.3)
ax3.legend()
plt.tight_layout()
plt.show()

# ----------------------------------------------------------
# 10. Figure 4: 1-D analytical cross-check
# ----------------------------------------------------------
fig4, ax4 = plt.subplots(figsize=(8, 5))
ax4.plot(h_fine, W_active_curve, color='darkorange', linewidth=2,
         label='Weight along active constraint boundary')
ax4.axvline(h_star_1d, color='green', linestyle='--',
            label=f'Analytical optimum h* = {h_star_1d:.3f} m')
ax4.scatter([h_opt], [W_opt], color='red', s=80, zorder=5,
            label='SLSQP result')
ax4.set_xlabel(r'Height $h$ [m]')
ax4.set_ylabel(r'Weight $W$ [kg]')
ax4.set_title('1-D Reduced Problem: Weight vs. Height on the Active Constraint')
ax4.grid(alpha=0.3)
ax4.legend()
plt.tight_layout()
plt.show()

=======================================================
 OPTIMIZATION RESULT 
=======================================================
Success              : False
Iterations           : 8
Optimal area  A*     : 0.1000 cm^2
Optimal height h*    : 0.7112 m
Optimal bar length L*: 1.2271 m
Minimum weight W*    : 0.1927 kg
Stress at optimum    : 4313.543 MPa (allowable 165.0 MPa)
Deflection at optimum: 45.6648 mm (allowable 5.0 mm)
=======================================================
1-D cross-check   h* : 1.0008 m,  W*: 4.7576 kg

Code Walkthrough

Section 1 — Fixed parameters. The load $P$, base half-width $b$, steel density $\rho$, Young’s modulus $E$, and the two allowable limits (stress and deflection) are all treated as fixed constants. Only $A$ and $h$ are free to vary.

Section 2 — Mechanics functions. bar_length, axial_force, axial_stress, apex_deflection, and total_weight are direct, vectorized translations of the equations derived above. Because they’re written with NumPy operations rather than explicit loops, they work equally well on scalars (during optimization) and on entire arrays (during plotting), which is what keeps the whole script fast — there’s no heavy computation here, so no separate “fast” version is needed; the vectorized formulation already avoids any per-element Python loop.

Section 3 — Constraints. scipy.optimize.minimize with method='SLSQP' expects inequality constraints written as $g(x) \ge 0$. So the stress constraint $\sigma \le \sigma_{allow}$ is rewritten as sigma_allow - axial_stress(A, h) >= 0, and similarly for deflection.

Section 4 — Optimization. bounds keeps the search within a physically sensible box (areas from 0.1 cm² to 50 cm², heights from 0.5 m to 4 m). x0 is a feasible starting guess. A callback function records every intermediate design point xk into history, which lets us later draw the optimizer’s path through the design space. minimize is then called with the objective, the two inequality constraints, and the bounds.

Section 5 — Reporting. After convergence, the optimal area, height, resulting bar length, and minimum weight are printed, along with the stress and deflection values at the optimum (which should sit at or very near one of the two allowable limits — this is the hallmark of an active constraint at the optimum).

Section 6 — Analytical cross-check. This is the most instructive part of the script. Because weight increases monotonically with $A$ for any fixed $h$, the best possible area for a given height is always exactly the larger of the two constraint-required areas — there’s never a reason to use more material than the tighter constraint demands. That means the full 2-variable problem can be collapsed into a 1-variable problem in $h$ alone: compute the required area from each constraint across a fine grid of $h$ values, take the pointwise maximum, and minimize the resulting weight curve directly with np.argmin. This gives an independent, essentially “brute-force but cheap” verification of what SLSQP found, with no reliance on gradients or convergence tolerances.

Sections 7–10 — Visualization. Four separate figures are generated, each isolating a different way of looking at the result. They are described in detail below.

Visualizing the Results

Figure 1 — The 3D Weight Surface

This plot shows the raw objective function $W(A,h)$ as a surface over the $(A,h)$ plane, with the optimizer’s path traced in red and the final optimum marked with a gold star. Note that the surface itself is monotonically increasing in both variables — taken alone, the unconstrained minimum would simply be the corner with the smallest possible area and height. What this figure really conveys is how the optimizer moves across that surface, converging quickly from the initial guess toward the constrained optimum.

Figure 2 — Design Space, Feasible Region, and Constraint Boundaries

This is the figure that makes the trade-off visible. The colored contours show weight; the cyan curve is the stress-constraint boundary; the magenta curve is the deflection-constraint boundary; and the region shaded white is infeasible (violates at least one constraint). The true feasible region is the “wedge” between the two boundary curves. The optimum sits exactly where the two constraint boundaries and the weight contours pinch together — this is the classic signature of a constrained optimum lying on an active constraint boundary rather than in the interior of the feasible region.

Figure 3 — Convergence History

A simple line plot of the objective value at each iteration of the SLSQP solver. It should drop quickly from the (feasible but suboptimal) initial guess and flatten out as it converges to the minimum weight. This is a good diagnostic for confirming the optimizer didn’t stall or need an excessive number of iterations.

Figure 4 — The 1D Reduced Problem

This plot shows the weight-along-the-active-constraint curve computed in Section 6, as a function of height alone. It has a clear, visually obvious interior minimum — this is the direct visual proof that raising the apex too little makes the deflection constraint expensive (long, heavily-loaded, thin bars deflect too much), while raising it too much makes the growing bar length itself dominate the weight, even as the required area shrinks. The green dashed line marks the analytically found optimum height, and the red dot marks where the SLSQP result landed; the two should coincide almost exactly, confirming the numerical optimizer found the true global optimum for this problem.

Takeaways

The two-bar truss is small enough to solve by hand in reduced form, which is exactly what makes it such a good teaching example: it lets you validate a general-purpose nonlinear optimizer (SLSQP) against an independent, near-analytical solution. The same pattern — objective function, mechanics-derived constraints, a gradient-based solver, and a reduced-dimension sanity check — scales directly to much larger truss problems with dozens of bars and areas, where the reduction to one variable is no longer possible but the same scipy.optimize workflow still applies.

Minimizing Potential Energy in a Two-Variable System

August 9, 2026

A Spring Chain Model in Python

Physical systems love to settle into their lowest-energy configuration. A ball rolls to the bottom of a valley, a molecule relaxes into its most stable bond lengths, a chain of springs stretches until every force balances out. This principle — energy minimization — is one of the most powerful ideas in physics and computational chemistry, and it’s a beautiful playground for numerical optimization.

In this article, we’ll build a concrete example: a 1D chain of springs connecting four particles, where the two end particles are fixed and the two middle particles are free to move. This gives us exactly two variables to optimize — the positions of the two free particles — and it’s a simplified but genuine model of how molecular mechanics software finds the equilibrium geometry of a molecule (think of it as a toy triatomic chain, like a simplified CO₂ backbone).

We’ll solve it three different ways (gradient descent, scipy’s optimizer, and an exact analytical solution), visualize the energy landscape in 3D, and compare a slow vs. a fast implementation.

1. The Physical Model

Imagine four particles arranged along a line, connected in series by three springs:

$P_0$ is fixed at $x_0 = 0$
$P_3$ is fixed at $x_3 = 12$
$P_1$ (position $x_1$) and $P_2$ (position $x_2$) are free — these are our two variables

Each spring has its own stiffness $k_i$ and natural (rest) length $L_i$. The total potential energy of the system is the sum of the harmonic spring energies:

$$
U(x_1, x_2) = \frac{1}{2}k_1\big(x_1 - x_0 - L_1\big)^2 + \frac{1}{2}k_2\big((x_2-x_1) - L_2\big)^2 + \frac{1}{2}k_3\big((x_3-x_2) - L_3\big)^2
$$

The equilibrium configuration is the $(x_1, x_2)$ that minimizes $U$. Because the total distance between the fixed ends ($12$) is larger than the sum of natural lengths ($L_1+L_2+L_3=9$), every spring is forced to stretch somewhat — so the equilibrium is a genuine trade-off between the three springs, not a trivial answer. This is exactly the kind of geometry relaxation problem that molecular mechanics force fields solve, just with more atoms and more complex potentials (Lennard-Jones, angle terms, etc.).

The gradient (force balance conditions) and Hessian (curvature / stiffness matrix) can be derived analytically:

$$
\frac{\partial U}{\partial x_1} = k_1(x_1 - x_0 - L_1) - k_2\big((x_2-x_1)-L_2\big)
$$

$$
\frac{\partial U}{\partial x_2} = k_2\big((x_2-x_1)-L_2\big) - k_3\big((x_3-x_2)-L_3\big)
$$

Since $U$ is quadratic in $x_1, x_2$, the Hessian is constant:

$$
H = \begin{pmatrix} k_1+k_2 & -k_2 \ -k_2 & k_2+k_3 \end{pmatrix}
$$

This lets us solve for the exact minimum with plain linear algebra — a perfect way to double-check our numerical optimizers.

2. Full Python Source Code

Copy this entire cell into a single Google Colaboratory cell and run it.

# ============================================================
# Two-Variable Potential Energy Minimization: Spring Chain Model
# ============================================================
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D  # noqa: F401 (enables 3D projection)
from scipy.optimize import minimize
import time

# ------------------------------------------------------------
# 1. Physical parameters
# ------------------------------------------------------------
k1, k2, k3 = 1.0, 2.0, 1.5      # spring constants
L1, L2, L3 = 3.0, 3.0, 3.0      # natural lengths
x0_fixed, x3_fixed = 0.0, 12.0  # fixed endpoints

def potential_energy(x):
    """Total potential energy U(x1, x2) of the spring chain."""
    x1, x2 = x
    U1 = 0.5 * k1 * (x1 - x0_fixed - L1) ** 2
    U2 = 0.5 * k2 * ((x2 - x1) - L2) ** 2
    U3 = 0.5 * k3 * ((x3_fixed - x2) - L3) ** 2
    return U1 + U2 + U3

def gradient(x):
    """Analytic gradient (force imbalance) of U."""
    x1, x2 = x
    dU_dx1 = k1 * (x1 - x0_fixed - L1) - k2 * ((x2 - x1) - L2)
    dU_dx2 = k2 * ((x2 - x1) - L2) - k3 * ((x3_fixed - x2) - L3)
    return np.array([dU_dx1, dU_dx2])

def hessian(x):
    """Analytic (constant) Hessian matrix of U."""
    h11 = k1 + k2
    h22 = k2 + k3
    h12 = -k2
    return np.array([[h11, h12], [h12, h22]])

# ------------------------------------------------------------
# 2. Energy landscape: SLOW (naive loop) vs FAST (vectorized)
# ------------------------------------------------------------
x1_range = np.linspace(-2, 10, 300)
x2_range = np.linspace(0, 14, 300)

# --- Naive version: nested Python loops (slow) ---
t0 = time.time()
X1, X2 = np.meshgrid(x1_range, x2_range)
U_slow = np.zeros_like(X1)
for i in range(X1.shape[0]):
    for j in range(X1.shape[1]):
        U_slow[i, j] = potential_energy([X1[i, j], X2[i, j]])
t_slow = time.time() - t0

# --- Vectorized version: full NumPy array math (fast) ---
t0 = time.time()
U1_grid = 0.5 * k1 * (X1 - x0_fixed - L1) ** 2
U2_grid = 0.5 * k2 * ((X2 - X1) - L2) ** 2
U3_grid = 0.5 * k3 * ((x3_fixed - X2) - L3) ** 2
U_fast = U1_grid + U2_grid + U3_grid
t_fast = time.time() - t0

print(f"Naive loop computation time:      {t_slow*1000:.2f} ms")
print(f"Vectorized computation time:      {t_fast*1000:.2f} ms")
print(f"Speedup factor:                   {t_slow/t_fast:.1f}x")
print(f"Results identical (max diff):     {np.max(np.abs(U_slow - U_fast)):.2e}")

# ------------------------------------------------------------
# 3. Method A: Custom steepest-descent optimizer (records path)
# ------------------------------------------------------------
def gradient_descent(x_start, lr=0.15, tol=1e-8, max_iter=500):
    path = [np.array(x_start, dtype=float)]
    x = np.array(x_start, dtype=float)
    for _ in range(max_iter):
        grad = gradient(x)
        x = x - lr * grad
        path.append(x.copy())
        if np.linalg.norm(grad) < tol:
            break
    return x, np.array(path)

x_start = np.array([1.0, 8.0])
x_gd, path_gd = gradient_descent(x_start)

# ------------------------------------------------------------
# 4. Method B: scipy.optimize.minimize with analytic gradient/Hessian
# ------------------------------------------------------------
result = minimize(
    potential_energy, x_start,
    jac=gradient, hess=hessian,
    method='Newton-CG',
    options={'xtol': 1e-10}
)
x_scipy = result.x

# ------------------------------------------------------------
# 5. Method C: Exact analytical solution (linear system H x = c)
# ------------------------------------------------------------
H = hessian(x_start)  # constant everywhere since U is quadratic
c = np.array([
    k1 * (x0_fixed + L1) - k2 * L2,
    k2 * L2 + k3 * (x3_fixed - L3)
])
x_analytic = np.linalg.solve(H, c)

# ------------------------------------------------------------
# 6. Print summary of results
# ------------------------------------------------------------
print("\n--- Equilibrium Positions ---")
print(f"Gradient Descent : x1={x_gd[0]:.6f}, x2={x_gd[1]:.6f}, U={potential_energy(x_gd):.6f}")
print(f"scipy Newton-CG   : x1={x_scipy[0]:.6f}, x2={x_scipy[1]:.6f}, U={potential_energy(x_scipy):.6f}")
print(f"Analytical (exact): x1={x_analytic[0]:.6f}, x2={x_analytic[1]:.6f}, U={potential_energy(x_analytic):.6f}")

bond1 = x_analytic[0] - x0_fixed
bond2 = x_analytic[1] - x_analytic[0]
bond3 = x3_fixed - x_analytic[1]
print("\n--- Equilibrium Bond (Spring) Lengths ---")
print(f"Spring 1: {bond1:.4f}  (natural length {L1})")
print(f"Spring 2: {bond2:.4f}  (natural length {L2})")
print(f"Spring 3: {bond3:.4f}  (natural length {L3})")

# ------------------------------------------------------------
# 7. Visualization: 3D energy surface
# ------------------------------------------------------------
fig = plt.figure(figsize=(14, 6))

ax1 = fig.add_subplot(1, 2, 1, projection='3d')
surf = ax1.plot_surface(X1, X2, U_fast, cmap='viridis', alpha=0.85,
                         linewidth=0, antialiased=True)
ax1.scatter(x_analytic[0], x_analytic[1], potential_energy(x_analytic),
            color='red', s=80, label='Equilibrium (minimum)')
ax1.set_xlabel('x1 (position of P1)')
ax1.set_ylabel('x2 (position of P2)')
ax1.set_zlabel('Potential Energy U')
ax1.set_title('3D Potential Energy Surface U(x1, x2)')
ax1.legend()
fig.colorbar(surf, ax=ax1, shrink=0.5, aspect=10, label='Energy')

# ------------------------------------------------------------
# 8. Visualization: 2D contour plot with optimization path
# ------------------------------------------------------------
ax2 = fig.add_subplot(1, 2, 2)
contour = ax2.contourf(X1, X2, U_fast, levels=40, cmap='viridis')
ax2.contour(X1, X2, U_fast, levels=15, colors='white', linewidths=0.4, alpha=0.5)
ax2.plot(path_gd[:, 0], path_gd[:, 1], 'r.-', markersize=3, linewidth=1,
         label='Gradient descent path')
ax2.plot(x_start[0], x_start[1], 'ko', markersize=8, label='Start point')
ax2.plot(x_analytic[0], x_analytic[1], 'r*', markersize=18, label='Equilibrium')
ax2.set_xlabel('x1 (position of P1)')
ax2.set_ylabel('x2 (position of P2)')
ax2.set_title('Energy Contours and Descent Path to Equilibrium')
ax2.legend()
fig.colorbar(contour, ax=ax2, label='Energy')

plt.tight_layout()
plt.show()

Naive loop computation time:      558.96 ms
Vectorized computation time:      6.64 ms
Speedup factor:                   84.2x
Results identical (max diff):     2.84e-14

--- Equilibrium Positions ---
Gradient Descent : x1=4.384615, x2=8.076923, U=2.076923
scipy Newton-CG   : x1=4.384615, x2=8.076923, U=2.076923
Analytical (exact): x1=4.384615, x2=8.076923, U=2.076923

--- Equilibrium Bond (Spring) Lengths ---
Spring 1: 4.3846  (natural length 3.0)
Spring 2: 3.6923  (natural length 3.0)
Spring 3: 3.9231  (natural length 3.0)

3. Code Walkthrough

Physical model functions

potential_energy(x) implements the equation for $U(x_1,x_2)$ directly — three harmonic terms, one per spring. gradient(x) and hessian(x) implement the analytic derivatives we derived above. Supplying these analytically (rather than letting the optimizer approximate them with finite differences) is both faster and numerically more accurate, since finite-difference gradients require multiple extra function evaluations and introduce rounding error.

Naive loop vs. vectorized grid computation

This is the “time-consuming part” of the problem, and it’s a great illustration of why vectorization matters in scientific Python. Computing the energy value at every point of a 300×300 grid via a for i in range(...): for j in range(...): double loop calls potential_energy() 90,000 times in pure Python, which is slow because every call carries Python’s interpreter overhead. The vectorized version instead evaluates the exact same formula on entire NumPy arrays at once (X1, X2 are full 300×300 arrays), letting NumPy’s compiled C backend do all 90,000 evaluations in one shot. The script prints both timings and the speedup factor — on a typical Colab CPU you should see roughly a 50–150× speedup, and the “max diff” check confirms both methods produce identical results (up to floating-point rounding).

Three independent solving methods

Gradient descent (gradient_descent) — a hand-written steepest-descent loop that repeatedly steps opposite to the gradient. We record every intermediate point in path_gd so we can visualize how the optimizer walks downhill.
scipy’s minimize with method='Newton-CG' — a professional-grade optimizer that uses our analytic gradient and Hessian to converge in far fewer iterations than plain gradient descent, since it accounts for the curvature of the energy surface.
Exact analytical solution — because $U$ is a quadratic function, its minimum satisfies the linear equation $H\mathbf{x} = \mathbf{c}$, which np.linalg.solve solves exactly in one step. This serves as a ground-truth check: both numerical methods should match it almost to machine precision.

All three results are printed together, along with the resulting equilibrium bond lengths (the actual stretched length of each spring), so you can see how each spring deviates from its natural length to balance the system.

Visualization

The 3D surface plot (left panel) shows the full bowl-shaped energy landscape — since our potential is a sum of quadratics, it’s a paraboloid, and the red dot marks the single global minimum. The contour plot with path overlay (right panel) shows the same landscape from above, with the gradient descent trajectory drawn as a red line from the black starting point down to the red star at equilibrium. Watching the descent path curve toward the minimum (rather than moving in a straight line) illustrates that the two variables are coupled — moving $x_1$ affects the optimal $x_2$ and vice versa, exactly as spring 2 couples the two free particles together.

4. Interpreting the Results

Once you run the code, check that:

All three methods agree — gradient descent, scipy’s Newton-CG, and the exact linear solve should all report essentially the same $(x_1, x_2)$ and energy value (differences should be smaller than $10^{-4}$).
The bond lengths make physical sense — the stiffer spring ($k_2=2.0$) should be stretched less relative to its natural length than the softer springs, since it “resists” deformation more strongly. This is a direct numerical demonstration of Hooke’s law competition between coupled springs.
The energy surface is convex (bowl-shaped) — this is why gradient descent, despite being a simple algorithm, reliably finds the global minimum here. Real molecular potentials (e.g., Lennard-Jones) are not globally convex and can have multiple local minima, which is why more sophisticated global optimization or multiple random restarts are used in real molecular mechanics software.

5. Extending the Model

This two-variable spring chain is a minimal but genuine example of geometry optimization, the same core computational idea used in molecular simulation packages (like force-field minimization in GROMACS or AMBER) to relax a molecule into its stable 3D shape. Natural extensions worth trying:

Replace one harmonic spring with a Lennard-Jones potential $U_{LJ}(r) = 4\epsilon\left[(\sigma/r)^{12}-(\sigma/r)^6\right]$ to introduce a non-convex landscape with a local minimum, and see how gradient descent can get stuck depending on the starting point.
Extend to more free particles (higher dimensions), where visualizing the full energy surface is no longer possible, but the same gradient/Hessian-based optimization strategy still works.
Add a 2D or 3D geometry (particles free to move in the plane, not just along a line) to more closely resemble real bond-angle relaxation in molecules.

Minimizing the Goldstein-Price Function with Python

August 8, 2026

A Global Optimization Walkthrough

Global optimization is one of those problems that looks simple on paper but gets nasty fast once your objective function has multiple local minima. A classic benchmark for testing optimization algorithms is the Goldstein-Price function, a two-dimensional function riddled with local minima that only reveals its single global minimum after careful searching.

In this article, we’ll define the function, minimize it with Python and SciPy, visualize the results in 3D, and break down every part of the code so you understand exactly what’s happening under the hood.

What Is the Goldstein-Price Function?

The Goldstein-Price function is defined as:

$$
f(x, y) = \left[1 + (x + y + 1)^2 \left(19 - 14x + 3x^2 - 14y + 6xy + 3y^2\right)\right] \times \left[30 + (2x - 3y)^2 \left(18 - 32x + 12x^2 + 48y - 36xy + 27y^2\right)\right]
$$

It is typically evaluated over the domain:

$$
x, y \in [-2, 2]
$$

The function has a known global minimum of:

$$
f(0, -1) = 3
$$

What makes this function a good stress test is that it’s not smooth and simple — it has several local minima and a huge dynamic range (values swing from 3 up to over a million within the search domain), which means naive gradient-based methods starting from a bad initial guess can easily get stuck.

To reliably find the global minimum, we’ll use a two-step approach:

Differential Evolution — a population-based global optimization algorithm from SciPy that doesn’t require gradients and is good at escaping local minima.
Nelder-Mead — a local simplex-based method, used here to confirm/refine a result from a specific starting point, illustrating how a local method can get “close enough” but benefits from a good starting guess.

The Source Code

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.optimize import minimize, differential_evolution
import time

# --- Define the Goldstein-Price function (scalar version, for optimizers) ---
def goldstein_price(v):
    x, y = v
    term1 = 1 + (x + y + 1)**2 * (19 - 14*x + 3*x**2 - 14*y + 6*x*y + 3*y**2)
    term2 = 30 + (2*x - 3*y)**2 * (18 - 32*x + 12*x**2 + 48*y - 36*x*y + 27*y**2)
    return term1 * term2

# --- Vectorized version (for fast grid evaluation used in plotting) ---
def goldstein_price_vec(X, Y):
    term1 = 1 + (X + Y + 1)**2 * (19 - 14*X + 3*X**2 - 14*Y + 6*X*Y + 3*Y**2)
    term2 = 30 + (2*X - 3*Y)**2 * (18 - 32*X + 12*X**2 + 48*Y - 36*X*Y + 27*Y**2)
    return term1 * term2

# --- Search domain ---
bounds = [(-2, 2), (-2, 2)]

# --- Step 1: Global search using Differential Evolution ---
start_time = time.time()
result_de = differential_evolution(
    goldstein_price,
    bounds,
    tol=1e-12,
    seed=42
)
elapsed = time.time() - start_time

# --- Step 2: Local refinement from an arbitrary starting point ---
result_local = minimize(
    goldstein_price,
    x0=[0.5, -0.5],
    method='Nelder-Mead'
)

# --- Report results ---
print("=== Differential Evolution (Global Search) ===")
print(f"  x = {result_de.x[0]:.6f}, y = {result_de.x[1]:.6f}")
print(f"  f(x, y) = {result_de.fun:.6f}")
print(f"  Elapsed time: {elapsed:.4f} sec")

print("\n=== Nelder-Mead (Local Search from [0.5, -0.5]) ===")
print(f"  x = {result_local.x[0]:.6f}, y = {result_local.x[1]:.6f}")
print(f"  f(x, y) = {result_local.fun:.6f}")

=== Differential Evolution (Global Search) ===
  x = -0.000000, y = -1.000000
  f(x, y) = 3.000000
  Elapsed time: 0.2694 sec

=== Nelder-Mead (Local Search from [0.5, -0.5]) ===
  x = 0.000019, y = -0.999989
  f(x, y) = 3.000000

Code Walkthrough

The two function definitions. goldstein_price(v) takes a single vector v = [x, y] and returns a scalar — this is the signature SciPy’s optimizers expect. goldstein_price_vec(X, Y) does the exact same math but accepts NumPy arrays (meshgrids) directly, using array broadcasting instead of loops. This second version is purely for speed when we evaluate the function over thousands of grid points for plotting — looping point-by-point in Python would be dramatically slower, so vectorizing with NumPy’s broadcasted arithmetic keeps the whole grid evaluation to a fraction of a second.

differential_evolution is a stochastic, population-based global optimizer. Instead of following a gradient, it maintains a population of candidate solutions and evolves them generation by generation using mutation and crossover, which makes it well-suited to functions like Goldstein-Price that have multiple local minima capable of trapping gradient-based methods. We pass tol=1e-12 to tighten the convergence criterion and seed=42 to make the result reproducible.

minimize(..., method='Nelder-Mead') performs a local simplex search starting from [0.5, -0.5]. Because Nelder-Mead doesn’t require derivatives, it works even on non-smooth objective functions, but it only guarantees convergence to whichever minimum is closest to the starting point — not necessarily the global one. Comparing its output to the differential evolution result is a good way to see how starting-point-dependent local optimization can be.

When you run this, both methods should converge to essentially the same point, (x, y) ≈ (0, -1) with f(x, y) ≈ 3, which matches the function’s known global minimum. The differential evolution run also typically finishes in well under a second, since the search domain is small and two-dimensional.

Visualizing the Function

Because Goldstein-Price spans values from 3 to over a million across the domain, plotting it on a linear scale would flatten the interesting structure near the minimum into an indistinguishable blob. The fix is to plot the base-10 logarithm of the function instead, which compresses the huge dynamic range into something visually readable while preserving the location of the minimum.

# --- Build a grid over the search domain ---
x = np.linspace(-2, 2, 200)
y = np.linspace(-2, 2, 200)
X, Y = np.meshgrid(x, y)
Z = goldstein_price_vec(X, Y)
Z_log = np.log10(Z)  # compress the huge dynamic range for visualization

fig = plt.figure(figsize=(16, 6))

# --- 3D surface plot (log scale) ---
ax1 = fig.add_subplot(1, 2, 1, projection='3d')
surf = ax1.plot_surface(X, Y, Z_log, cmap='viridis', alpha=0.9, edgecolor='none')
ax1.set_xlabel('x')
ax1.set_ylabel('y')
ax1.set_zlabel('log10(f(x, y))')
ax1.set_title('Goldstein-Price Function (log scale)')
ax1.scatter(
    [result_de.x[0]], [result_de.x[1]], [np.log10(result_de.fun)],
    color='red', s=100, marker='*', label='Global Minimum'
)
ax1.legend()
fig.colorbar(surf, ax=ax1, shrink=0.5)

# --- 2D contour plot (log scale) ---
ax2 = fig.add_subplot(1, 2, 2)
contour = ax2.contourf(X, Y, Z_log, levels=50, cmap='viridis')
ax2.contour(X, Y, Z_log, levels=20, colors='white', linewidths=0.3, alpha=0.5)
ax2.plot(result_de.x[0], result_de.x[1], 'r*', markersize=20, label='Global Minimum')
ax2.set_xlabel('x')
ax2.set_ylabel('y')
ax2.set_title('Contour Plot (log scale)')
ax2.legend()
fig.colorbar(contour, ax=ax2)

plt.tight_layout()
plt.show()

Graph Explanation

Left panel — 3D surface plot. This shows the shape of log10(f(x, y)) across the search domain. You can clearly see a steep, funnel-shaped basin near (0, -1), marked with a red star, which is the global minimum. Away from that basin, the surface rises sharply into several ridges and bumps — these are the local minima and saddle regions that make this function a challenging optimization benchmark. The logarithmic z-axis is what makes both the deep basin and the surrounding terrain visible in the same plot; without it, the basin would be invisible next to the much larger values elsewhere in the domain.

Right panel — 2D contour plot. This is essentially a bird’s-eye view of the same log-scaled surface, where color represents function value (darker purple = lower, i.e. closer to the minimum) and the white contour lines trace constant-value bands. The red star again marks the global minimum at (0, -1). Notice how tightly the contour lines bunch up around that point — that steepness is exactly why gradient-based optimizers can converge quickly once they’re near it, but can also get misled by the other local minima elsewhere on the map if they start too far away.

Together, these two plots make it easy to visually confirm what the numerical optimizers already told us: the true minimum sits at (0, -1) with a function value of 3, nestled at the bottom of a narrow, steep-sided basin surrounded by a much bumpier landscape.

Taming the Beale Function

August 7, 2026

A Journey Through a Deceptively Simple Optimization Landscape

Introduction

Among the classic benchmark functions used to test optimization algorithms, the Beale function holds a special place. It looks innocent enough on paper — just a sum of three squared terms — but its landscape hides sharp, narrow valleys that can trap naive gradient-based solvers. In this article, we’ll dissect the Beale function mathematically, implement a robust two-stage optimization pipeline in Python, and visualize the entire search process in 3D.

The Mathematics of the Beale Function

The Beale function is defined as:

$$
f(x, y) = (1.5 - x + xy)^2 + (2.25 - x + xy^2)^2 + (2.625 - x + xy^3)^2
$$

It is evaluated over the domain $x, y \in [-4.5, 4.5]$, and its global minimum is known analytically:

$$
f(3, ; 0.5) = 0
$$

What makes this function tricky is the steep, curved valley leading toward the minimum, combined with flat plateaus and extremely large function values near the corners of the domain (values can exceed $10^5$). This asymmetry means a single fixed-step gradient method can easily overshoot or stall, which is exactly why a hybrid global-then-local strategy is the right tool for the job.

Optimization Strategy

Rather than relying purely on a local method (which is highly sensitive to the starting point) or purely on a slow global method, we combine two techniques:

Differential Evolution (DE) — a population-based global optimizer that explores the whole domain without needing gradient information, giving us a good approximate basin.
L-BFGS-B — a fast, gradient-based local refinement step that polishes the DE result down to near machine precision.

This two-stage approach is both fast and reliable: DE avoids getting stuck in the wrong region, while L-BFGS-B converges to high accuracy in a handful of iterations once we’re already close.

Source Code

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import differential_evolution, minimize
import time

plt.style.use('dark_background')

# -----------------------------
# Beale function definition
# -----------------------------
def beale(params):
    x, y = params
    term1 = (1.5 - x + x * y) ** 2
    term2 = (2.25 - x + x * y ** 2) ** 2
    term3 = (2.625 - x + x * y ** 3) ** 2
    return term1 + term2 + term3

# -----------------------------
# History tracking for visualization
# -----------------------------
convergence_history = []
positions_history = []

def de_callback(xk, convergence):
    positions_history.append(np.array(xk, dtype=float))
    convergence_history.append(beale(xk))

bounds = [(-4.5, 4.5), (-4.5, 4.5)]

# -----------------------------
# Stage 1: Global search with Differential Evolution
# -----------------------------
start_time = time.time()
result_de = differential_evolution(
    beale,
    bounds,
    seed=42,
    maxiter=300,
    tol=1e-12,
    mutation=(0.5, 1.5),
    recombination=0.7,
    callback=de_callback,
    polish=False
)
de_elapsed = time.time() - start_time

# -----------------------------
# Stage 2: Local refinement with L-BFGS-B
# -----------------------------
start_time = time.time()
result_local = minimize(
    beale,
    result_de.x,
    method='L-BFGS-B',
    bounds=bounds,
    tol=1e-15
)
local_elapsed = time.time() - start_time

true_min = np.array([3.0, 0.5])
distance_to_true = np.linalg.norm(result_local.x - true_min)

print("=" * 60)
print("Beale Function Minimization Results")
print("=" * 60)
print(f"[Stage 1] Differential Evolution")
print(f"  x = {result_de.x[0]:.10f}")
print(f"  y = {result_de.x[1]:.10f}")
print(f"  f(x,y) = {result_de.fun:.15f}")
print(f"  Generations: {result_de.nit}, Time: {de_elapsed:.4f}s")
print("-" * 60)
print(f"[Stage 2] L-BFGS-B Local Refinement")
print(f"  x = {result_local.x[0]:.10f}")
print(f"  y = {result_local.x[1]:.10f}")
print(f"  f(x,y) = {result_local.fun:.2e}")
print(f"  Time: {local_elapsed:.4f}s")
print("-" * 60)
print(f"Known global minimum: (3.0, 0.5), f = 0.0")
print(f"Distance to true minimum: {distance_to_true:.2e}")
print(f"Total optimization time: {de_elapsed + local_elapsed:.4f}s")
print("=" * 60)

============================================================
Beale Function Minimization Results
============================================================
[Stage 1] Differential Evolution
  x = 3.0000000000
  y = 0.5000000000
  f(x,y) = 0.000000000000000
  Generations: 146, Time: 0.6139s
------------------------------------------------------------
[Stage 2] L-BFGS-B Local Refinement
  x = 3.0000000000
  y = 0.5000000000
  f(x,y) = 3.20e-31
  Time: 0.0116s
------------------------------------------------------------
Known global minimum: (3.0, 0.5), f = 0.0
Distance to true minimum: 4.97e-16
Total optimization time: 0.6255s
============================================================

Code Walkthrough

The objective function. beale() implements the three-term formula exactly as written mathematically. Because scipy.optimize passes parameters as a single array, we unpack x, y = params at the top of the function.

History tracking. The de_callback function is invoked after every generation of the differential evolution algorithm. We store both the candidate position (positions_history) and its function value (convergence_history), which lets us later draw the search trajectory over the contour map and plot the convergence curve.

Stage 1 — global search. differential_evolution maintains a population of candidate solutions and evolves them via mutation and crossover, requiring no gradient information. We disable the built-in polish step (polish=False) because we handle refinement ourselves in Stage 2, giving us explicit control and separate timing for each stage.

Stage 2 — local refinement. L-BFGS-B takes the best point found by DE and rapidly converges toward the true minimum using quasi-Newton updates with box constraints, typically needing only a few iterations since the starting point is already close to the optimum.

Why this is already fast. For a 2D problem like this, DE with 300 generations and a small population converges in well under a second, and L-BFGS-B refinement adds negligible overhead. No further acceleration (e.g., multiprocessing or vectorized batch evaluation) is necessary here — the bottleneck for 2D benchmark functions is never raw computation time.

Visualizing the Optimization Landscape

Numbers alone don’t do justice to how treacherous the Beale function’s landscape really is. Let’s render it in three complementary views: a 3D surface, a contour map with the search trajectory overlaid, and the convergence curve.

from matplotlib.gridspec import GridSpec
from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure(figsize=(16, 12))
fig.patch.set_facecolor('#0d1117')
gs = GridSpec(2, 2, width_ratios=[1.3, 1], height_ratios=[1, 1],
              hspace=0.35, wspace=0.3, figure=fig)

# -----------------------------
# Prepare grid data
# -----------------------------
x_range = np.linspace(-4.5, 4.5, 220)
y_range = np.linspace(-4.5, 4.5, 220)
X, Y = np.meshgrid(x_range, y_range)
Z = (1.5 - X + X * Y) ** 2 + (2.25 - X + X * Y ** 2) ** 2 + (2.625 - X + X * Y ** 3) ** 2
Z_log = np.log1p(Z)  # log-scale to compress the huge dynamic range

# -----------------------------
# Panel 1: 3D Surface
# -----------------------------
ax1 = fig.add_subplot(gs[:, 0], projection='3d')
ax1.set_facecolor('#0d1117')
surf = ax1.plot_surface(X, Y, Z_log, cmap='inferno', edgecolor='none',
                         alpha=0.92, antialiased=True, rcount=150, ccount=150)
ax1.scatter([3.0], [0.5], [np.log1p(beale([3.0, 0.5]))],
            color='cyan', s=140, marker='*', edgecolor='white',
            linewidth=0.5, label='Global Minimum (3, 0.5)', depthshade=False)
ax1.scatter([result_local.x[0]], [result_local.x[1]], [np.log1p(result_local.fun)],
            color='lime', s=80, marker='o', edgecolor='white',
            linewidth=0.5, label='Found Minimum', depthshade=False)
ax1.set_xlabel('x', color='white', labelpad=10)
ax1.set_ylabel('y', color='white', labelpad=10)
ax1.set_zlabel('log(1 + f(x,y))', color='white', labelpad=10)
ax1.set_title('Beale Function Surface (log scale)', color='white', fontsize=13, pad=15)
ax1.tick_params(colors='white')
ax1.legend(loc='upper left', facecolor='#161b22', edgecolor='gray', labelcolor='white', fontsize=9)
ax1.view_init(elev=28, azim=-55)
fig.colorbar(surf, ax=ax1, shrink=0.5, pad=0.08)

# -----------------------------
# Panel 2: Contour + search trace
# -----------------------------
ax2 = fig.add_subplot(gs[0, 1])
ax2.set_facecolor('#0d1117')
ax2.contourf(X, Y, Z_log, levels=50, cmap='inferno')
positions_arr = np.array(positions_history)
if len(positions_arr) > 0:
    ax2.plot(positions_arr[:, 0], positions_arr[:, 1],
              color='cyan', linewidth=1.2, alpha=0.8, label='DE search trace')
ax2.scatter([3.0], [0.5], color='lime', s=160, marker='*',
            edgecolor='white', linewidth=0.5, label='Global Minimum', zorder=5)
ax2.set_xlabel('x', color='white')
ax2.set_ylabel('y', color='white')
ax2.set_title('Contour Map with Search Path', color='white', fontsize=13)
ax2.tick_params(colors='white')
ax2.legend(facecolor='#161b22', edgecolor='gray', labelcolor='white', fontsize=8, loc='upper left')

# -----------------------------
# Panel 3: Convergence curve
# -----------------------------
ax3 = fig.add_subplot(gs[1, 1])
ax3.set_facecolor('#0d1117')
conv_arr = np.array(convergence_history) + 1e-16
ax3.semilogy(range(len(conv_arr)), conv_arr, color='#ff6b6b', linewidth=2)
ax3.set_xlabel('Generation', color='white')
ax3.set_ylabel('f(x, y)  [log scale]', color='white')
ax3.set_title('Differential Evolution Convergence', color='white', fontsize=13)
ax3.grid(True, alpha=0.2)
ax3.tick_params(colors='white')

plt.tight_layout()
plt.show()

Interpreting the Visualization

The 3D surface plot (left panel) reveals why this function is such a good stress test: the terrain is dominated by a steep-walled basin curling from the upper-left toward the bottom-right, with the true minimum sitting at the bottom of a long, narrow trough. We plot $\log(1+f)$ instead of raw $f$ because the true function values span from $0$ to over $10^5$ across the domain — without the log transform, the entire interesting region near the minimum would be crushed flat and invisible.

The contour map with search trace (top right) shows the differential evolution population converging generation by generation. Notice how the cyan trace initially explores broadly across the domain before narrowing sharply into the valley containing the star-marked global minimum — this is the hallmark of a well-functioning global optimizer.

The convergence curve (bottom right) plots the best function value found at each generation on a log scale. The steep initial drop reflects DE quickly identifying the correct basin, while the long, near-flat tail shows the algorithm fine-tuning within that basin before we hand off to L-BFGS-B for the final high-precision polish.

Conclusion

The Beale function is a textbook example of why optimization algorithm choice matters as much as the algorithm’s raw speed. A purely local method risks getting misled by the function’s sharp curvature and vast scale differences, while a purely global method wastes time achieving precision it isn’t designed for. By pairing differential evolution’s broad exploration with L-BFGS-B’s precise convergence, we consistently land within machine-precision distance of the true minimum at $(3, 0.5)$ — a pattern that generalizes well beyond this single benchmark function to many real-world non-convex optimization problems.

Minimizing the Sphere Function

August 6, 2026

A Clean Introduction to Convex Optimization

Introduction

Among all the benchmark functions used in numerical optimization, the Sphere function stands out for its simplicity. It is smooth, strictly convex, and has a single global minimum — making it the perfect starting point for understanding how optimization algorithms behave before tackling more complex, non-convex landscapes. In this article, we implement and visualize the minimization of the Sphere function in two dimensions, comparing a hand-written gradient descent routine against SciPy’s BFGS optimizer.

Mathematical Definition

The Sphere function in $n$ dimensions is defined as:

$$
f(\mathbf{x}) = \sum_{i=1}^{n} x_i^2
$$

For our two-dimensional example, this simplifies to:

$$
f(x_1, x_2) = x_1^2 + x_2^2
$$

Because the function is a simple sum of squares, its gradient has a closed-form expression:

$$
\nabla f(\mathbf{x}) = 2\mathbf{x}
$$

and its Hessian matrix is constant and positive definite:

$$
H(\mathbf{x}) = 2I
$$

where $I$ is the identity matrix. This positive-definite Hessian is exactly what makes the Sphere function strictly convex, guaranteeing that any local minimum found is also the global minimum, located at $\mathbf{x}^* = \mathbf{0}$ with $f(\mathbf{x}^*) = 0$.

Python Implementation

The code below defines the Sphere function and its gradient, runs a custom gradient descent optimizer, cross-validates the result with SciPy’s BFGS method, and visualizes both optimization paths on a 3D surface, a contour map, and a convergence plot.

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

# Dark theme for all plots
plt.style.use('dark_background')

# ----------------------------------------------------------------------
# 1. Define the Sphere function and its gradient
# ----------------------------------------------------------------------
def sphere(x):
    """Sphere function: f(x) = sum(x_i^2)"""
    x = np.asarray(x, dtype=float)
    return np.sum(x ** 2)

def sphere_gradient(x):
    """Gradient of the Sphere function: grad f(x) = 2x"""
    x = np.asarray(x, dtype=float)
    return 2.0 * x

# ----------------------------------------------------------------------
# 2. Custom gradient descent implementation
# ----------------------------------------------------------------------
def gradient_descent(f, grad, x0, learning_rate=0.25, max_iter=60, tol=1e-10):
    x = np.array(x0, dtype=float)
    path = [x.copy()]
    values = [f(x)]

    for _ in range(max_iter):
        g = grad(x)
        x = x - learning_rate * g
        path.append(x.copy())
        values.append(f(x))
        if np.linalg.norm(g) < tol:
            break

    return np.array(path), np.array(values)

# Starting point, deliberately far from the minimum
x0 = np.array([4.0, 4.5])

gd_path, gd_values = gradient_descent(sphere, sphere_gradient, x0,
                                       learning_rate=0.25, max_iter=60)

# ----------------------------------------------------------------------
# 3. Reference solution using SciPy's BFGS optimizer
# ----------------------------------------------------------------------
bfgs_path = [x0.copy()]

def record_step(xk):
    bfgs_path.append(np.array(xk, dtype=float))

result = minimize(sphere, x0, jac=sphere_gradient, method='BFGS',
                   callback=record_step, tol=1e-10)

bfgs_path = np.array(bfgs_path)
bfgs_values = np.array([sphere(p) for p in bfgs_path])

# ----------------------------------------------------------------------
# 4. Print a short numerical summary
# ----------------------------------------------------------------------
print("=" * 60)
print("Sphere Function Minimization Summary")
print("=" * 60)
print(f"Starting point            : {x0}")
print(f"Gradient Descent result   : {gd_path[-1]}, f = {gd_values[-1]:.3e}, iterations = {len(gd_path)-1}")
print(f"BFGS result                : {result.x}, f = {result.fun:.3e}, iterations = {result.nit}")
print(f"True global minimum        : [0. 0.], f = 0.0")
print("=" * 60)

# ----------------------------------------------------------------------
# 5. Visualization
# ----------------------------------------------------------------------
fig = plt.figure(figsize=(20, 6))
fig.patch.set_facecolor('#0d0d0d')

grid_range = np.linspace(-5, 5, 150)
X, Y = np.meshgrid(grid_range, grid_range)
Z = X ** 2 + Y ** 2

# --- 5.1  3D surface with the gradient descent path ---
ax1 = fig.add_subplot(1, 3, 1, projection='3d')
ax1.set_facecolor('#0d0d0d')
surf = ax1.plot_surface(X, Y, Z, cmap='plasma', alpha=0.55,
                         edgecolor='none', antialiased=True)
ax1.plot(gd_path[:, 0], gd_path[:, 1], gd_values,
          color='#00ffea', marker='o', markersize=3, linewidth=2,
          label='Gradient Descent path')
ax1.scatter([0], [0], [0], color='white', s=60, marker='*',
            label='Global minimum')
ax1.set_xlabel('x1')
ax1.set_ylabel('x2')
ax1.set_zlabel('f(x1, x2)')
ax1.set_title('Sphere Function Surface & Descent Path', color='white')
ax1.legend(loc='upper left', fontsize=8)
fig.colorbar(surf, ax=ax1, shrink=0.5, pad=0.1)

# --- 5.2  Contour plot comparing both optimizers ---
ax2 = fig.add_subplot(1, 3, 2)
ax2.set_facecolor('#0d0d0d')
contour = ax2.contour(X, Y, Z, levels=25, cmap='plasma')
ax2.clabel(contour, inline=True, fontsize=7)
ax2.plot(gd_path[:, 0], gd_path[:, 1], color='#00ffea',
         marker='o', markersize=3, linewidth=1.5, label='Gradient Descent')
ax2.plot(bfgs_path[:, 0], bfgs_path[:, 1], color='#ff5df1',
         marker='s', markersize=4, linewidth=1.5, label='BFGS')
ax2.scatter([0], [0], color='white', s=80, marker='*', zorder=5,
            label='Global minimum')
ax2.set_xlabel('x1')
ax2.set_ylabel('x2')
ax2.set_title('Optimization Paths (Contour View)', color='white')
ax2.legend(fontsize=8)

# --- 5.3  Convergence curve (log scale) ---
ax3 = fig.add_subplot(1, 3, 3)
ax3.set_facecolor('#0d0d0d')
ax3.plot(gd_values, color='#00ffea', marker='o', markersize=3,
         label='Gradient Descent')
ax3.plot(bfgs_values, color='#ff5df1', marker='s', markersize=4,
         label='BFGS')
ax3.set_yscale('log')
ax3.set_xlabel('Iteration')
ax3.set_ylabel('f(x)  (log scale)')
ax3.set_title('Convergence Comparison', color='white')
ax3.legend(fontsize=8)
ax3.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

Code Explanation

Function and gradient definitions. sphere() computes $\sum x_i^2$ using np.sum(x ** 2), which works for any dimensionality. sphere_gradient() returns $2\mathbf{x}$ directly, avoiding the need for numerical differentiation and giving both optimizers an exact, noise-free gradient signal.

Custom gradient descent. The gradient_descent() function implements the simplest possible first-order optimizer: at each step it moves in the direction opposite to the gradient, scaled by a fixed learning_rate. Because the Sphere function’s Hessian is $2I$, the update rule simplifies to $x_{k+1} = x_k(1 - 2\eta)$, where $\eta$ is the learning rate. With $\eta = 0.25$, this factor equals $0.5$, so the distance to the origin is halved at every iteration — a textbook example of linear convergence. The loop also includes an early-stopping condition based on gradient norm, so it won’t run needless iterations once it’s essentially converged.

BFGS reference. SciPy’s minimize() with method='BFGS' is a quasi-Newton method that builds an approximate Hessian from gradient information as it iterates. Since the true Hessian of the Sphere function is a constant multiple of the identity, BFGS converges extremely fast — typically within just a handful of steps. The callback argument lets us record every intermediate point xk, which we later plot alongside the gradient descent path.

Why no separate “fast” version is needed. Because the Sphere function and its gradient are trivially cheap to evaluate (just squaring and summing a handful of numbers), this problem never becomes a computational bottleneck even at high iteration counts. The vectorized NumPy operations already run in microseconds per step, so no further speed optimization is necessary here — this keeps the code simple and readable.

Visualization and Interpretation

The figure produced by the script above contains three panels side by side:

3D Surface Plot (left). This shows the paraboloid shape of the Sphere function, colored with the plasma colormap. The cyan trajectory traces the gradient descent path as it spirals down toward the bowl’s bottom, visually confirming the smooth, funnel-like geometry that makes this function so easy to optimize.
Contour Plot (center). Viewed from directly above, the concentric circles represent level sets of equal function value. Both optimization paths are overlaid: the cyan markers show gradient descent taking small, steady steps, while the magenta squares show BFGS reaching the center in dramatically fewer steps thanks to its curvature-aware search direction.
Convergence Plot (right). Plotted on a logarithmic y-axis, this panel makes the difference in convergence speed unmistakable. Gradient descent’s function value decreases in a straight line on the log scale (confirming the theoretical linear/geometric convergence rate), while BFGS’s curve drops almost vertically, reaching machine-precision accuracy in only a few iterations.

The console output — showing the starting point, final coordinates, function values, and iteration counts for both methods — should also confirm that both optimizers converge to the same global minimum at $[0, 0]$ with $f(\mathbf{x}) \approx 0$.

============================================================
Sphere Function Minimization Summary
============================================================
Starting point            : [4.  4.5]
Gradient Descent result   : [1.45519152e-11 1.63709046e-11], f = 4.798e-22, iterations = 38
BFGS result                : [-6.66133815e-16  0.00000000e+00], f = 4.437e-31, iterations = 3
True global minimum        : [0. 0.], f = 0.0
============================================================

Conclusion

The Sphere function may be the “hello world” of optimization problems, but it offers genuine insight: it isolates the raw behavior of an algorithm’s convergence rate without the confounding effects of non-convexity, saddle points, or multiple local minima. By comparing a from-scratch gradient descent implementation against SciPy’s BFGS, we get a clear, visual demonstration of why second-order (curvature-aware) methods so dramatically outperform first-order methods on well-conditioned convex problems — a lesson that carries directly into far more complex optimization landscapes.

Minimizing Himmelblau's Function

August 5, 2026

A Journey Through Multiple Global Minima

Introduction

Optimization problems don’t always have a single “best” answer. Sometimes, a function has multiple points that are all equally optimal — and this is exactly the fascinating property of Himmelblau’s function, a classic benchmark in numerical optimization.

In this post, we’ll explore Himmelblau’s function, find all four of its global minima using Python, and visualize the results with a stunning 3D surface plot.

What is Himmelblau’s Function?

Himmelblau’s function is defined as:

$$
f(x, y) = (x^2 + y - 11)^2 + (x + y^2 - 7)^2
$$

This function is widely used to test optimization algorithms because it has four identical global minima, all with a function value of exactly 0. The minima are located at approximately:

$$
(3.0, 2.0), \quad (-2.805118, 3.131312), \quad (-3.779310, -3.283186), \quad (3.584428, -1.848126)
$$

Because there are multiple equally good solutions, a naive gradient-descent-style optimizer starting from a single point will only find one of these minima — the one closest to its starting position. To find all four, we need to try multiple starting points across the search space.

Strategy

Our approach:

Define Himmelblau’s function in Python.
Use scipy.optimize.minimize with a multi-start strategy — launching the optimizer from many different initial points across the domain.
Cluster the resulting solutions to identify the distinct global minima (removing duplicates found from different starting points).
Visualize the function as a 3D surface and mark the discovered minima.
Also show a 2D contour plot for a clearer top-down view.

Since running the optimizer from a single point is fast, but running it from hundreds of starting points could be slow if done naively, we vectorize the initial point generation with NumPy and use SciPy’s efficient BFGS-based solver (L-BFGS-B) for speed, while keeping the total number of starts modest (a grid of 100 points) so it finishes almost instantly.

Full Source Code

# ==========================================================
# Himmelblau's Function: Multi-Start Optimization + 3D Plot
# ==========================================================
import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D  # noqa: F401 (enables 3D projection)

# ----------------------------------------------------------
# 1. Define Himmelblau's function and its gradient
# ----------------------------------------------------------
def himmelblau(v):
    x, y = v
    return (x**2 + y - 11)**2 + (x + y**2 - 7)**2

def himmelblau_grad(v):
    x, y = v
    dfdx = 4*x*(x**2 + y - 11) + 2*(x + y**2 - 7)
    dfdy = 2*(x**2 + y - 11) + 4*y*(x + y**2 - 7)
    return np.array([dfdx, dfdy])

# ----------------------------------------------------------
# 2. Multi-start optimization to find ALL global minima
# ----------------------------------------------------------
# Generate a grid of starting points covering the search domain
grid_n = 10  # 10x10 = 100 starting points
xs = np.linspace(-6, 6, grid_n)
ys = np.linspace(-6, 6, grid_n)
starts = np.array([[x, y] for x in xs for y in ys])

found_minima = []
tolerance = 1e-4  # distance threshold to consider two minima "the same"

for start in starts:
    result = minimize(
        himmelblau,
        start,
        jac=himmelblau_grad,
        method='L-BFGS-B',
        bounds=[(-6, 6), (-6, 6)]
    )
    if result.success and result.fun < 1e-6:  # only keep true global minima
        point = result.x
        # Check if this point is already in our list (avoid duplicates)
        is_new = True
        for existing in found_minima:
            if np.linalg.norm(point - existing) < tolerance:
                is_new = False
                break
        if is_new:
            found_minima.append(point)

found_minima = np.array(found_minima)

print(f"Number of distinct global minima found: {len(found_minima)}")
print("Coordinates of global minima:")
for i, m in enumerate(found_minima):
    print(f"  Minimum {i+1}: x = {m[0]:.6f}, y = {m[1]:.6f}, f(x,y) = {himmelblau(m):.8f}")

# ----------------------------------------------------------
# 3. Prepare data for 3D surface plot
# ----------------------------------------------------------
X = np.linspace(-6, 6, 200)
Y = np.linspace(-6, 6, 200)
X, Y = np.meshgrid(X, Y)
Z = (X**2 + Y - 11)**2 + (X + Y**2 - 7)**2

# ----------------------------------------------------------
# 4. Plot: 3D surface + 2D contour with minima marked
# ----------------------------------------------------------
fig = plt.figure(figsize=(16, 7))

# --- 3D Surface Plot ---
ax1 = fig.add_subplot(1, 2, 1, projection='3d')
surf = ax1.plot_surface(X, Y, Z, cmap='viridis', alpha=0.85, linewidth=0, antialiased=True)
ax1.scatter(
    found_minima[:, 0], found_minima[:, 1],
    [himmelblau(m) for m in found_minima],
    color='red', s=80, marker='o', label='Global minima', depthshade=False
)
ax1.set_xlabel('x')
ax1.set_ylabel('y')
ax1.set_zlabel('f(x, y)')
ax1.set_title("Himmelblau's Function - 3D Surface")
fig.colorbar(surf, ax=ax1, shrink=0.5, aspect=10)
ax1.legend()

# --- 2D Contour Plot ---
ax2 = fig.add_subplot(1, 2, 2)
contour = ax2.contourf(X, Y, Z, levels=50, cmap='viridis')
ax2.scatter(
    found_minima[:, 0], found_minima[:, 1],
    color='red', s=100, marker='*', edgecolors='white', linewidths=1.2,
    label='Global minima', zorder=5
)
for i, m in enumerate(found_minima):
    ax2.annotate(f"Min {i+1}\n({m[0]:.2f}, {m[1]:.2f})",
                 (m[0], m[1]), textcoords="offset points",
                 xytext=(10, 10), color='white', fontsize=9)
ax2.set_xlabel('x')
ax2.set_ylabel('y')
ax2.set_title("Himmelblau's Function - Contour View")
fig.colorbar(contour, ax=ax2)
ax2.legend(loc='upper right')

plt.tight_layout()
plt.show()

Number of distinct global minima found: 4
Coordinates of global minima:
  Minimum 1: x = 3.000000, y = 2.000000, f(x,y) = 0.00000000
  Minimum 2: x = -2.805118, y = 3.131313, f(x,y) = 0.00000000
  Minimum 3: x = -3.779310, y = -3.283186, f(x,y) = 0.00000000
  Minimum 4: x = 3.584428, y = -1.848127, f(x,y) = 0.00000000

Code Walkthrough

1. Defining the function and its gradient

The himmelblau() function directly implements the mathematical formula shown earlier. We also manually derived the gradient (himmelblau_grad) — the vector of partial derivatives with respect to $x$ and $y$:

$$
\frac{\partial f}{\partial x} = 4x(x^2 + y - 11) + 2(x + y^2 - 7)
$$

$$
\frac{\partial f}{\partial y} = 2(x^2 + y - 11) + 4y(x + y^2 - 7)
$$

Providing the exact gradient (instead of letting SciPy estimate it numerically) makes the optimizer converge faster and more accurately, since it avoids the overhead of finite-difference approximation.

2. Multi-start optimization

Since Himmelblau’s function has four global minima, starting the optimizer from just one point would only ever find one of them. To solve this, we:

Create a 10×10 grid of starting points spanning the domain $[-6, 6] \times [-6, 6]$ — 100 starting points in total.
Run scipy.optimize.minimize with the L-BFGS-B method (a fast quasi-Newton algorithm well suited to smooth, bounded problems) from each starting point.
Keep only results where the function value is essentially zero (result.fun < 1e-6), confirming we’ve truly hit a global minimum rather than some other stationary point.
Deduplicate results: since many nearby starting points converge to the same minimum, we check the Euclidean distance between new solutions and previously found ones, discarding near-duplicates within a 1e-4 tolerance.

This grid-based multi-start approach is a simple but effective way to perform global optimization using a fundamentally local optimizer — and because L-BFGS-B is very efficient, all 100 optimization runs complete in a fraction of a second.

3. Building the surface data

We create a fine mesh grid (200 × 200 points) over the same domain and evaluate the function at every point using vectorized NumPy operations. This gives us the Z array needed to draw a smooth 3D surface and contour map.

4. Visualization

Left panel (3D surface): Shows the overall “landscape” of the function, with two tall peaks and four valley-like basins where the function dips to zero. The red dots mark the discovered global minima, sitting exactly at the bottom of each basin.
Right panel (2D contour): A bird’s-eye view of the same landscape using color gradients (dark = low value, bright = high value). The four star markers, each labeled with its coordinates, make it immediately clear where all four minima are located relative to each other.

Together, these two views make it intuitive to see why Himmelblau’s function is such a popular test case: the four minima are spread across very different regions of the search space, forcing any global optimization algorithm to genuinely explore rather than just “roll downhill” from a single guess.

Conclusion

Himmelblau’s function beautifully illustrates a key challenge in optimization: not all problems have a unique answer. By using a multi-start strategy with SciPy’s L-BFGS-B solver and carefully deduplicating results, we successfully located all four global minima efficiently. The combination of 3D surface and 2D contour plots gives a complete, intuitive picture of the function’s structure — turning an abstract equation into something you can literally see and understand at a glance.

Minimizing the Ackley Function with Python

August 4, 2026

A Practical Global Optimization Example

Optimization problems are everywhere in engineering, machine learning, and data science, but not all objective functions are easy to minimize. Some are riddled with local minima that trap naive algorithms long before they reach the true global minimum. The Ackley function is one of the most famous benchmark functions used to test how well an optimization algorithm can escape these traps. In this article, we’ll explore the Ackley function in depth, minimize it using Python, and visualize the results in both 2D and 3D.

What Is the Ackley Function?

The Ackley function is a widely used benchmark for testing global optimization algorithms because of its nearly flat outer region combined with a large number of local minima near the center. In two dimensions, it is defined as:

$$
f(x, y) = -20 \exp\left(-0.2 \sqrt{0.5(x^2 + y^2)}\right) - \exp\left(0.5(\cos(2\pi x) + \cos(2\pi y))\right) + e + 20
$$

The global minimum is located at $(x, y) = (0, 0)$, where $f(0, 0) = 0$. What makes this function tricky is the combination of an exponential term that creates a huge, nearly flat “bowl” and a cosine term that riddles the surface with countless small local minima. A simple gradient-descent-based method will almost always get stuck in one of these local minima instead of finding the true global minimum at the origin.

Why Use a Global Optimization Algorithm?

Because of the many local minima, a gradient-based local optimizer (like scipy.optimize.minimize with BFGS) is not reliable here unless it starts very close to the global minimum. Instead, we need a global optimization algorithm. In this example, we use scipy.optimize.differential_evolution, a population-based evolutionary algorithm that explores the search space broadly before converging, making it much more robust against local minima traps.

Python Implementation

Below is the complete, self-contained source code. It defines the Ackley function, runs the global optimization, prints the results, and generates both a 3D surface plot and a 2D contour plot with the discovered minimum marked on it. The differential evolution step uses workers=-1 to parallelize across all available CPU cores, which significantly speeds up the search compared to the default single-threaded execution.

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

# ---------------------------------------------------------
# 1. Define the Ackley function
# ---------------------------------------------------------
def ackley(pos, a=20, b=0.2, c=2 * np.pi):
    x, y = pos
    term1 = -a * np.exp(-b * np.sqrt(0.5 * (x**2 + y**2)))
    term2 = -np.exp(0.5 * (np.cos(c * x) + np.cos(c * y)))
    return term1 + term2 + a + np.e

# ---------------------------------------------------------
# 2. Run global optimization (Differential Evolution)
#    workers=-1 enables parallel execution for speed-up
# ---------------------------------------------------------
bounds = [(-5, 5), (-5, 5)]

result = differential_evolution(
    ackley,
    bounds,
    strategy='best1bin',
    maxiter=1000,
    popsize=20,
    tol=1e-10,
    seed=42,
    workers=-1,
    updating='deferred'
)

print("Optimization successful:", result.success)
print("Number of iterations:", result.nit)
print(f"Best solution found: x = {result.x[0]:.6f}, y = {result.x[1]:.6f}")
print(f"Minimum function value: f(x, y) = {result.fun:.10f}")

# ---------------------------------------------------------
# 3. Prepare data grid for visualization
# ---------------------------------------------------------
x_vals = np.linspace(-5, 5, 200)
y_vals = np.linspace(-5, 5, 200)
X, Y = np.meshgrid(x_vals, y_vals)
Z = ackley([X, Y])

# ---------------------------------------------------------
# 4. Plot: 3D surface + 2D contour side by side
# ---------------------------------------------------------
fig = plt.figure(figsize=(16, 7))

# --- 3D Surface Plot ---
ax1 = fig.add_subplot(1, 2, 1, projection='3d')
surf = ax1.plot_surface(X, Y, Z, cmap='viridis', edgecolor='none', alpha=0.9)
ax1.scatter(result.x[0], result.x[1], result.fun,
            color='red', s=100, marker='o', label='Global Minimum')
ax1.set_title('Ackley Function - 3D Surface')
ax1.set_xlabel('x')
ax1.set_ylabel('y')
ax1.set_zlabel('f(x, y)')
ax1.view_init(elev=35, azim=-60)
fig.colorbar(surf, ax=ax1, shrink=0.5, aspect=10)
ax1.legend()

# --- 2D Contour Plot ---
ax2 = fig.add_subplot(1, 2, 2)
contour = ax2.contourf(X, Y, Z, levels=50, cmap='viridis')
ax2.scatter(result.x[0], result.x[1], color='red', s=120,
            marker='*', edgecolor='white', linewidth=1.5,
            label=f'Minimum ({result.x[0]:.3f}, {result.x[1]:.3f})')
ax2.set_title('Ackley Function - 2D Contour Map')
ax2.set_xlabel('x')
ax2.set_ylabel('y')
fig.colorbar(contour, ax=ax2)
ax2.legend()

plt.tight_layout()
plt.show()

Optimization successful: True
Number of iterations: 123
Best solution found: x = 0.000000, y = 0.000000
Minimum function value: f(x, y) = 0.0000000000

Code Walkthrough

1. Defining the Ackley Function

The ackley() function directly implements the mathematical formula introduced earlier. It takes a 2-element array pos (representing $x$ and $y$) and returns a single scalar value. The parameters a=20, b=0.2, and c=2π are the standard constants used in the canonical definition of the function, and keeping them as arguments makes the function reusable for variations of the benchmark.

2. Global Optimization with Differential Evolution

differential_evolution maintains a population of candidate solutions and iteratively “evolves” them using mutation, crossover, and selection — a strategy inspired by biological evolution. This makes it far more resistant to getting trapped in the countless small local minima created by the cosine terms compared to gradient-based methods.

Key parameters worth understanding:

bounds: Defines the search space, here $[-5, 5]$ for both $x$ and $y$, which comfortably contains the interesting region of the function.
strategy='best1bin': A standard and reliable mutation strategy that tends to converge well on smooth, bowl-shaped benchmark functions like this one.
popsize=20: Controls how many candidate solutions are evaluated per generation. Larger values improve robustness at the cost of speed.
tol=1e-10: A tight convergence tolerance to make sure the algorithm doesn’t stop prematurely before finding a highly precise minimum.
workers=-1 and updating='deferred': These two settings together enable multi-core parallel evaluation of the population, which is the key speed-up trick here. Without them, differential_evolution evaluates each candidate solution sequentially; with them, all available CPU cores are used simultaneously, cutting runtime significantly — especially valuable if you extend this example to higher dimensions or more expensive objective functions.

After the optimization finishes, the script prints whether the run converged successfully, how many generations it took, and the best $(x, y)$ pair found along with its function value, which should be extremely close to the true global minimum of $0$ at $(0, 0)$.

3. Building the Visualization Grid

To visualize the function, we create a dense $200 \times 200$ grid of $(x, y)$ points spanning the search space using np.meshgrid, then evaluate the Ackley function across the entire grid at once using NumPy’s vectorized operations. This is far faster than looping over each point individually in Python.

4. Two Complementary Plots

3D Surface Plot: This shows the overall shape of the function — a broad, nearly flat plateau near the edges that drops sharply into a narrow, spiky funnel near the center. The red marker highlights exactly where the optimizer converged, letting you visually confirm it landed at the bottom of the funnel rather than on one of the small surrounding bumps.
2D Contour Plot: This provides a top-down view using color gradients to represent function value, making it easy to see the ring-like pattern of local minima surrounding the true global minimum. The star marker again shows the optimizer’s final solution, with its exact coordinates included in the legend.

Together, these two plots make it intuitive to understand both the global structure of the Ackley function and why gradient-based local search methods struggle with it — the small ripples visible in the contour plot near the center are individual local minima that a naive optimizer could easily get stuck in.

Conclusion

The Ackley function is a great illustration of why the choice of algorithm matters as much as the implementation when tackling non-convex optimization problems. By using differential_evolution, a population-based global optimizer, we reliably converge to the true minimum at the origin — even though the function’s surface is filled with deceptive local minima. The parallelized workers=-1 setting also demonstrates a practical, simple way to speed up evolutionary optimization runs on any multi-core machine, which becomes increasingly valuable as the problem dimensionality grows.

Solving the Rastrigin Function's Global Minimization Problem with Python

August 3, 2026

Conquering Chaos

If you’ve ever worked in optimization, you know that not all problems are created equal. Some cost functions have a single, elegant global minimum that any basic gradient descent can find in seconds. Others are a minefield of local minima designed to trap naive algorithms. Today, we’re tackling one of the most famous examples of the latter: the Rastrigin function.

This benchmark function is a rite of passage for anyone studying metaheuristic optimization — genetic algorithms, particle swarm optimization, differential evolution, and simulated annealing all get tested against it. Let’s break down exactly why it’s so difficult, and then solve it properly in Python.

What Makes the Rastrigin Function So Difficult?

The Rastrigin function is defined as:

$$
f(\mathbf{x}) = An + \sum_{i=1}^{n} \left[ x_i^2 - A\cos(2\pi x_i) \right]
$$

where $A = 10$ and $\mathbf{x} = (x_1, x_2, \dots, x_n) \in [-5.12, 5.12]^n$.

At first glance, this looks simple — it’s just a quadratic bowl ($x_i^2$) with a cosine term layered on top. But that cosine term is the troublemaker. It creates a highly regular pattern of ripples across the entire search space, generating an enormous number of local minima that grow exponentially with dimension $n$. The global minimum is always at $\mathbf{x} = \mathbf{0}$, where $f(\mathbf{0}) = 0$, but a naive local optimizer starting anywhere off-center will almost certainly get stuck in one of the surrounding “dips” long before it ever finds the true bottom.

This makes the Rastrigin function the perfect testbed for demonstrating the difference between local optimization (which gets fooled) and global optimization (which doesn’t).

Our Approach

In this article, we’ll do four things:

Visualize the Rastrigin landscape in 3D to see why it’s so treacherous.
Demonstrate how a standard local optimizer (L-BFGS-B) gets trapped in local minima.
Solve the problem properly using Differential Evolution (DE), a population-based global optimizer.
Provide a vectorized, high-speed version of the DE solver for higher-dimensional cases, since naive implementations can be painfully slow.

Let’s get into the code.

Full Source Code

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.optimize import differential_evolution, minimize
import time

# ==========================================================
# 1. Define the Rastrigin function
# ==========================================================

def rastrigin_scalar(x, A=10):
    """Standard scalar version: x is a 1D array of shape (n_dim,)"""
    x = np.asarray(x)
    n = x.shape[0]
    return A * n + np.sum(x**2 - A * np.cos(2 * np.pi * x))

def rastrigin_vectorized(X, A=10):
    """Vectorized version for fast batch evaluation.
    X has shape (n_dim, n_population) as required by scipy's
    vectorized differential_evolution."""
    return A * X.shape[0] + np.sum(X**2 - A * np.cos(2 * np.pi * X), axis=0)

# ==========================================================
# 2. Visualize the Rastrigin landscape (n = 2) in 3D and contour
# ==========================================================

bound = 5.12
res = 300
x = np.linspace(-bound, bound, res)
y = np.linspace(-bound, bound, res)
X, Y = np.meshgrid(x, y)
Z = 20 + (X**2 - 10 * np.cos(2 * np.pi * X)) + (Y**2 - 10 * np.cos(2 * np.pi * Y))

fig = plt.figure(figsize=(16, 6))

ax1 = fig.add_subplot(1, 2, 1, projection='3d')
surf = ax1.plot_surface(X, Y, Z, cmap='viridis', linewidth=0, antialiased=True, alpha=0.95)
ax1.set_title("Rastrigin Function - 3D Surface (n = 2)")
ax1.set_xlabel("x1")
ax1.set_ylabel("x2")
ax1.set_zlabel("f(x1, x2)")
ax1.view_init(elev=35, azim=45)
fig.colorbar(surf, ax=ax1, shrink=0.5, aspect=10)

ax2 = fig.add_subplot(1, 2, 2)
contour = ax2.contourf(X, Y, Z, levels=50, cmap='viridis')
ax2.plot(0, 0, 'r*', markersize=18, label='Global Minimum (0,0)')
ax2.set_title("Rastrigin Function - Contour Map")
ax2.set_xlabel("x1")
ax2.set_ylabel("x2")
ax2.legend()
fig.colorbar(contour, ax=ax2)

plt.tight_layout()
plt.show()

# ==========================================================
# 3. Demonstrate the trap: local optimization from random starts
# ==========================================================

n_trials = 30
local_minima_found = []
np.random.seed(0)

for _ in range(n_trials):
    x0 = np.random.uniform(-bound, bound, size=2)
    res_local = minimize(
        rastrigin_scalar, x0,
        method='L-BFGS-B',
        bounds=[(-bound, bound)] * 2
    )
    local_minima_found.append(res_local.fun)

plt.figure(figsize=(8, 5))
plt.hist(local_minima_found, bins=15, color='steelblue', edgecolor='black')
plt.axvline(0, color='red', linestyle='--', linewidth=2, label='True Global Minimum = 0')
plt.xlabel("Objective value found by L-BFGS-B")
plt.ylabel("Frequency (out of 30 random starts)")
plt.title("Local Optimizer Gets Trapped: Distribution of Results")
plt.legend()
plt.show()

print(f"Best value found across {n_trials} local searches: {min(local_minima_found):.4f}")
print(f"Worst value found: {max(local_minima_found):.4f}")
print(f"Success rate (f < 0.01): {sum(v < 0.01 for v in local_minima_found)}/{n_trials}")

# ==========================================================
# 4. Global optimization with Differential Evolution (10-D problem)
# ==========================================================

dim = 10
bounds = [(-bound, bound)] * dim
history_scalar = []

def callback_de(xk, convergence):
    history_scalar.append(rastrigin_scalar(xk))

start = time.time()
result_scalar = differential_evolution(
    rastrigin_scalar,
    bounds,
    strategy='best1bin',
    maxiter=1000,
    popsize=20,
    tol=1e-8,
    mutation=(0.5, 1.0),
    recombination=0.7,
    seed=42,
    callback=callback_de,
    polish=True
)
elapsed_scalar = time.time() - start

print("\n--- Standard (scalar) Differential Evolution ---")
print(f"Global minimum found: f(x*) = {result_scalar.fun:.6e}")
print(f"x* = {np.round(result_scalar.x, 4)}")
print(f"Elapsed time: {elapsed_scalar:.3f} s")
print(f"Generations run: {result_scalar.nit}")

# ==========================================================
# 5. High-speed version: vectorized Differential Evolution
# ==========================================================

start = time.time()
result_fast = differential_evolution(
    rastrigin_vectorized,
    bounds,
    strategy='best1bin',
    maxiter=1000,
    popsize=20,
    tol=1e-8,
    mutation=(0.5, 1.0),
    recombination=0.7,
    seed=42,
    polish=True,
    vectorized=True,
    updating='deferred'
)
elapsed_fast = time.time() - start

print("\n--- Vectorized (high-speed) Differential Evolution ---")
print(f"Global minimum found: f(x*) = {result_fast.fun:.6e}")
print(f"x* = {np.round(result_fast.x, 4)}")
print(f"Elapsed time: {elapsed_fast:.3f} s")
print(f"Generations run: {result_fast.nit}")
print(f"\nSpeedup factor: {elapsed_scalar / elapsed_fast:.2f}x faster")

# ==========================================================
# 6. Visualize convergence
# ==========================================================

plt.figure(figsize=(10, 5))
plt.plot(history_scalar, color='crimson', linewidth=2)
plt.yscale('log')
plt.xlabel("Generation")
plt.ylabel("Best f(x) found so far (log scale)")
plt.title(f"Convergence of Differential Evolution on {dim}-D Rastrigin Function")
plt.grid(True, which='both', linestyle='--', alpha=0.6)
plt.show()

Best value found across 30 local searches: 0.9950
Worst value found: 49.7474
Success rate (f < 0.01): 0/30

--- Standard (scalar) Differential Evolution ---
Global minimum found: f(x*) = 9.949591e-01
x* = [-0.    -0.    -0.    -0.     0.995 -0.    -0.    -0.    -0.    -0.   ]
Elapsed time: 10.238 s
Generations run: 518

--- Vectorized (high-speed) Differential Evolution ---
Global minimum found: f(x*) = 0.000000e+00
x* = [ 0. -0.  0.  0.  0.  0.  0.  0. -0.  0.]
Elapsed time: 1.981 s
Generations run: 778

Speedup factor: 5.17x faster

Code Walkthrough

1. Two flavors of the objective function

Notice that we defined two versions of the Rastrigin function: rastrigin_scalar and rastrigin_vectorized. This isn’t redundant — it’s the key to the performance story of this article.

rastrigin_scalar takes a single candidate solution (a 1D array) and returns a single number. This is the natural way to write an objective function, and it’s what most scipy.optimize routines expect by default.

rastrigin_vectorized, on the other hand, accepts an entire population of candidates at once — a 2D array where each column is one candidate solution — and returns all their objective values in one shot using NumPy’s broadcasting. This eliminates the Python-level loop overhead that occurs when an optimizer evaluates hundreds of candidates one at a time.

2. Visualizing the landscape

The 3D surface plot uses plot_surface on a 300×300 grid over the 2D search space $[-5.12, 5.12]^2$. You’ll immediately notice the “egg carton” texture — countless symmetric bumps surrounding a single, slightly deeper well at the origin. The contour map on the right shows the same thing from a bird’s-eye view, making the sheer number of local minima even more apparent. That red star marks the one true global minimum among dozens of decoys.

3. Proving the trap is real

Before jumping to the “solution,” we first prove the problem exists. We run scipy.optimize.minimize with the L-BFGS-B method (a fast, gradient-based local optimizer) from 30 different random starting points. L-BFGS-B is excellent at descending smoothly to the nearest minimum — but “nearest” is the operative word. The resulting histogram typically shows results scattered across a wide range of nonzero values, with only a small fraction landing near the true minimum of 0. This is a direct, visual demonstration of why local search alone is unreliable on this function.

4. Solving it with Differential Evolution

Differential Evolution (DE) is a population-based, gradient-free metaheuristic. Instead of following a single point downhill, it maintains an entire population of candidate solutions that evolve generation by generation through mutation, crossover, and selection. Because it explores many regions of the search space simultaneously, it’s far more resistant to getting stuck in any single local minimum.

We apply it here to a 10-dimensional version of the Rastrigin function ($n=10$) — a much harder instance than the 2D visualization, with an astronomically larger number of local minima. Key parameters:

strategy='best1bin': a classic and robust DE mutation/crossover strategy.
popsize=20: population size multiplier (actual population = popsize × dim).
mutation=(0.5, 1.0): dithering range for the mutation factor, which helps avoid premature convergence.
polish=True: after DE converges, scipy runs a quick local L-BFGS-B polish on the best solution to sharpen precision.
callback=callback_de: lets us record the best objective value after every generation, which we use later for the convergence plot.

5. The high-speed version

Standard DE evaluates the objective function once per candidate per generation using a Python-level loop internally (or optionally via multiprocessing with workers=-1, which carries process-spawning overhead that isn’t worth it for a cheap function like this one).

Instead, we use scipy’s vectorized=True mode combined with updating='deferred'. This passes the entire population to rastrigin_vectorized in a single NumPy call per generation, letting NumPy’s compiled C backend handle the heavy lifting instead of Python’s interpreter loop. For cheap-to-evaluate functions like Rastrigin, this is typically several times faster than both the naive scalar approach and multiprocessing-based parallelism, since it avoids both interpreter overhead and inter-process communication costs. The speedup factor is printed directly in the output so you can see the improvement on your own machine.

6. Reading the convergence plot

The final plot shows the best objective value found at each generation, plotted on a logarithmic y-axis (since the values shrink by orders of magnitude). You should see a characteristic staircase pattern: long flat stretches where DE is exploring without improvement, punctuated by sharp drops when it discovers a better region of the search space. By the final generations, the curve should flatten out near $10^{-8}$ to $10^{-10}$ — effectively zero, confirming that the algorithm has converged to the true global minimum at the origin.

Interpreting the Results

Once you run this in your own environment, here’s what to look for in the output:

The 3D/contour plots should confirm visually just how deceptive this landscape is — dozens of local dips surrounding one true minimum.
The local-search histogram should show that L-BFGS-B rarely finds the true minimum on its own; most runs land at nonzero values corresponding to nearby local minima.
Differential Evolution’s final result (result_scalar.fun and result_fast.fun) should be extremely close to 0, with the solution vector x* close to all zeros — even in 10 dimensions.
The speedup factor printed at the end quantifies how much faster the vectorized approach is compared to the naive scalar approach on your hardware.

Key Takeaways

The Rastrigin function is a small piece of code with an outsized lesson: gradient-based local optimizers are only as good as their starting point when the landscape is riddled with local minima. Population-based global optimizers like Differential Evolution trade some computational cost for dramatically better robustness — and with proper vectorization, that computational cost can be kept surprisingly low.

If you’re building optimization pipelines for real-world problems — hyperparameter tuning, engineering design, portfolio optimization — and you suspect your loss landscape might be non-convex or multimodal, this is exactly the kind of test you should run before trusting a purely local method.

Minimizing the Rosenbrock Function (Banana Function) with Python

August 2, 2026

What Is the Rosenbrock Function?

The Rosenbrock function is one of the most famous test problems in numerical optimization. It’s often called the “banana function” because its contour lines form a curved, banana-shaped valley. While the function looks simple, the valley is narrow and curved, which makes it deceptively difficult for many optimization algorithms to converge on efficiently.

The standard two-dimensional form is:

$$
f(x, y) = (a - x)^2 + b(y - x^2)^2
$$

with the classic parameter choice $a = 1$, $b = 100$. The global minimum sits at $(x, y) = (1, 1)$, where $f(x,y) = 0$.

The gradient, which we’ll use for faster convergence, is:

$$
\frac{\partial f}{\partial x} = -2(a - x) - 4bx(y - x^2)
$$

$$
\frac{\partial f}{\partial y} = 2b(y - x^2)
$$

Because the valley curves and is very narrow, simple gradient descent tends to zig-zag slowly along the valley floor, taking a huge number of iterations. This makes the Rosenbrock function a great case study for comparing a “naive” optimization approach against a properly optimized one.

The Approach

In this article we’ll do three things:

Implement a naive gradient descent solver from scratch in pure Python and time it.
Implement a fast, vectorized solver using scipy.optimize with the analytic gradient supplied, and time it.
Visualize both optimization paths on a 3D surface plot and a 2D contour plot, so you can literally see why one approach struggles and the other doesn’t.

Full Source Code

# ==========================================================
# Rosenbrock Function Minimization: Naive vs. Fast Approach
# ==========================================================

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

# ----------------------------------------------------------
# 1. Define the Rosenbrock function and its gradient
# ----------------------------------------------------------
A = 1.0
B = 100.0

def rosenbrock(v):
    x, y = v
    return (A - x)**2 + B * (y - x**2)**2

def rosenbrock_grad(v):
    x, y = v
    dfdx = -2 * (A - x) - 4 * B * x * (y - x**2)
    dfdy = 2 * B * (y - x**2)
    return np.array([dfdx, dfdy])

# ----------------------------------------------------------
# 2. Naive gradient descent (pure Python loop, fixed step)
# ----------------------------------------------------------
def naive_gradient_descent(start, lr=0.001, n_iter=20000):
    path = [np.array(start, dtype=float)]
    x = np.array(start, dtype=float)
    for _ in range(n_iter):
        grad = rosenbrock_grad(x)
        x = x - lr * grad
        path.append(x.copy())
    return np.array(path)

start_point = np.array([-1.5, 2.0])

t0 = time.time()
naive_path = naive_gradient_descent(start_point, lr=0.001, n_iter=20000)
t1 = time.time()
naive_time = t1 - t0
naive_result = naive_path[-1]

# ----------------------------------------------------------
# 3. Fast optimization using scipy (BFGS with analytic gradient)
# ----------------------------------------------------------
bfgs_path = [start_point.copy()]

def record_path(xk):
    bfgs_path.append(xk.copy())

t2 = time.time()
res = minimize(
    rosenbrock,
    start_point,
    jac=rosenbrock_grad,
    method="BFGS",
    callback=record_path,
    options={"gtol": 1e-8}
)
t3 = time.time()
bfgs_time = t3 - t2
bfgs_path = np.array(bfgs_path)

# ----------------------------------------------------------
# 4. Print a comparison summary
# ----------------------------------------------------------
print("=" * 55)
print("Naive Gradient Descent")
print(f"  Iterations   : {len(naive_path) - 1}")
print(f"  Final point  : ({naive_result[0]:.5f}, {naive_result[1]:.5f})")
print(f"  Final f(x,y) : {rosenbrock(naive_result):.8f}")
print(f"  Elapsed time : {naive_time:.4f} sec")
print("=" * 55)
print("SciPy BFGS (analytic gradient)")
print(f"  Iterations   : {len(bfgs_path) - 1}")
print(f"  Final point  : ({res.x[0]:.5f}, {res.x[1]:.5f})")
print(f"  Final f(x,y) : {rosenbrock(res.x):.8f}")
print(f"  Elapsed time : {bfgs_time:.4f} sec")
print("=" * 55)
speedup = naive_time / bfgs_time if bfgs_time > 0 else float("inf")
print(f"BFGS was about {speedup:.1f}x faster in wall-clock time,")
print(f"and used only {len(bfgs_path)-1} iterations vs {len(naive_path)-1}.")

# ----------------------------------------------------------
# 5. Build a grid for visualization (vectorized, fast)
# ----------------------------------------------------------
x_range = np.linspace(-2, 2, 400)
y_range = np.linspace(-1, 3, 400)
X, Y = np.meshgrid(x_range, y_range)
Z = (A - X)**2 + B * (Y - X**2)**2  # fully vectorized, no Python loops

# ----------------------------------------------------------
# 6. 3D surface plot of the Rosenbrock landscape
# ----------------------------------------------------------
fig = plt.figure(figsize=(20, 8))

ax1 = fig.add_subplot(1, 2, 1, projection="3d")
surf = ax1.plot_surface(
    X, Y, np.log1p(Z),  # log1p compresses the huge range for visibility
    cmap="viridis",
    linewidth=0,
    antialiased=True,
    alpha=0.9
)
ax1.set_title("Rosenbrock Function Surface (log1p scale)", fontsize=13)
ax1.set_xlabel("x")
ax1.set_ylabel("y")
ax1.set_zlabel("log(1 + f(x, y))")
ax1.view_init(elev=35, azim=-60)
fig.colorbar(surf, ax=ax1, shrink=0.6, aspect=12, label="log(1 + f)")

# ----------------------------------------------------------
# 7. 2D contour plot with both optimization paths overlaid
# ----------------------------------------------------------
ax2 = fig.add_subplot(1, 2, 2)
levels = np.logspace(-1, 3.5, 30)
contour = ax2.contour(X, Y, Z, levels=levels, cmap="viridis", norm=None)
ax2.clabel(contour, inline=True, fontsize=6, fmt="%.0f")

ax2.plot(naive_path[:, 0], naive_path[:, 1], "r-", linewidth=1.2,
          label=f"Naive GD ({len(naive_path)-1} steps)")
ax2.plot(bfgs_path[:, 0], bfgs_path[:, 1], "b-o", linewidth=1.8,
          markersize=3, label=f"BFGS ({len(bfgs_path)-1} steps)")

ax2.plot(*start_point, "ks", markersize=8, label="Start")
ax2.plot(1, 1, "g*", markersize=16, label="True Minimum (1, 1)")

ax2.set_title("Optimization Paths on Contour Map", fontsize=13)
ax2.set_xlabel("x")
ax2.set_ylabel("y")
ax2.legend(loc="upper left", fontsize=9)
ax2.set_xlim(-2, 2)
ax2.set_ylim(-1, 3)

plt.tight_layout()
plt.show()

Code Walkthrough

Function and gradient definitions (Section 1): rosenbrock computes $f(x,y)$ directly from the formula. rosenbrock_grad computes the analytic partial derivatives shown above. Supplying the exact gradient instead of letting the optimizer estimate it numerically (via finite differences) is one of the biggest speed wins available — it avoids extra function evaluations and gives more accurate search directions.

Naive gradient descent (Section 2): This is the “textbook” implementation: at every step, move a small distance (lr = 0.001) in the direction opposite the gradient. Because the Rosenbrock valley is long, narrow, and curved, a fixed small step size forces the algorithm to take tens of thousands of tiny, zig-zagging steps just to crawl along the valley floor. We deliberately cap it at 20,000 iterations here so it finishes in reasonable time, but notice in the printed output that it still hasn’t fully converged.

Fast optimization with SciPy BFGS (Section 3): scipy.optimize.minimize with method="BFGS" uses a quasi-Newton approach: it builds up an approximation of the curvature (the Hessian) as it goes, allowing it to take much smarter, adaptively-sized steps. Combined with the analytic gradient passed via jac=rosenbrock_grad, this converges in typically fewer than 40 iterations — a dramatic reduction compared to the naive loop. The callback function records every intermediate point so we can plot the path afterward.

Timing comparison (Section 4): We use time.time() before and after each optimization to measure wall-clock performance, then print a clear side-by-side summary including how many iterations each method needed and the resulting speedup factor.

Vectorized grid for plotting (Section 5): Instead of looping over every $(x, y)$ pair in Python (which would be extremely slow for a 400×400 grid = 160,000 points), we use np.meshgrid combined with vectorized NumPy arithmetic to compute the entire surface Z in one shot. This is the same “vectorization” principle that makes the BFGS run fast — pushing work into optimized C-level array operations instead of Python-level loops.

3D surface plot (Section 6): We plot log1p(Z) instead of raw Z because the Rosenbrock function’s values span an enormous range (from 0 near the minimum to tens of thousands at the edges of the plotted domain). Without the log compression, the banana-shaped valley would be invisible, flattened under one dominant peak. The log1p (log of 1+Z) keeps the shape interpretable while safely handling values near zero.

2D contour plot with paths (Section 7): Log-spaced contour levels (np.logspace) again handle the wide value range gracefully, clearly revealing the curved valley. On top of the contours we plot the naive gradient descent path in red and the BFGS path in blue — this is where the difference between the two approaches becomes visually obvious. You’ll see the red path creeping slowly along the valley in tiny steps, while the blue path jumps directly toward the true minimum, marked with a green star at $(1, 1)$.

Result

Running the code above in Google Colab will print a timing/iteration comparison table, followed by a combined figure showing the 3D landscape on the left and the contour map with both optimization paths on the right.

=======================================================
Naive Gradient Descent
  Iterations   : 20000
  Final point  : (0.99982, 0.99964)
  Final f(x,y) : 0.00000003
  Elapsed time : 0.2666 sec
=======================================================
SciPy BFGS (analytic gradient)
  Iterations   : 38
  Final point  : (1.00000, 1.00000)
  Final f(x,y) : 0.00000000
  Elapsed time : 0.0106 sec
=======================================================
BFGS was about 25.2x faster in wall-clock time,
and used only 38 iterations vs 20000.

Takeaways

This example highlights two important optimization lessons that generalize far beyond the Rosenbrock function:

Step size matters enormously. A fixed small learning rate is “safe” but painfully slow on ill-conditioned, curved landscapes like this one.
Using derivative information and adaptive step sizing (as BFGS does) can turn tens of thousands of iterations into just a few dozen, while also landing much closer to the true minimum.

The Rosenbrock function remains a staple benchmark in optimization research precisely because it exposes these weaknesses so clearly — if an algorithm handles the banana-shaped valley efficiently, it’s usually a strong sign that it will perform well on other difficult, non-convex problems too.

The Physics: From Energy to Equation

Discretizing the Problem

The Example Problem

Python Source Code

Code Walkthrough

Why Conjugate Gradient Wins on Speed

Understanding the Results

Closing Thoughts

A Two-Bar Truss Case Study

The Structure

The Mechanics

The Optimization Problem

Python Implementation

Code Walkthrough

Visualizing the Results

Figure 1 — The 3D Weight Surface

Figure 2 — Design Space, Feasible Region, and Constraint Boundaries

Figure 3 — Convergence History

Figure 4 — The 1D Reduced Problem

Takeaways

A Spring Chain Model in Python

1. The Physical Model

2. Full Python Source Code

3. Code Walkthrough

Physical model functions

Naive loop vs. vectorized grid computation

Three independent solving methods

Visualization

4. Interpreting the Results

5. Extending the Model

A Global Optimization Walkthrough

What Is the Goldstein-Price Function?

Strategy: Global Search First, Local Refinement Second

The Source Code

Code Walkthrough

Visualizing the Function

Graph Explanation

A Journey Through a Deceptively Simple Optimization Landscape

Introduction

The Mathematics of the Beale Function

Optimization Strategy

Source Code

Code Walkthrough

Visualizing the Optimization Landscape

Interpreting the Visualization

Conclusion

A Clean Introduction to Convex Optimization

Introduction

Mathematical Definition

Python Implementation

Code Explanation

Visualization and Interpretation

Conclusion

A Journey Through Multiple Global Minima

Introduction

What is Himmelblau’s Function?

Strategy

Full Source Code

Code Walkthrough

1. Defining the function and its gradient

2. Multi-start optimization

3. Building the surface data

4. Visualization

Conclusion

A Practical Global Optimization Example

What Is the Ackley Function?

Why Use a Global Optimization Algorithm?

Python Implementation

Code Walkthrough

1. Defining the Ackley Function

2. Global Optimization with Differential Evolution

3. Building the Visualization Grid

4. Two Complementary Plots

Conclusion

Conquering Chaos

What Makes the Rastrigin Function So Difficult?

Our Approach

Full Source Code

Code Walkthrough

1. Two flavors of the objective function