Files
numethods/tutorials/tutorial4_root_finding.ipynb

343 lines
14 KiB
Plaintext
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

{
"cells": [
{
"cell_type": "markdown",
"id": "3c51341a",
"metadata": {},
"source": [
"# Tutorial 4: Root-Finding Methods\n",
"---\n",
"In this tutorial, we study classical **root-finding algorithms** for nonlinear equations. We will:\n",
"\n",
"- Define the root-finding problem mathematically\n",
"- Derive several algorithms (bisection, fixed-point, Newton, secant)\n",
"- Discuss convergence conditions and error behavior\n",
"- Compare methods with worked examples using the `numethods` package.\n"
]
},
{
"cell_type": "markdown",
"id": "2999aa4f",
"metadata": {},
"source": [
"## 1. Problem Setup and Notation\n",
"\n",
"We seek to solve a nonlinear scalar equation:\n",
"$$\n",
"f(x) = 0, \\quad f: \\mathbb{R} \\to \\mathbb{R}.\n",
"$$\n",
"with a continuously differentiable function.\n",
"\n",
"\n",
"### Root, residual, and error\n",
"- A **root** $r$ satisfies $f(r)=0$.\n",
"- **Absolute error:** $(e_k = |x_k - x^\\star|)$.\n",
"- **Residual:** $(r_k = |f(x_k)|)$. \n",
"Note that small residual does not always imply small error.\n",
"\n",
"### Multiplicity\n",
"A root $r$ has **multiplicity** $m$ if\n",
"$$\n",
"f(r) = f'(r) = \\dots = f^{(m-1)}(r) = 0, \\quad f^{(m)}(r) \\neq 0.\n",
"$$\n",
"If $x^\\star$ satisfies $f(x^\\star)=0$ and $f'(x^\\star)\\ne 0$, we say the root is **simple** (multiplicity 1).\n",
"\n",
"If $f'(x^\\star)=\\cdots=f^{(m-1)}(x^\\star)=0$ and $f^{(m)}(x^\\star)\\ne 0$, we say the root has **multiplicity** (m).\n",
"\n",
"- **Simple roots** ($m=1$): most methods converge rapidly.\n",
"- **Multiple roots** ($m>1$): convergence often slows.\n"
]
},
{
"cell_type": "markdown",
"id": "b610cea2",
"metadata": {},
"source": [
"## 2. Bisection Method\n",
"\n",
"**Assumption (Intermediate Value Theorem):** If f is continuous on ([a,b]) and (f(a),f(b) < 0),\n",
"then there exists $x^\\star$ in (a,b) with $f(x^\\star)=0$.\n",
"\n",
"- Assumes $f$ is continuous on $[a,b]$ with $f(a)f(b)<0$.\n",
"- Repeatedly bisect interval and select subinterval containing the root.\n",
"\n",
"**Iteration:**\n",
"$$\n",
"c_k = \\frac{a_k+b_k}{2}, \\quad f(c_k).\n",
"$$\n",
"\n",
"**Error bound:** interval length halves each step:\n",
"$$\n",
"|c_k-r| \\le \\frac{b-a}{2^k}.\n",
"$$\n",
"- Convergence: **linear**, guaranteed.\n"
]
},
{
"cell_type": "markdown",
"id": "3be4a510",
"metadata": {},
"source": [
"## 3. Fixed-Point Iteration\n",
"- Rewrite equation as $x=g(x)$.\n",
"- Iterate:\n",
"$$\n",
"x_{k+1} = g(x_k).\n",
"$$\n",
"\n",
"**Convergence theorem (Banach fixed-point):** If g is continuously differentiable near $(x^\\star)$ and\n",
"$$\n",
"|g'(x_\\star)| < 1,\n",
"$$\n",
"then for initial guesses $x_0$ sufficiently close to $x^\\star$, the iterates converge **linearly** to $x^\\star$ with asymptotic rate $|g'(x^\\star)|$.\n",
"\n",
"**Choice of g.** Different rearrangements yield different g's with different convergence properties.\n",
"A poor choice (with $(|g'|\\ge 1))$ can diverge.\n",
"\n",
"- Rate: linear with factor $|g'(r)|$.\n"
]
},
{
"cell_type": "markdown",
"id": "40e66a30",
"metadata": {},
"source": [
"## 4. Newtons Method\n",
"From Taylor expansion:\n",
"$$\n",
"f(x) \\approx f(x_k) + f'(x_k)(x-x_k).\n",
"$$\n",
"Set $f(x)=0$ to solve for next iterate:\n",
"$$\n",
"x_{k+1} = x_k - \\frac{f(x_k)}{f'(x_k)}.\n",
"$$\n",
"\n",
"- **Quadratic convergence** for simple roots.\n",
"- For multiple roots: drops to linear.\n",
"- Requires derivative, sensitive to initial guess.\n"
]
},
{
"cell_type": "markdown",
"id": "84888305",
"metadata": {},
"source": [
"## 5. Secant Method\n",
"- Avoids derivative by approximating slope with finite difference:\n",
"$$\n",
"x_{k+1} = x_k - f(x_k) \\frac{x_k - x_{k-1}}{f(x_k)-f(x_{k-1})}.\n",
"$$\n",
"\n",
"- Convergence order: $\\approx 1.618$ (superlinear).\n",
"- More efficient than Newton if derivative expensive.\n"
]
},
{
"cell_type": "markdown",
"id": "7ed8a0b5",
"metadata": {},
"source": [
"## 6. Stopping Criteria\n",
"We stop iteration when:\n",
"- $|f(x_k)| \\le \\varepsilon$ (residual small), or\n",
"- $|x_{k+1}-x_k| \\le \\varepsilon(1+|x_{k+1}|)$ (relative step small).\n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "4c9f93d9",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Bisection root: 1.4142135605216026\n",
"Newton root: 1.4142135623730951\n",
"Secant root: 1.414213562373095\n",
"Fixed-point failed: Fixed-point iteration did not converge\n"
]
}
],
"source": [
"from numethods import Bisection, FixedPoint, NewtonRoot, Secant, print_trace\n",
"import math\n",
"\n",
"# Example function: f(x) = x^2 - 2\n",
"f = lambda x: x**2 - 2\n",
"df = lambda x: 2*x\n",
"\n",
"# Bisection\n",
"bisect = Bisection(f, 0, 2, tol=1e-8)\n",
"root_b = bisect.solve()\n",
"print('Bisection root:', root_b)\n",
"\n",
"# Newton\n",
"newton = NewtonRoot(f, df, 1.0, tol=1e-12)\n",
"root_n = newton.solve()\n",
"print('Newton root:', root_n)\n",
"\n",
"# Secant\n",
"sec = Secant(f, 0, 2, tol=1e-12)\n",
"root_s = sec.solve()\n",
"print('Secant root:', root_s)\n",
"\n",
"# Fixed point: g(x)=sqrt(2) ~ rewriting\n",
"g = lambda x: (2/x) # not always convergent, but demonstrates\n",
"try:\n",
" fp = FixedPoint(g, 1.0, tol=1e-8)\n",
" root_fp = fp.solve()\n",
" print('Fixed-point root:', root_fp)\n",
"except Exception as e:\n",
" print('Fixed-point failed:', e)\n"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "7f8fe13e",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Newton Method Trace (x^2 - 2):\n",
" iter | x | f(x) | df(x) | x_new | error\n",
"------------------------------------------------------------------------------\n",
" 0 | 1 | -1 | 2 | 1.5 | 0.5\n",
" 1 | 1.5 | 0.25 | 3 | 1.41667 | 0.0833333\n",
" 2 | 1.41667 | 0.00694444 | 2.83333 | 1.41422 | 0.00245098\n",
" 3 | 1.41422 | 6.0073e-06 | 2.82843 | 1.41421 | 2.1239e-06\n",
" 4 | 1.41421 | 4.51061e-12 | 2.82843 | 1.41421 | 1.59472e-12\n",
"\n",
"Secant Method Trace (x^2 - 2):\n",
" iter | x0 | x1 | x2 | f(x0) | f(x1) | error\n",
"-------------------------------------------------------------------------------------------\n",
" 0 | 0 | 2 | 1 | -2 | 2 | 1\n",
" 1 | 2 | 1 | 1.33333 | 2 | -1 | 0.333333\n",
" 2 | 1 | 1.33333 | 1.42857 | -1 | -0.222222 | 0.0952381\n",
" 3 | 1.33333 | 1.42857 | 1.41379 | -0.222222 | 0.0408163 | 0.0147783\n",
" 4 | 1.42857 | 1.41379 | 1.41421 | 0.0408163 | -0.00118906 | 0.000418335\n",
" 5 | 1.41379 | 1.41421 | 1.41421 | -0.00118906 | -6.00729e-06 | 2.12421e-06\n",
" 6 | 1.41421 | 1.41421 | 1.41421 | -6.00729e-06 | 8.93146e-10 | 3.15775e-10\n",
" 7 | 1.41421 | 1.41421 | 1.41421 | 8.93146e-10 | -8.88178e-16 | 2.22045e-16\n",
"\n",
"Bisection Method Trace (x^2 - 2):\n",
" iter | a | b | c | f(a) | f(b) | f(c) | interval\n",
"--------------------------------------------------------------------------------------------------------\n",
" 0 | 0 | 2 | 1 | -2 | 2 | -1 | 2\n",
" 1 | 1 | 2 | 1.5 | -1 | 2 | 0.25 | 1\n",
" 2 | 1 | 1.5 | 1.25 | -1 | 0.25 | -0.4375 | 0.5\n",
" 3 | 1.25 | 1.5 | 1.375 | -0.4375 | 0.25 | -0.109375 | 0.25\n",
" 4 | 1.375 | 1.5 | 1.4375 | -0.109375 | 0.25 | 0.0664062 | 0.125\n",
" 5 | 1.375 | 1.4375 | 1.40625 | -0.109375 | 0.0664062 | -0.0224609 | 0.0625\n",
" 6 | 1.40625 | 1.4375 | 1.42188 | -0.0224609 | 0.0664062 | 0.0217285 | 0.03125\n",
" 7 | 1.40625 | 1.42188 | 1.41406 | -0.0224609 | 0.0217285 | -0.000427246 | 0.015625\n",
" 8 | 1.41406 | 1.42188 | 1.41797 | -0.000427246 | 0.0217285 | 0.0106354 | 0.0078125\n",
" 9 | 1.41406 | 1.41797 | 1.41602 | -0.000427246 | 0.0106354 | 0.00510025 | 0.00390625\n",
" 10 | 1.41406 | 1.41602 | 1.41504 | -0.000427246 | 0.00510025 | 0.00233555 | 0.00195312\n",
" 11 | 1.41406 | 1.41504 | 1.41455 | -0.000427246 | 0.00233555 | 0.000953913 | 0.000976562\n",
" 12 | 1.41406 | 1.41455 | 1.41431 | -0.000427246 | 0.000953913 | 0.000263274 | 0.000488281\n",
" 13 | 1.41406 | 1.41431 | 1.41418 | -0.000427246 | 0.000263274 | -8.20011e-05 | 0.000244141\n",
" 14 | 1.41418 | 1.41431 | 1.41425 | -8.20011e-05 | 0.000263274 | 9.06326e-05 | 0.00012207\n",
" 15 | 1.41418 | 1.41425 | 1.41422 | -8.20011e-05 | 9.06326e-05 | 4.31482e-06 | 6.10352e-05\n",
" 16 | 1.41418 | 1.41422 | 1.4142 | -8.20011e-05 | 4.31482e-06 | -3.88434e-05 | 3.05176e-05\n",
" 17 | 1.4142 | 1.41422 | 1.41421 | -3.88434e-05 | 4.31482e-06 | -1.72643e-05 | 1.52588e-05\n",
" 18 | 1.41421 | 1.41422 | 1.41421 | -1.72643e-05 | 4.31482e-06 | -6.47477e-06 | 7.62939e-06\n",
" 19 | 1.41421 | 1.41422 | 1.41421 | -6.47477e-06 | 4.31482e-06 | -1.07998e-06 | 3.8147e-06\n",
" 20 | 1.41421 | 1.41422 | 1.41421 | -1.07998e-06 | 4.31482e-06 | 1.61742e-06 | 1.90735e-06\n",
" 21 | 1.41421 | 1.41421 | 1.41421 | -1.07998e-06 | 1.61742e-06 | 2.68718e-07 | 9.53674e-07\n",
" 22 | 1.41421 | 1.41421 | 1.41421 | -1.07998e-06 | 2.68718e-07 | -4.05632e-07 | 4.76837e-07\n",
" 23 | 1.41421 | 1.41421 | 1.41421 | -4.05632e-07 | 2.68718e-07 | -6.84571e-08 | 2.38419e-07\n",
" 24 | 1.41421 | 1.41421 | 1.41421 | -6.84571e-08 | 2.68718e-07 | 1.0013e-07 | 1.19209e-07\n",
" 25 | 1.41421 | 1.41421 | 1.41421 | -6.84571e-08 | 1.0013e-07 | 1.58366e-08 | 5.96046e-08\n",
" 26 | 1.41421 | 1.41421 | 1.41421 | -6.84571e-08 | 1.58366e-08 | -2.63102e-08 | 2.98023e-08\n",
" 27 | 1.41421 | 1.41421 | 1.41421 | -2.63102e-08 | 1.58366e-08 | -5.23681e-09 | 1.49012e-08\n",
" 28 | 1.41421 | 1.41421 | 1.41421 | -5.23681e-09 | 1.58366e-08 | 5.2999e-09 | 7.45058e-09\n",
" 29 | 1.41421 | 1.41421 | 1.41421 | -5.23681e-09 | 5.2999e-09 | 3.15445e-11 | 3.72529e-09\n",
"\n",
"Fixed-Point Iteration Trace (x^2 - 2):\n",
" iter | x | x_new | error\n",
"----------------------------------------------------\n",
" 0 | 1 | 1.5 | 0.5\n",
" 1 | 1.5 | 1.41667 | 0.0833333\n",
" 2 | 1.41667 | 1.41422 | 0.00245098\n",
" 3 | 1.41422 | 1.41421 | 2.1239e-06\n",
" 4 | 1.41421 | 1.41421 | 1.59495e-12\n"
]
}
],
"source": [
"# Newton\n",
"steps = NewtonRoot(f, df, x0=1.0).trace()\n",
"print(\"Newton Method Trace (x^2 - 2):\")\n",
"print_trace(steps)\n",
"\n",
"# Secant\n",
"steps = Secant(f, 0, 2).trace()\n",
"print(\"\\nSecant Method Trace (x^2 - 2):\")\n",
"print_trace(steps)\n",
"\n",
"# Bisection\n",
"steps = Bisection(f, 0, 2).trace()\n",
"print(\"\\nBisection Method Trace (x^2 - 2):\")\n",
"print_trace(steps)\n",
"\n",
"# Fixed-point: solve\n",
"g = lambda x: 0.5 * (x + 2 / x)\n",
"steps = FixedPoint(g, 1.0).trace()\n",
"print(\"\\nFixed-Point Iteration Trace (x^2 - 2):\")\n",
"print_trace(steps)"
]
},
{
"cell_type": "markdown",
"id": "833b538c",
"metadata": {},
"source": [
"## 7. Comparison of Methods\n",
"| Method | Requires derivative | Convergence rate | Guarantee? |\n",
"|--------|---------------------|------------------|------------|\n",
"| Bisection | No | Linear | Yes (if sign change) |\n",
"| Fixed-Point | No | Linear | Not always |\n",
"| Newton | Yes | Quadratic | Locally (good guess) |\n",
"| Secant | No | ~1.618 (superlinear) | Locally (good guess) |\n"
]
},
{
"cell_type": "markdown",
"id": "0d372aa0",
"metadata": {},
"source": [
"## 8. Exercises\n",
"1. Apply all four methods to $f(x)=\\cos x - x$.\n",
"2. Try Newtons method on $f(x)=(x-1)^2$ and compare convergence rate with $f(x)=x^2-2$.\n",
"3. Modify Secant to stop if denominator too small.\n"
]
}
],
"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.7"
}
},
"nbformat": 4,
"nbformat_minor": 5
}