{
"cells": [
{
"cell_type": "markdown",
"id": "0f8f26bb",
"metadata": {},
"source": [
"# 6 · Why Line Search Matters\n",
"\n",
"A Newton *direction* $-H^{-1}\\nabla E$ points the right way, but the *step\n",
"length* still matters. Too large a step overshoots, inverts elements, and the\n",
"energy **explodes** instead of decreasing.\n",
"\n",
"**Backtracking line search** starts from a full step and shrinks it until the\n",
"energy actually decreases (the Armijo condition). To make the failure easy to\n",
"watch, we move the handle **slowly right and back** and take a few Newton\n",
"iterations per frame — once with a fixed, too-large step, once with line\n",
"search."
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "eb10e2dc",
"metadata": {
"execution": {
"iopub.execute_input": "2026-06-09T05:22:32.682889Z",
"iopub.status.busy": "2026-06-09T05:22:32.682704Z",
"iopub.status.idle": "2026-06-09T05:22:33.192540Z",
"shell.execute_reply": "2026-06-09T05:22:33.192245Z"
}
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import numpy as np\n",
"import scipy as sp\n",
"import matplotlib.pyplot as plt\n",
"import simkit\n",
"from simkit import backtracking_line_search\n",
"import utils"
]
},
{
"cell_type": "markdown",
"id": "bb7247ad",
"metadata": {},
"source": [
"## A thick, stiff, nearly-incompressible beam\n",
"\n",
"Same objective as tutorial 5 but **thicker** (so bending stores more energy) and\n",
"**nearly incompressible** ($\\nu = 0.49$ — resisting volume change makes the\n",
"energy much stiffer), so an over-eager step bites fast and the value of line\n",
"search is obvious. The handle (the full right edge) follows a there-and-back\n",
"ramp.\n",
"\n",
"`BeamObjective` composes the now-familiar terms — the Neo-Hookean\n",
"material, a fixed **pin** spring on the left edge, and a movable **handle**\n",
"spring on the right edge — and sums them **one per line** in\n",
"`energy / gradient / hessian`."
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "9e184383",
"metadata": {
"execution": {
"iopub.execute_input": "2026-06-09T05:22:33.193793Z",
"iopub.status.busy": "2026-06-09T05:22:33.193707Z",
"iopub.status.idle": "2026-06-09T05:22:33.195966Z",
"shell.execute_reply": "2026-06-09T05:22:33.195747Z"
}
},
"outputs": [],
"source": [
"class BeamObjective:\n",
" \"\"\"Minimize elastic(x) + pin spring + handle spring for a thick beam.\"\"\"\n",
"\n",
" def __init__(self, X, T, K=1e6):\n",
" mu, lam = simkit.ympr_to_lame(50.0, 0.49) # high Poisson ratio\n",
" self.p = utils.precompute(X, T, mu=mu, lam=lam) # J, vol, masses, ...\n",
" self.psi = utils.make_material(\"Neo-Hookean\") # elastic term\n",
" self.pin = utils.PenaltySpring(self.p.n, self.p.dim, K)\n",
" self.handle = utils.PenaltySpring(self.p.n, self.p.dim, K)\n",
"\n",
" def set_pin(self, idx, targets): self.pin.set(idx, targets); return self\n",
" def set_handle(self, idx, targets): self.handle.set(idx, targets); return self\n",
"\n",
" def energy(self, x):\n",
" E_elastic = self.psi.energy(x, self.p)\n",
" E_pin = self.pin.energy(x)\n",
" E_handle = self.handle.energy(x)\n",
" return E_elastic + E_pin + E_handle\n",
"\n",
" def gradient(self, x):\n",
" g_elastic = self.psi.gradient(x, self.p)\n",
" g_pin = self.pin.gradient(x)\n",
" g_handle = self.handle.gradient(x)\n",
" return g_elastic + g_pin + g_handle\n",
"\n",
" def hessian(self, x):\n",
" H_elastic = self.psi.hessian(x, self.p)\n",
" H_pin = self.pin.hessian(x)\n",
" H_handle = self.handle.hessian(x)\n",
" return H_elastic + H_pin + H_handle"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "f1b9edf0",
"metadata": {
"execution": {
"iopub.execute_input": "2026-06-09T05:22:33.196890Z",
"iopub.status.busy": "2026-06-09T05:22:33.196837Z",
"iopub.status.idle": "2026-06-09T05:22:33.198585Z",
"shell.execute_reply": "2026-06-09T05:22:33.198383Z"
}
},
"outputs": [],
"source": [
"X, T = utils.triangulated_grid(nx=16, ny=7, width=2.0, height=0.6) # thick beam\n",
"n, dim = X.shape\n",
"\n",
"pin_idx = np.where(X[:, 0] <= X[:, 0].min() + 1e-6)[0]\n",
"right_idx = np.where(X[:, 0] >= X[:, 0].max() - 1e-6)[0]\n",
"\n",
"ramp = np.concatenate([np.linspace(0, 1, 22), np.linspace(1, 0, 22)]) # right then back\n",
"offsets = [np.array([1.0 * r, 0.0]) for r in ramp]"
]
},
{
"cell_type": "markdown",
"id": "2afeeac2",
"metadata": {},
"source": [
"## The quasi-static driver\n",
"\n",
"At each handle position we take a few Newton iterations. The only difference\n",
"between the two runs is the step: a fixed `STEP` (deliberately too big) versus\n",
"the backtracking step that guarantees the energy decreases.\n",
"\n",
"We solve the linear system $H\\,\\Delta x = -g$ for the Newton direction, then\n",
"either take a fixed multiple of it or let `backtracking_line_search` pick a safe\n",
"step length $\\alpha$."
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "df45bdc4",
"metadata": {
"execution": {
"iopub.execute_input": "2026-06-09T05:22:33.199490Z",
"iopub.status.busy": "2026-06-09T05:22:33.199436Z",
"iopub.status.idle": "2026-06-09T05:22:33.318631Z",
"shell.execute_reply": "2026-06-09T05:22:33.318412Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"fixed step (no line search): blew up at frame 9 of 44\n",
"line search: completed all 44 frames\n"
]
}
],
"source": [
"DIVERGE = 1e3\n",
"\n",
"def drive(use_line_search, STEP=2.0, iters_per_frame=2):\n",
" beam = BeamObjective(X, T).set_pin(pin_idx, X[pin_idx])\n",
" x = X.flatten().reshape(-1, 1).copy()\n",
" states, blew_at = [], None\n",
" for k, off in enumerate(offsets):\n",
" beam.set_handle(right_idx, X[right_idx] + off)\n",
" for _ in range(iters_per_frame):\n",
" g = beam.gradient(x)\n",
" H = beam.hessian(x)\n",
" dx = sp.sparse.linalg.spsolve(H.tocsc(), -g).reshape(-1, 1)\n",
" if use_line_search:\n",
" alpha = backtracking_line_search(beam.energy, x, g, dx)[0]\n",
" else:\n",
" alpha = STEP\n",
" x = x + alpha * dx\n",
" states.append(x.reshape(n, dim).copy())\n",
" if not np.isfinite(np.abs(x).max()) or np.abs(x).max() > DIVERGE:\n",
" blew_at = k\n",
" break\n",
" return states, blew_at\n",
"\n",
"fixed_states, blew = drive(use_line_search=False, STEP=2.0)\n",
"ls_states, _ = drive(use_line_search=True)\n",
"print(\"fixed step (no line search): blew up at frame\", blew, \"of\", len(offsets))\n",
"print(\"line search: completed all\", len(ls_states), \"frames\")"
]
},
{
"cell_type": "markdown",
"id": "187991e3",
"metadata": {},
"source": [
"## Without line search: the solver explodes (fixed step = 2x the Newton step)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "707d8387",
"metadata": {
"execution": {
"iopub.execute_input": "2026-06-09T05:22:33.319628Z",
"iopub.status.busy": "2026-06-09T05:22:33.319573Z",
"iopub.status.idle": "2026-06-09T05:22:33.481412Z",
"shell.execute_reply": "2026-06-09T05:22:33.481166Z"
}
},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
""
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"lims = ((-1.3, 2.2), (-1.6, 1.6))\n",
"fig, anim = utils.animate_mesh(fixed_states, T, lims=lims,\n",
" title=\"Fixed step, NO line search -> explodes\", pin_pts=X[pin_idx],\n",
" handle_traj=[s[right_idx] for s in fixed_states], fps=12)\n",
"utils.show_video(fig, anim, \"media/06_no_line_search.mp4\", fps=12)"
]
},
{
"cell_type": "markdown",
"id": "d61ec2a8",
"metadata": {},
"source": [
"## With line search: smooth and stable"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "905b8958",
"metadata": {
"execution": {
"iopub.execute_input": "2026-06-09T05:22:33.482496Z",
"iopub.status.busy": "2026-06-09T05:22:33.482427Z",
"iopub.status.idle": "2026-06-09T05:22:34.082428Z",
"shell.execute_reply": "2026-06-09T05:22:34.082178Z"
}
},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
""
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"fig, anim = utils.animate_mesh(ls_states, T, lims=lims,\n",
" title=\"Backtracking line search -> stable\", pin_pts=X[pin_idx],\n",
" handle_traj=[s[right_idx] for s in ls_states], fps=20)\n",
"utils.show_video(fig, anim, \"media/06_line_search.mp4\", fps=20)"
]
},
{
"cell_type": "markdown",
"id": "6ba8ccb1",
"metadata": {},
"source": [
"### Takeaways\n",
"* A good Newton *direction* isn't enough — an unchecked **step length** can\n",
" diverge.\n",
"* **Backtracking line search** shrinks the step until the energy decreases.\n",
"* Cheap insurance: a few extra energy evaluations buy a solver that won't blow up."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.13.9"
}
},
"nbformat": 4,
"nbformat_minor": 5
}