{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Bisection method"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"%config InlineBackend.figure_format = 'svg'\n",
"from numpy import exp,sin,linspace,sign,abs\n",
"from matplotlib.pyplot import figure,plot,grid,xlabel,ylabel"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let us define three different functions for which we will find the root."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"def f1(x):\n",
" f = exp(x) - sin(x)\n",
" return f"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"def f2(x):\n",
" f = x**2 - 4.0*x*sin(x) + (2.0*sin(x))**2\n",
" return f"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"def f3(x):\n",
" f = x**2 - 4.0*x*sin(x) + (2.0*sin(x))**2 - 0.5\n",
" return f"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The following function implements the bisection method. Note that it takes some default arguments."
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"def bisect(fun,a,b,N=100,eps=1.0e-4,delta=1.0e-4,debug=False):\n",
" fa, fb = fun(a), fun(b)\n",
" sa, sb = sign(fa), sign(fb)\n",
" \n",
" if abs(fa) < delta:\n",
" return (a,0)\n",
"\n",
" if abs(fb) < delta:\n",
" return (b,0)\n",
"\n",
" # check if interval is admissible\n",
" if fa*fb > 0.0:\n",
" if debug:\n",
" print(\"Interval is not admissible\\n\")\n",
" return (0,1)\n",
"\n",
" for i in range(N):\n",
" e = b-a\n",
" c = a + 0.5*e\n",
" if abs(e) < eps*abs(c):\n",
" if debug:\n",
" print(\"Interval size is below tolerance\\n\")\n",
" return (c,0)\n",
" fc = fun(c)\n",
" if abs(fc) < delta:\n",
" if debug:\n",
" print(\"Function value is below tolerance\\n\")\n",
" return (c,0)\n",
" sc = sign(fc)\n",
" if sa != sc:\n",
" b = c\n",
" fb= fc\n",
" sb= sc\n",
" else:\n",
" a = c\n",
" fa= fc\n",
" sa= sc\n",
" if debug:\n",
" print(\"%5d %16.8e %16.8e %16.8e\"%(i+1,c,abs(b-a),abs(fc)))\n",
" \n",
" # If we reached here, then there is no convergence\n",
" print(\"No convergence in %d iterations !!!\" % N)\n",
" return (0,2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## First function"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 1 -3.00000000e+00 1.00000000e+00 1.90907076e-01\n",
" 2 -3.50000000e+00 5.00000000e-01 3.20585844e-01\n",
" 3 -3.25000000e+00 2.50000000e-01 6.94209267e-02\n",
" 4 -3.12500000e+00 1.25000000e-01 6.05288259e-02\n",
" 5 -3.18750000e+00 6.25000000e-02 4.61629389e-03\n",
" 6 -3.15625000e+00 3.12500000e-02 2.79283147e-02\n",
" 7 -3.17187500e+00 1.56250000e-02 1.16471966e-02\n",
" 8 -3.17968750e+00 7.81250000e-03 3.51301957e-03\n",
" 9 -3.18359375e+00 3.90625000e-03 5.52273640e-04\n",
" 10 -3.18164062e+00 1.95312500e-03 1.48021741e-03\n",
" 11 -3.18261719e+00 9.76562500e-04 4.63932552e-04\n",
"Function value is below tolerance\n",
"\n"
]
}
],
"source": [
"M = 100 # Maximum number of iterations\n",
"eps = 1.0e-4 # Tolerance on the interval\n",
"delta = 1.0e-4 # Tolerance on the function\n",
"a, b = -4.0, -2.0\n",
"r,status = bisect(f1,a,b,M,eps,delta,True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Second function"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Interval is not admissible\n",
"\n"
]
}
],
"source": [
"M = 100 # Maximum number of iterations\n",
"eps = 1.0e-4 # Tolerance on the interval\n",
"delta = 1.0e-4 # Tolerance on the function\n",
"a, b = -4.0, -2.0\n",
"r,status = bisect(f2,a,b,M,eps,delta,True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Lets visualize this function."
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"x = linspace(-3,3,500)\n",
"y = f2(x)\n",
"plot(x,y)\n",
"grid(True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"It looks like there are double roots which bisection method cannot compute such roots since the function value does not change around the root !!!"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Third function"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 1 -5.00000000e-01 2.50000000e+00 2.89455689e-01\n",
" 2 -1.75000000e+00 1.25000000e+00 4.52488254e-01\n",
" 3 -2.37500000e+00 6.25000000e-01 4.75412891e-01\n",
" 4 -2.06250000e+00 3.12500000e-01 4.10335430e-01\n",
" 5 -2.21875000e+00 1.56250000e-01 1.10488056e-01\n",
" 6 -2.29687500e+00 7.81250000e-02 1.42093933e-01\n",
" 7 -2.25781250e+00 3.90625000e-02 6.27310346e-03\n",
" 8 -2.23828125e+00 1.95312500e-02 5.44180845e-02\n",
" 9 -2.24804688e+00 9.76562500e-03 2.46591833e-02\n",
" 10 -2.25292969e+00 4.88281250e-03 9.34083833e-03\n",
" 11 -2.25537109e+00 2.44140625e-03 1.57095736e-03\n",
" 12 -2.25659180e+00 1.22070312e-03 2.34178305e-03\n",
" 13 -2.25598145e+00 6.10351562e-04 3.83092531e-04\n",
" 14 -2.25567627e+00 3.05175781e-04 5.94512219e-04\n",
" 15 -2.25582886e+00 1.52587891e-04 1.05854829e-04\n",
"Interval size is below tolerance\n",
"\n"
]
}
],
"source": [
"M = 100 # Maximum number of iterations\n",
"eps = 1.0e-4 # Tolerance on the interval\n",
"delta = 1.0e-4 # Tolerance on the function\n",
"a, b = -3.0, +2.0\n",
"r,status = bisect(f3,a,b,M,eps,delta,True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We dont need to specify all arguments to the function as some of them have default values."
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 1 -5.00000000e-01 2.50000000e+00 2.89455689e-01\n",
" 2 -1.75000000e+00 1.25000000e+00 4.52488254e-01\n",
" 3 -2.37500000e+00 6.25000000e-01 4.75412891e-01\n",
" 4 -2.06250000e+00 3.12500000e-01 4.10335430e-01\n",
" 5 -2.21875000e+00 1.56250000e-01 1.10488056e-01\n",
" 6 -2.29687500e+00 7.81250000e-02 1.42093933e-01\n",
" 7 -2.25781250e+00 3.90625000e-02 6.27310346e-03\n",
" 8 -2.23828125e+00 1.95312500e-02 5.44180845e-02\n",
" 9 -2.24804688e+00 9.76562500e-03 2.46591833e-02\n",
" 10 -2.25292969e+00 4.88281250e-03 9.34083833e-03\n",
" 11 -2.25537109e+00 2.44140625e-03 1.57095736e-03\n",
" 12 -2.25659180e+00 1.22070312e-03 2.34178305e-03\n",
" 13 -2.25598145e+00 6.10351562e-04 3.83092531e-04\n",
" 14 -2.25567627e+00 3.05175781e-04 5.94512219e-04\n",
" 15 -2.25582886e+00 1.52587891e-04 1.05854829e-04\n",
"Interval size is below tolerance\n",
"\n"
]
}
],
"source": [
"r,status = bisect(f3,a,b,debug=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Exercise\n",
"\n",
"We have used a `for` loop to perform the iterations. Modify the code to use `while` loop instead."
]
}
],
"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.7"
}
},
"nbformat": 4,
"nbformat_minor": 4
}