{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Newton method" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Consider the function\n", "$$\n", "f(x) = \\exp(x) - \\frac{3}{2} - \\arctan(x)\n", "$$" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%config InlineBackend.figure_format = 'svg'\n", "from numpy import linspace,abs\n", "from sympy import Symbol,exp,atan,diff,lambdify\n", "from matplotlib.pyplot import figure,subplot,plot,grid,title" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We define the function using sympy and automatically compute its derivative. Then we make functions out of these." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "f = exp(x) - atan(x) - 1.5\n", "df= exp(x) - 1/(x**2 + 1)\n" ] }, { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", " \n", " 2023-01-17T12:27:04.339117\n", " image/svg+xml\n", " \n", " \n", " Matplotlib v3.6.2, https://matplotlib.org/\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "x = Symbol('x')\n", "\n", "# Define the function\n", "f = exp(x)-3.0/2.0-atan(x)\n", "\n", "# Compute its derivative\n", "df = diff(f,x)\n", "print(\"f = \", f)\n", "print(\"df= \", df)\n", "\n", "# Make functions\n", "F = lambdify(x, f, modules=['numpy'])\n", "DF = lambdify(x, df, modules=['numpy'])\n", "\n", "# Plot the functions\n", "X = linspace(-1,1,100)\n", "figure(figsize=(9,4))\n", "subplot(1,2,1)\n", "plot(X,F(X))\n", "grid(True)\n", "title('Function')\n", "subplot(1,2,2)\n", "plot(X,DF(X))\n", "grid(True)\n", "title('Derivative');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we implement the Newton method. We have to specify the maximum number of iterations, an initial guess and a tolerance to decide when to stop." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 0 8.71059792151655e-01 3.71059792151655e-01 1.72847808769943e-01\n", " 1 7.76133043351671e-01 9.49267487999836e-02 1.30353403123729e-02\n", " 2 7.67717620302437e-01 8.41542304923425e-03 9.81774270505387e-05\n", " 3 7.67653269950858e-01 6.43503515784806e-05 5.71995606435394e-09\n", " 4 7.67653266201279e-01 3.74957960683442e-09 2.22044604925031e-16\n", " 5 7.67653266201279e-01 1.45556000556775e-16 0.00000000000000e+00\n" ] } ], "source": [ "M = 20 # maximum iterations\n", "x = 0.5 # initial guess\n", "eps = 1e-14 # relative tolerance on root\n", "\n", "f = F(x)\n", "for i in range(M):\n", " df = DF(x)\n", " dx = -f/df\n", " x = x + dx\n", " e = abs(dx)\n", " f = F(x)\n", " print(\"%6d %22.14e %22.14e %22.14e\" % (i,x,e,abs(f)))\n", " if e < eps * abs(x):\n", " break" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercise\n", "\n", "Modify the above Newton method implemention to use `while` loop." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.10.8" } }, "nbformat": 4, "nbformat_minor": 4 }