{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Newton's Method"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Given an approximation $x_0$ for a root of $f$, we compute the next approximation $x_1$ as the point where the tangent line of $f$ at $x_0$ intersects the horizontal axes."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"\n",
"(process:65500): GLib-GIO-WARNING **: 11:56:44.657: Unexpectedly, UWP app `Evernote.Evernote_10.43.7.0_x64__q4d96b2w5wcc2' (AUMId `Evernote.Evernote_q4d96b2w5wcc2!Evernote') supports 1 extensions but has no verbs\n"
]
}
],
"source": [
"using Plots"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f(x) = x*(x-1.5)\n",
"plot(f, -0.5:0.1:3.5, label=\"f\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let $x_0 = 3$. The derivative of $f = x^2 - 1.5 x$ is $f' = 2 x - 1.5$."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4.5"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df(x) = 2*x - 1.5\n",
"slope = df(3.0)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The equation for the tangent line at $(3, f(3))$ is $y - f(3) = f'(3)(x - 3)$. To find $x_1$ solve this equation for $y = 0$:\n",
"\n",
"$$x_1 = 3 - \\frac{f(3)}{f'(3)}$$."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2.0"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x1 = 3 - f(3.0)/slope"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we can complete the plot."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot!([3], [f(3)], line=(:sticks, \"black\"), label=\"f(x0)\")\n",
"plot!([3], [0], line=(:scatter, :red), label=\"\")\n",
"plot!([3], [f(3)], line=(:scatter, :red), label=\"\")\n",
"plot!([2], [0], line=(:scatter, :red), label=\"\")\n",
"plot!([2, 3], [0, f(3)], line=(:red), label=\"tangent\")\n",
"annotate!(3, -0.4, text(\"x0\"))\n",
"annotate!(2, -0.4, text(\"x1\"))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# The Newton Operator"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Newton's method is a fixed-point iteration of the form $x = g(x)$, defined by\n",
"\n",
"$$\n",
" g(x) = x - \\frac{f(x)}{f'(x)}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We will use SymPy to compute derivatives of expressions in x."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$x$"
],
"text/plain": [
"x"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using SymPy\n",
"x = Sym(\"x\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Strings are convenient to store and represent expressions. The string representations of expression are the input and output of the following function."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"SymPyDerivative"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"\"\"\"\n",
"Returns the string representation of the derivative of the expression in x,\n",
"given in the string s.\n",
"\"\"\"\n",
"function SymPyDerivative(s::String)\n",
" evaluated = sympify(s) # evaluates the string as a SymPy expression in x\n",
" derivative = diff(evaluated, x) # the derivative is a SymPy expression in x\n",
" return string(derivative) # returns the string representation of the derivative\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"\"2*x*exp(x^2)*sin(x) + exp(x^2)*cos(x)\""
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"example = SymPyDerivative(\"exp(x^2)*sin(x)\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With ``sympify`` we turn a string into a SymPy expression."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$2 x e^{x^{2}} \\sin{\\left(x \\right)} + e^{x^{2}} \\cos{\\left(x \\right)}$"
],
"text/plain": [
" ⎛ 2⎞ ⎛ 2⎞ \n",
" ⎝x ⎠ ⎝x ⎠ \n",
"2⋅x⋅ℯ ⋅sin(x) + ℯ ⋅cos(x)"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"symfun = sympify(example)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With ``lambdify`` we turn a SymPy expression into a Julia function. "
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"#118 (generic function with 1 method)"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"exafun = lambdify(symfun)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-1160.948672607985"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"exafun(3.0)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For the iterative algorithm, instead of the entire\n",
"$$\n",
" g(x) = x - \\frac{f(x)}{f'(x)}\n",
"$$\n",
"we will define the update function\n",
"$$\n",
" \\Delta x = - \\frac{f(x)}{f'(x)}\n",
"$$\n",
"used to compute $x_1 = x_0 + \\Delta x$."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"NewtonUpdate"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"\"\"\"\n",
"Given a string representation of an expression in x,\n",
"returns the Julia function for the update -f(x)/f'(x),\n",
"where f(x) is the function defined by the input string.\n",
"\"\"\"\n",
"function NewtonUpdate(f::String, verbose::Bool=true)\n",
" sdf = SymPyDerivative(f)\n",
" sdx = string(\"-(\", f, \")/(\", sdf, \")\")\n",
" if verbose\n",
" println(\"update : $sdx\")\n",
" end\n",
" edx = sympify(sdx)\n",
" return lambdify(edx)\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"update : -(x^3 + x - 1)/(3*x^2 + 1)\n"
]
},
{
"data": {
"text/plain": [
"-0.25"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"dx = NewtonUpdate(\"x^3 + x - 1\")\n",
"dx(1.0)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Running Newton's Method"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we are ready to define a function to run Newton's method."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Newton"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using Printf # to format floating-point numbers\n",
"\n",
"\"\"\"\n",
"Runs Newton's method for the function f(x) = 0.\n",
"\n",
"ON ENTRY :\n",
"\n",
" f a function in one variable\n",
"\n",
" x0 initial approximation for the root\n",
"\n",
" N maximum number of iterations\n",
"\n",
" dxtol tolerance on the magnitude of the update\n",
"\n",
" fxtol tolerance on the residual\n",
"\n",
"ON RETURN :\n",
"\n",
" (root, |dx|, |f(root)|, failure)\n",
" a tuple with the computed approximation for the root\n",
" and a boolean which is true if the method did not converge.\n",
"\"\"\"\n",
"function Newton(f::Function, update::Function, x0::Float64,\n",
" N::Int64=6, dxtol::Float64=1.0e-8, fxtol::Float64=1.0e-8)\n",
" myroot = x0; dx = 1; y = 1.0\n",
" stri = @sprintf(\"%3d\", 0)\n",
" strx = @sprintf(\"%.16e\", myroot)\n",
" println(\"step root |dx| |f(x)|\")\n",
" println(\"$stri $strx\")\n",
" for i=1:N\n",
" y = abs(Float64(f(myroot)))\n",
" dx = Float64(update(myroot))\n",
" myroot = myroot + dx\n",
" stri = @sprintf(\"%3d\", i)\n",
" strx = @sprintf(\"%.16e\", myroot)\n",
" stry = @sprintf(\"%.2e\", y)\n",
" strdx = @sprintf(\"%.2e\", abs(dx))\n",
" println(\"$stri $strx $strdx $stry\")\n",
" if((abs(dx) < dxtol) | (y < fxtol))\n",
" return (myroot, abs(dx), y, false)\n",
" end\n",
" end\n",
" return (myroot, abs(dx), y, true)\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"update : -(x^3 + x - 1)/(3*x^2 + 1)\n",
"step root |dx| |f(x)|\n",
" 0 5.0000000000000000e-01\n",
" 1 7.1428571428571430e-01 2.14e-01 3.75e-01\n",
" 2 6.8317972350230416e-01 3.11e-02 7.87e-02\n",
" 3 6.8232842330457821e-01 8.51e-04 2.04e-03\n",
" 4 6.8232780382834712e-01 6.19e-07 1.48e-06\n",
" 5 6.8232780382801939e-01 3.28e-13 7.86e-13\n"
]
},
{
"data": {
"text/plain": [
"(0.6823278038280194, 3.2777958153976343e-13, 7.855938122247608e-13, false)"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sf = \"x^3 + x - 1\"\n",
"dx = NewtonUpdate(sf)\n",
"fx = lambdify(sympify(sf))\n",
"Newton(fx,dx,0.5)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.8.0",
"language": "julia",
"name": "julia-1.8"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.8.0"
}
},
"nbformat": 4,
"nbformat_minor": 4
}