{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This notebook illustrates the computations of Lecture 3 of MCS 471."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 1. The Bisection Method"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We introduce the bisection method with some plots."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"\n",
"(process:52560): GLib-GIO-WARNING **: 17:46:35.620: 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": [
"plot(cos, 0.0:0.01:3.0, label=\"cos\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we extend the plot, using ``plot!``. We look for a root in the interval $[\\pi/4, 2 \\pi/3]$."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot!([pi/4], [cos(pi/4)], line=(:sticks), label=\"cos(pi/4)\")\n",
"plot!([2*pi/3], [cos(2*pi/3)], line=(:sticks), label=\"cos(2*pi/3)\")\n",
"plot!([pi/4], [cos(pi/4)], line=(:scatter, :blue), label=\"\")\n",
"plot!([2*pi/3], [cos(2*pi/3)], line=(:scatter, :blue), label=\"\")\n",
"plot!([pi/2], [0], line=(:scatter, :blue), label=\"\")\n",
"plot!([2*pi/3], [0], line=(:scatter, :green), label=\"\")\n",
"plot!([pi/4], [0], line=(:scatter, :red), label=\"\")\n",
"annotate!(pi/2, 0.1, text(\"r\"))\n",
"annotate!(2*pi/3, 0.1, text(\"b\"))\n",
"annotate!(pi/4, -0.1, text(\"a\"))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Because $\\cos(\\pi/4) \\cos(2 \\pi/3) < 0$ and $\\cos()$ is a nice continuous function, we know there is a root inside $[\\pi/4, 2 \\pi/3]$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To approximate the root, we repeatedly bisect the interval $[a, b]$, cutting the interval in half each time and keeping the invariant $f(a) f(b) < 0$ when cutting the interval."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"using Printf # to format numbers"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"bisect"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"\"\"\"\n",
"Applies the bisection method to the function f on [a,b],\n",
"f is assumed to be continuous and f(a)*f(b) < 0.\n",
"Stops when |f(a)| < eps or |f(b)| < eps or |b-a| < eps.\n",
"Returns the left and right bound of the final interval\n",
"enclosing the root, and a Boolean reporting failure.\n",
"Failure is reported when the accuracy requirement\n",
"is not satisfied in N steps; otherwise fail = 0 on return.\n",
"\n",
"Example : \n",
"\n",
" (a,b,fail) = bisect(cos,pi/4,2*pi/3,1.0e-4,100)\n",
"\"\"\"\n",
"function bisect(f::Function,\n",
" a::Float64,b::Float64,eps::Float64,N::Int)\n",
" println(\"running the bisection method...\")\n",
" println(\"-------------------------------------------\")\n",
" println(\" a b m |f(m)| |b-a| \")\n",
" println(\"-------------------------------------------\")\n",
" for i = 1:N \n",
" m = (a+b)/2\n",
" fm = f(m)\n",
" if fm*f(a) < 0\n",
" b = m\n",
" else\n",
" a = m\n",
" end\n",
" stra = @sprintf(\"%4.3f\", a)\n",
" strb = @sprintf(\"%4.3f\", b)\n",
" strm = @sprintf(\"%4.3f\", m)\n",
" strafm = @sprintf(\"%.2e\", abs(fm))\n",
" strabma = @sprintf(\"%.2e\", abs(b-a))\n",
" println(\" $stra $strb $strm $strafm $strabma\")\n",
" if (abs(fm) < eps) | ((b-a) < eps)\n",
" fail = false;\n",
" stri = string(i)\n",
" println(\"succeeded after $stri steps\")\n",
" return (a, b, fail)\n",
" end\n",
" end\n",
" strN = string(N)\n",
" println(\"failed requirements after $strN steps\")\n",
" fail = true\n",
" return (a, b, fail)\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let us run the example."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"running the bisection method...\n",
"-------------------------------------------\n",
" a b m |f(m)| |b-a| \n",
"-------------------------------------------\n",
" 1.440 2.094 1.440 1.31e-01 6.54e-01\n",
" 1.440 1.767 1.767 1.95e-01 3.27e-01\n",
" 1.440 1.604 1.604 3.27e-02 1.64e-01\n",
" 1.522 1.604 1.522 4.91e-02 8.18e-02\n",
" 1.563 1.604 1.563 8.18e-03 4.09e-02\n",
" 1.563 1.583 1.583 1.23e-02 2.05e-02\n",
" 1.563 1.573 1.573 2.05e-03 1.02e-02\n",
" 1.568 1.573 1.568 3.07e-03 5.11e-03\n",
" 1.570 1.573 1.570 5.11e-04 2.56e-03\n",
" 1.570 1.572 1.572 7.67e-04 1.28e-03\n",
" 1.570 1.571 1.571 1.28e-04 6.39e-04\n",
" 1.571 1.571 1.571 1.92e-04 3.20e-04\n",
" 1.571 1.571 1.571 3.20e-05 1.60e-04\n",
"succeeded after 13 steps\n"
]
},
{
"data": {
"text/plain": [
"(1.5707643688618154, 1.57092415852722, false)"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"(a,b,fail) = bisect(cos,pi/4,2*pi/3,1.0e-4,100)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"It is important to understand why 13 steps were needed, given the accuracy requirment of 4 decimal places, starting with the error bound of $|b - a|/2$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 2. Cobweb Diagrams"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To illustrate fixed-point iterations, we make some plots."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"using Plots"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Consider the function below to make a cobweb diagram."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"cobweb"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"\"\"\"\n",
"Applies the fixed-point iteration to a given function g\n",
"and makes a cobweb plot.\n",
"\n",
"ON ENTRY :\n",
" g a function in one variable\n",
" x0 initial guess for the fixed-point iteration\n",
" maxit upper bound on the number of iterations\n",
" tol tolerance on the abs(g(x) - x) where x is\n",
" the current approximation for the fixed point\n",
" a left bound for the range of the plot\n",
" b right bound for the range of the plot\n",
"\n",
"ON RETURN :\n",
" x the current approximation for the fixed point\n",
" numit the number of iterations done\n",
" fail true if the accuracy requirement was not met,\n",
" false otherwise;\n",
" plt a plot of the cobweb diagram.\n",
"\n",
"EXAMPLE :\n",
" g(x) = 1-x^3\n",
" (x, numit, fail, plt) = cobweb(g,0.5,10,1.0e-4,-0.1,1.1)\n",
"\"\"\"\n",
"function cobweb(g::Function,x0::Float64,maxit::Int,tol::Float64,\n",
" a::Float64,b::Float64)\n",
" r = a:0.01:b\n",
" plt = plot(g, r, label=\"g\")\n",
" ylims!((a, b))\n",
" diagonal(x) = x\n",
" plot!(diagonal, r, label=\"\")\n",
" println(\"running a fixed-point iteration ...\")\n",
" strit = @sprintf(\"%3d\", 0)\n",
" strx0 = @sprintf(\"%.4e\", x0)\n",
" println(\"$strit : $strx0\")\n",
" xprevious = x0\n",
" xnext = xprevious\n",
" for i=1:maxit\n",
" xnext = g(xprevious)\n",
" if i == 1\n",
" plot!([xprevious], [xnext], line=(:sticks, :black), label=\"\")\n",
" else\n",
" if xprevious < xnext\n",
" rng = xprevious:0.01:xnext; xr = [xprevious for i in rng]\n",
" plot!(xr, rng, line=(:black), label=\"\")\n",
" else\n",
" rng = xnext:0.01:xprevious; xr = [xprevious for i in rng]\n",
" plot!(xr, rng, line=(:black), label=\"\")\n",
" end\n",
" end\n",
" if xprevious < xnext\n",
" rng = xprevious:0.01:xnext; yr = [xnext for i in rng]\n",
" plot!(rng, yr, line=(:black), label=\"\")\n",
" else\n",
" rng = xnext:0.01:xprevious; yr = [xnext for i in rng]\n",
" plot!(rng, yr, line=(:black), label=\"\")\n",
" end\n",
" strit = @sprintf(\"%3d\", i)\n",
" strxi = @sprintf(\"%.4e\", xnext)\n",
" error = abs(xnext - xprevious)\n",
" strerr = @sprintf(\"%.2e\", error)\n",
" println(\"$strit : $strxi : $strerr\" )\n",
" if error < tol\n",
" return (xnext, i, false, plt)\n",
" end\n",
" xprevious = xnext\n",
" end\n",
" return (xnext, maxit, true, plt)\n",
"end\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let us do one step of the cobweb function ..."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"running a fixed-point iteration ...\n",
" 0 : 5.0000e-01\n",
" 1 : 8.7500e-01 : 3.75e-01\n"
]
}
],
"source": [
"g1(x) = 1 - x^3\n",
"(endpt, numit, fail, cobweb1) = cobweb(g1, 0.5, 1, 1.0e-8, -0.1, 1.1);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"and now look carefully at the plot"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"display(cobweb1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"What if we do two steps?"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"running a fixed-point iteration ...\n",
" 0 : 5.0000e-01\n",
" 1 : 8.7500e-01 : 3.75e-01\n",
" 2 : 3.3008e-01 : 5.45e-01\n"
]
}
],
"source": [
"(endpt, numit, fail, cobweb2) = cobweb(g1, 0.5, 2, 1.0e-8, -0.1, 1.1);"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"display(cobweb2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let us do 10 steps ..."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"running a fixed-point iteration ...\n",
" 0 : 5.0000e-01\n",
" 1 : 8.7500e-01 : 3.75e-01\n",
" 2 : 3.3008e-01 : 5.45e-01\n",
" 3 : 9.6404e-01 : 6.34e-01\n",
" 4 : 1.0405e-01 : 8.60e-01\n",
" 5 : 9.9887e-01 : 8.95e-01\n",
" 6 : 3.3761e-03 : 9.95e-01\n",
" 7 : 1.0000e+00 : 9.97e-01\n",
" 8 : 1.1544e-07 : 1.00e+00\n",
" 9 : 1.0000e+00 : 1.00e+00\n",
" 10 : 0.0000e+00 : 1.00e+00\n"
]
}
],
"source": [
"(endpt, numit, fail, cobweb10) = cobweb(g1, 0.5, 10, 1.0e-8, -0.1, 1.1);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"What was already noticeable in the second step should now be obvious..."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"display(cobweb10)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The shape of the picture justifies the name cobweb diagram.\n",
"*Observe that we are NOT converging to the fixed point!*"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"running a fixed-point iteration ...\n",
" 0 : 5.0000e-01\n",
" 1 : 7.9370e-01 : 2.94e-01\n",
" 2 : 5.9088e-01 : 2.03e-01\n",
" 3 : 7.4236e-01 : 1.51e-01\n",
" 4 : 6.3631e-01 : 1.06e-01\n",
" 5 : 7.1380e-01 : 7.75e-02\n",
" 6 : 6.5901e-01 : 5.48e-02\n",
" 7 : 6.9863e-01 : 3.96e-02\n",
" 8 : 6.7045e-01 : 2.82e-02\n",
" 9 : 6.9073e-01 : 2.03e-02\n",
" 10 : 6.7626e-01 : 1.45e-02\n"
]
},
{
"data": {
"text/plain": [
"(0.6762589249268274, 10, true, Plot{Plots.GRBackend() n=22})"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"g2(x) = (1 - x)^(1/3)\n",
"(endpt, numit, fail, cobwebg2) = cobweb(g2, 0.5, 10, 1.0e-8, 0.1, 0.9)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We observe that there is convergence, although the convergence is *really slow!* It makes for a good plot though."
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"display(cobwebg2)"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"running a fixed-point iteration ...\n",
" 0 : 5.0000e-01\n",
" 1 : 7.1429e-01 : 2.14e-01\n",
" 2 : 6.8318e-01 : 3.11e-02\n",
" 3 : 6.8233e-01 : 8.51e-04\n",
" 4 : 6.8233e-01 : 6.19e-07\n",
" 5 : 6.8233e-01 : 3.28e-13\n"
]
},
{
"data": {
"text/plain": [
"(0.6823278038280193, 5, false, Plot{Plots.GRBackend() n=12})"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"g3(x) = (1 + 2*x^3)/(1 + 3*x^2)\n",
"(endpt, numit, fail, cobwebg3) = cobweb(g3, 0.5, 10, 1.0e-8, -0.1, 1.1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Although we allowed 10 steps, the iterations stopped after five steps, reporting an error of ``3.28e-13``, which is really tiny. The plot reflects this very fast convergence."
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"display(cobwebg3)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To understand the convergence (or divergence), we consider the derivative of the fixed-point function at the fixed point."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.7.2",
"language": "julia",
"name": "julia-1.7"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.7.2"
}
},
"nbformat": 4,
"nbformat_minor": 2
}