{ "cells": [ { "cell_type": "markdown", "id": "0e81a916", "metadata": {}, "source": [ "In lecture 6 of MCS 471, we consider the secant method and the method of false position (a.k.a. the regula falsi method)." ] }, { "cell_type": "markdown", "id": "a1212356", "metadata": {}, "source": [ "# 1. The Secant Method" ] }, { "cell_type": "markdown", "id": "d7731d67", "metadata": {}, "source": [ "The secant method is similar to Newton's method, but uses the secant instead of the tangent line, thus avoiding the calculation of the derivative." ] }, { "cell_type": "markdown", "id": "94f4ccf0", "metadata": {}, "source": [ "A simple implementation of the secant method is in the function below." ] }, { "cell_type": "code", "execution_count": 1, "id": "79ec6de4", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "secant" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using Printf\n", "\"\"\"\n", "Applies the secant method to a function.\n", "\n", "ON ENTRY :\n", "\n", " f a function in one variable\n", " x1,x2 two start points\n", " dxtol tolerance on the forward error\n", " fxtol tolerance on the backward error\n", " N maximum number of iterations\n", "\n", "ON RETURN : (x, absdx, absfx, nbrit, fail)\n", "\n", " x approximation for the root\n", " absdx estimated forward error\n", " absfx estimated backward error\n", " nbrit number of iterations\n", " fail true if tolerances not reached,\n", " false otherwise.\n", "\n", "EXAMPLE:\n", "\n", " (root, err, res, nit, fail) = secant(cos,pi/4,2*pi/3)\n", "\"\"\"\n", "function secant(f::Function,x1::Float64,x2::Float64,\n", " dxtol::Float64=1.0e-8,fxtol::Float64=1.0e-8,N::Int64=10)\n", " dx = 1\n", " println(\"running the secant method...\")\n", " title = \" root |dx| |f(x)|\"\n", " println(\"step : $title\")\n", " fx1 = f(x1)\n", " fx2 = f(x2)\n", " for i = 1:N \n", " dx = fx2*(x2 - x1)/(fx2 - fx1)\n", " x1 = x2 # store x2 for next step\n", " x2 = x2 - dx\n", " fx1 = fx2\n", " fx2 = f(x2)\n", " stri = @sprintf(\"%3d\", i)\n", " strx = @sprintf(\"%.16e\", x2)\n", " strdx = @sprintf(\"%.2e\", abs(dx))\n", " strfx = @sprintf(\"%.2e\", abs(fx2))\n", " println(\"$stri : $strx $strdx $strfx\")\n", " if((abs(fx2) < dxtol) | (abs(dx) < fxtol))\n", " stri = string(i)\n", " println(\"succeeded after $stri steps\")\n", " return (x2, abs(dx), abs(fx2), i, false)\n", " end\n", " end\n", " strN = string(N)\n", " println(\"failed requirements after $strN steps\")\n", " return (x2, abs(dx), abs(fx2), N, true)\n", "end" ] }, { "cell_type": "code", "execution_count": 2, "id": "d2962bac", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "running the secant method...\n", "step : root |dx| |f(x)|\n", " 1 : 1.5521908171562901e+00 5.42e-01 1.86e-02\n", " 2 : 1.5716418753086792e+00 1.95e-02 8.46e-04\n", " 3 : 1.5707962802269275e+00 8.46e-04 4.66e-08\n", " 4 : 1.5707963267949021e+00 4.66e-08 5.49e-15\n", "succeeded after 4 steps" ] }, { "data": { "text/plain": [ "(1.570796326794902, 4.656797464636609e-8, 5.489882783168415e-15, 4, false)" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" }, { "name": "stdout", "output_type": "stream", "text": [ "\n" ] } ], "source": [ "result = secant(cos,pi/4,2*pi/3)" ] }, { "cell_type": "markdown", "id": "28ed979b", "metadata": {}, "source": [ "Observe that the method converges fast, as fast as the quadratic convergence of Newton's method, although the convergence is in theory superlinear." ] }, { "cell_type": "markdown", "id": "40989276", "metadata": {}, "source": [ "# 2. A Bad Case for the Secant Method" ] }, { "cell_type": "code", "execution_count": 3, "id": "8405a14f", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "\n", "(process:71484): GLib-GIO-WARNING **: 19:10:04.270: Unexpectedly, UWP app `Evernote.Evernote_10.44.8.0_x64__q4d96b2w5wcc2' (AUMId `Evernote.Evernote_q4d96b2w5wcc2!Evernote') supports 1 extensions but has no verbs\n" ] } ], "source": [ "using Plots" ] }, { "cell_type": "code", "execution_count": 4, "id": "cb8e2353", "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "f(x) = -(x-2)*(x+2)\n", "plot(f,-3:0.01:3, label=\"f\")\n", "plot!([1], [f(1)], line=(:scatter, :black), label=\"\")\n", "plot!([1], [f(1)], line=(:sticks, :black), label=\"\")\n", "plot!([-1], [f(-1)], line=(:scatter, :black), label=\"\")\n", "plot!([-1], [f(-1)], line=(:sticks, :black), label=\"\")\n", "plot!([-1], [0], line=(:scatter, :black), label=\"\")\n", "annotate!(-1,-0.5,text(\"x3\"))\n", "plot!([1], [0], line=(:scatter, :black), label=\"\")\n", "annotate!(1,-0.5,text(\"x2\"))\n", "plot!([-1, 1], [0,f(1)], line=(:green), label=\"secant1\")\n", "plot!([-5/2,-1],[-9/4,0], line=(:red), label=\"\")\n", "plot!([-5/2], [-9/4], line=(:scatter, :black), label=\"\")\n", "plot!([-5/2], [-9/4], line=(:sticks, :black), label=\"\")\n", "plot!([-5/2], [0], line=(:scatter, :black), label=\"\")\n", "annotate!(-5/2,0.5,text(\"x1\"))\n", "plot!([-3,3],[3,3],line=(:red), label=\"secant2\")" ] }, { "cell_type": "markdown", "id": "12c16864", "metadata": {}, "source": [ "The plot above shows that the horizontal red line has no intersection with the horizontal axis and thus the next point in the secant method is undefined." ] }, { "cell_type": "markdown", "id": "7dedb31a", "metadata": {}, "source": [ "# 3. The Method of False Position" ] }, { "cell_type": "markdown", "id": "e2419b14", "metadata": {}, "source": [ "The secant method fails if the function value at two consecutive approximations is the same." ] }, { "cell_type": "code", "execution_count": 5, "id": "0047c644", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "regulafalsi" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using Printf # to format floating-point numbers\n", "\n", "\"\"\"\n", "Applies regula falsi to the function f on [a,b],\n", "\n", "ON ENTRY :\n", "\n", " f continuous function on [a,b]\n", " a, b [a, b] brackets a root, f(a)*f(b) < 0\n", " dxtol tolerance on the forward error,\n", " stops when |b - a| < dxtol\n", " fxtol tolerance on the backward error,\n", " stops when |min(f(a),f(b))| < fxtol\n", " N maximum number of steps\n", "\n", "ON RETURN : (A, B, steps, fail)\n", "\n", " root approximation of the root\n", " A, B [A, B] brackets a root\n", " steps number of steps\n", " fail true if neither of the stop criteria\n", " involving dxtol and fxtol are met.\n", "\n", "Example : \n", "\n", " f(x) = -(x-2)*(x+2)\n", " (a,b,n,fail) = regularfalsi(f,-3.0,1.0)\n", "\"\"\"\n", "function regulafalsi(f::Function,a::Float64,b::Float64,\n", " dxtol::Float64=1.0e-4,\n", " fxtol::Float64=1.0e-4,N::Int=10)\n", " println(\"running the method of false position...\")\n", " println(\"--------------------------------------------\")\n", " println(\" c |f(c)| |b-a| \")\n", " println(\"--------------------------------------------\")\n", " c = 0.0\n", " for i = 1:N \n", " c = (b*f(a) - a*f(b))/(f(a) - f(b))\n", " fc = f(c)\n", " if f(a)*fc < 0\n", " b = c\n", " else\n", " a = c\n", " end\n", " strc = @sprintf(\"%.16e\", c)\n", " strfc = @sprintf(\"%.2e\", fc)\n", " strba = @sprintf(\"%.2e\", b-a)\n", " println(\"$strc $strfc $strba\")\n", " if (abs(fc) < fxtol) | ((b-a) < dxtol)\n", " fail = false\n", " stri = string(i)\n", " println(\"succeeded after $stri steps\")\n", " return (c, a, b, i, fail)\n", " end\n", " end\n", " strN = string(N)\n", " println(\"failed requirements after $strN steps\")\n", " fail = true\n", " return (c, a, b, N, fail)\n", "end" ] }, { "cell_type": "code", "execution_count": 6, "id": "c5679bff", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "running the method of false position...\n", "--------------------------------------------\n", " c |f(c)| |b-a| \n", "--------------------------------------------\n", "-5.0000000000000000e-01 3.75e+00 2.50e+00\n", "-1.5714285714285714e+00 1.53e+00 1.43e+00\n", "-1.9062499999999998e+00 3.66e-01 1.09e+00\n", "-1.9808917197452227e+00 7.61e-02 1.02e+00\n", "-1.9961636828644502e+00 1.53e-02 1.00e+00\n", "-1.9992321474276937e+00 3.07e-03 1.00e+00\n", "-1.9998464058980137e+00 6.14e-04 1.00e+00\n", "-1.9999692802359281e+00 1.23e-04 1.00e+00\n", "-1.9999938560094370e+00 2.46e-05 1.00e+00\n", "succeeded after 9 steps\n" ] }, { "data": { "text/plain": [ "(-1.999993856009437, -3.0, -1.999993856009437, 9, false)" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "f(x) = -(x-2)*(x+2)\n", "(a,b,n,fail) = regulafalsi(f,-3.0,1.0)" ] }, { "cell_type": "markdown", "id": "9b782a66", "metadata": {}, "source": [ "The above experiment shows that the method of false position succeeds." ] }, { "cell_type": "code", "execution_count": 7, "id": "6b574bbf", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "running the method of false position...\n", "--------------------------------------------\n", " c |f(c)| |b-a| \n", "--------------------------------------------\n", "8.0000000000000004e-01 4.32e-01 1.80e+00\n", "6.4233576642335766e-01 4.03e-01 1.64e+00\n", "5.0724081625369188e-01 3.77e-01 1.51e+00\n", "3.9079015186250143e-01 3.40e-01 1.39e+00\n", "2.9297471287509347e-01 2.93e-01 1.29e+00\n", "2.1394906994743529e-01 2.39e-01 1.21e+00\n", "1.5268548979574564e-01 1.86e-01 1.15e+00\n", "1.0694125264175376e-01 1.39e-01 1.11e+00\n", "7.3828661942057047e-02 1.00e-01 1.07e+00\n", "5.0428827117398115e-02 7.07e-02 1.05e+00\n", "failed requirements after 10 steps\n" ] }, { "data": { "text/plain": [ "(0.050428827117398115, -1.0, 0.050428827117398115, 10, true)" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "f(x) = x^3 - 2*x^2 + 1.5*x\n", "(a,b,n,fail) = regulafalsi(f,-1.0,1.0)" ] }, { "cell_type": "markdown", "id": "7f371eb3", "metadata": {}, "source": [ "The above experiment shows that the convergence of the method of false position can be really slow." ] } ], "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": 5 }