{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# \n", "\n", "# Newton’s method\n", "\n", "We get started by loading our package that brings in plotting and other\n", "features, including those provided by the `Roots` package:" ], "id": "b8cbbdd1-de64-4b59-aba0-1d32bee0b835" }, { "cell_type": "code", "execution_count": 0, "metadata": {}, "outputs": [], "source": [ "using MTH229\n", "using Plots\n", "plotly()" ], "id": "2" }, { "cell_type": "markdown", "metadata": {}, "source": [ "------------------------------------------------------------------------\n", "\n", "### Quick background\n", "\n", "Read about this material here: [Newton’s\n", "Method](http://mth229.github.io/newton.html).\n", "\n", "Symbolic math - as is done behind the scenes at the Wolfram alpha web\n", "site - is pretty nice. For so many problems it can easily do what is\n", "tedious work. However, for some questions, only numeric solutions are\n", "possible. For example, there is no general formula to solve a fifth\n", "order polynomial the way there is a quadratic formula for solving\n", "quadratic polynomials. Even an innocuous polynomial like\n", "$f(x) = x^5 - x - 1$ has no easy algebraic solution for is one root.\n", "\n", "A graph shows what looks like just one answer between $1$ and $2$,\n", "closer to $1$" ], "id": "81a37c8c-7213-40f1-8caa-27248cf27243" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/html": [ "" ] } } ], "source": [ "f(x) = x^5 - x - 1\n", "plot(f, -2, 2)\n", "plot!(zero, -2, 2)" ], "id": "6" }, { "cell_type": "markdown", "metadata": {}, "source": [ "From this, we can call the bisection method to identify the value:" ], "id": "c5d1bcc3-34eb-4b24-a0ce-5bee8988d05e" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/plain": [ "1.1673039782614187" ] } } ], "source": [ "fzero(f, -2, 2)" ], "id": "8" }, { "cell_type": "markdown", "metadata": {}, "source": [ "We’ve seen the bisection method previously to find a root, but this is\n", "somewhat cumbersome to use as it needs a\n", "[bracketing](https://en.wikipedia.org/wiki/Bisection_method#The_method)\n", "interval to begin. Moreover, bisection can be computationally slow (by\n", "comparison).\n", "\n", "Here we discuss Newton’s method. Like the bisection method it is an\n", "*iterative algorithm*. However instead of identifying a bracketing\n", "interval, we only need to identify a **nearby** *initial* guess, $x_0$.\n", "\n", "Starting with $x_0$ the algorithm to produce $x_1$ is easy to describe:\n", "\n", "- form the tangent line at $(x_0, f(x_0))$.\n", "\n", "- let $x_1$ be the intersection point of this tangent line with the\n", " $x$ axis.\n", "\n", "If we can go from $x_0$ to $x_1$ we can repeat the update step to get\n", "$x_2$ and then $x_3$, …\n", "\n", "Graphically, this figure illustrates the process until the guesses get\n", "so close they crowd the figure. (Ignore the code to produce the figure.)" ], "id": "9ff9a5a2-a7a9-4e30-94bc-40771140ded1" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/html": [ "" ] } } ], "source": [ "f(x) = x^5 - x - 1\n", "plot(f, 0.9, 1.3; legend=false, linewidth=2)\n", "plot!(zero)\n", "newton_plot!(f, 1.0; annotate_steps=3)" ], "id": "10" }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the figure, the sequence of guesses can be seen, basically $x_0=1$,\n", "$x_1=1.25$, $x_2=1.178\\dots$, …, $1.16730\\dots$.\n", "\n", "To find these values numerically, we first need an algebraic\n", "representation. For this problem, we can describe the tangent line’s\n", "slope by *either* $f'(x_0)$ *or* by using “rise over run”:\n", "\n", "$$\n", "f'(x_0) = \\frac{f(x_0) - 0}{x_0 - x_1}.\n", "$$\n", "\n", "Solving for $x_0$, this yields the update formula:\n", "\n", "$$\n", "x_1 = x_0 - f(x_0) / f'(x_0).\n", "$$\n", "\n", "That is, the new guess shifts the old guess by an increment\n", "$f(x_0)/f'(x_0)$.\n", "\n", "In `Julia`, we can do one step with:" ], "id": "5d51810b-7bf1-4c05-ad69-8b33cf403e95" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/plain": [ "1.25" ] } } ], "source": [ "f(x) = x^5 - x - 1\n", "x = 1\n", "x = x - f(x) / f'(x)" ], "id": "12" }, { "cell_type": "markdown", "metadata": {}, "source": [ "(We don’t use indexing, but rather update our binding for the `x`\n", "variable.)\n", "\n", "Is `x` close to being the zero? We don’t know the actual zero - we are\n", "trying to approximate it - but we do know the function’s value at the\n", "actual zero. For this new guess the function value is" ], "id": "ab096154-8613-4b34-b1ee-7f5b80f79456" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/plain": [ "0.8017578125" ] } } ], "source": [ "f(x)" ], "id": "14" }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is much closer to $0$ than $f(1)$, the value at our initial guess,\n", "but not nearly as close as we can get using Newton’s method. We just\n", "need to **iterate** - run a few more steps.\n", "\n", "We do another step just by running the last line. For example, we run 5\n", "more steps by copying and pasting the same expression:" ], "id": "1d881930-fe59-4af3-902f-db26af41d40c" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/plain": [ "1.1673039782614187" ] } } ], "source": [ "x = x - f(x) / f'(x)\n", "x = x - f(x) / f'(x)\n", "x = x - f(x) / f'(x)\n", "x = x - f(x) / f'(x)\n", "x = x - f(x) / f'(x)" ], "id": "16" }, { "cell_type": "markdown", "metadata": {}, "source": [ "The value of `x` updates. But is it getting closer to a *zero*? If so,\n", "then $f(x)$ should be close to zero. We can see both values with:" ], "id": "db2972ce-5948-473e-8441-5b3bbeb93a09" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/plain": [ "(1.1673039782614187, 6.661338147750939e-16)" ] } } ], "source": [ "x, f(x)" ], "id": "18" }, { "cell_type": "markdown", "metadata": {}, "source": [ "This shows $f(x)$ is not exactly $0.0$ but it is as close as we can get.\n", "Repeating the algorithm does not change the value of `x`. (On a\n", "computer, floating point issues creep in when values are close to 0, and\n", "these prevent values being mathematically exact.) As we can’t improve,\n", "we stop. Our value of `x` is an *approximate* zero and `f(x)` is within\n", "machine tolerance of being `0`.\n", "\n", "How do we know how to stop? When the algorithm works, we will stop when\n", "the `x` value *basically* stops updating, as `f(x)` is basically `0`.\n", "However, the algorithm need not work, so any implementation must keep\n", "track of how many steps are taken and stop when this gets out of hand.\n", "\n", "For convenience, the `newton` function in the `MTH229` package (using\n", "the `Roots` package) will iterate until convergence. If we pass in the\n", "optional argument `verbose=true` we will see the sequence of steps.\n", "\n", "For example, for $f(x) = x^3 - 2x - 5$, a function that Newton himself\n", "considered, a solution near $2$, is found with:" ], "id": "7b3c216c-6f6c-4a45-a645-7f099610f30c" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/plain": [ "2.0945514815423265" ] } } ], "source": [ "x = 2\n", "f(x) = x^3 - 2x -5\n", "xstar = newton(f, 2)" ], "id": "20" }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can see the approximate zero and the function value, as follows:" ], "id": "e2465c81-0400-4867-91ef-331d783ef7dd" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/plain": [ "(2.0945514815423265, -8.881784197001252e-16)" ] } } ], "source": [ "xstar, f(xstar)" ], "id": "22" }, { "cell_type": "markdown", "metadata": {}, "source": [ "------------------------------------------------------------------------\n", "\n", "### Using fzero from the Roots package\n", "\n", "As mentioned, the `newton` function in the `Roots` package implements\n", "Newton’s method. The `Roots` package also provides the `fzero` function\n", "for finding roots. This function will:\n", "\n", "- use bisection when called as `fzero(f, a, b)` with the zero\n", " **between** the bracketing interval `[a,b]`\n", "- use a method like Newton’s method when called as `fzero(f, c)` with\n", " the zero **near** `c`.\n", "\n", "For example:" ], "id": "45d701d1-5d33-4423-818b-916d9c699548" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "output_type": "display_data", "metadata": {}, "data": { "text/plain": [ "3.141592653589793" ] } } ], "source": [ "fzero(sin, 3) # start with initial guess of 3, returns 3.141592653589793" ], "id": "24" }, { "cell_type": "markdown", "metadata": {}, "source": [ "The utility of this function is that it does not require a derivative to\n", "be taken and it is a little less sensitive than Newton’s method to the\n", "initial guess. The use of `fzero` is recommended: just plot to find a\n", "nearby guess and use `fzero(f,c)` to improve this guess.\n", "\n", "### When Newton’s method fails\n", "\n", "The error in the $n$th step using Newton’s method at a simple zero\n", "follows a formula: $|e_{n+1}| \\leq (1/2) |f''(a)/f'(b)| \\cdot |e_n|^2$,\n", "for some $a$ and $b$. Generally this ensures that the error at step\n", "$n+1$ is smaller than the error at step $n$ squared. But this can fail\n", "due to various cases:\n", "\n", "- the initial guess is not *near* the zero,\n", "\n", "- the derivative, $|f'(x)|$, is too small, or\n", "\n", "- the second derivative, $|f''(x)|$, is too big, or possibly\n", " undefined.\n", "\n", "### Quadratic convergence\n", "\n", "When Newton’s method converges to a *simple zero* it is said to have\n", "*quadratic convergence*. A simple zero is one with multiplicity 1 and\n", "quadratic convergence says basically that the error at the $i+1$st step\n", "is like the error for $i$th step squared. In particular, if the error is\n", "like $10^{-3}$ on one step, it will be like $10^{-6}$, then $10^{-12}$\n", "then $10^{-24}$ on subsequent steps. (Which is typically beyond the\n", "limit of a floating point approximation.) This is why one can *usually*\n", "take just 5, or so, steps to get to an answer.\n", "\n", "Not so for multiple roots and some simple roots.\n", "\n", "------------------------------------------------------------------------" ], "id": "3eb7ec93-b5f0-4b4d-99aa-3944e8d560b3" }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# Your commands go here\n" ], "id": "26" } ], "nbformat": 4, "nbformat_minor": 5, "metadata": { "kernel_info": { "name": "julia" }, "kernelspec": { "name": "julia", "display_name": "Julia", "language": "julia" }, "language_info": { "name": "julia", "codemirror_mode": "julia", "version": "1.10.0" } } }