{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"# Fixed-Point Iteration\n",
"\n",
"### Modules - Root Finding\n",
"\n",
"By Eilif Sommer Øyre, Jonas Themsland and Jon Andreas Støvneng \n",
"\n",
"\n",
"Last edited: March 11th 2018\n",
"\n",
"---"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"## Introduction\n",
"\n",
"An equation $f(x) = 0$ can always be written as $g(x) = x$. Fixed-point iteration can be applied to approximate the fixed number $r = g(r)$.\n",
"\n",
"After selecting an initial guess $x_0$, the fixed-point iteration is $x_{n+1} = g(x_n)$. That is,\n",
"\n",
"\\begin{equation*}\n",
"x_0 = \\textrm{initial guess} \\\\\n",
"x_1 = g(x_0) \\\\\n",
"x_2 = g(x_1) \\\\\n",
"x_3 = g(x_2) \\\\\n",
"\\vdots\n",
"\\end{equation*}\n",
"\n",
"The algorithm is repeated until either a testing condition $e_{i+1} = \\left|\\:x_{i+1} - x_{i}\\:\\right| < \\alpha$, where $\\alpha$ is some tolerance limit is met, or until a fixed number of iterations $N$ is reached. Note that the method may or may not converge. Note that the choice of $g(x)$ is in general not unique."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Example 1\n",
"\n",
"As a simple introductory example we use the fixed-point iteration to solve the equation\n",
"\n",
"$$\\frac{1}{2}\\sin(x) - x + 1 = 0.$$\n",
"\n",
"The most natural is to use $g(x)= \\sin(x)/2 + 1$ and then use the fixed-point iteration to solve $g(x)=x$."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"x1:\t1.84147\n",
"x2:\t1.96359\n",
"x3:\t1.92384\n",
"x4:\t1.93832\n",
"x5:\t1.93322\n",
"x6:\t1.93504\n",
"x7:\t1.93439\n",
"x8:\t1.93462\n",
"x9:\t1.93454\n"
]
}
],
"source": [
"import math\n",
"\n",
"x = 1 # initial guess\n",
"N = 10 # iterations\n",
"\n",
"for i in range(1, N):\n",
" x = math.sin(x) + 1\n",
" print(\"x%i:\\t%.5f\"%(i, x))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Convergence\n",
"\n",
"For the iteration scheme to return a fixed point $r$, it needs to converge. A criterion for convergence is that the error, $e_i = \\left|\\:r - x_i \\:\\right|$, decreases for each iteration step. This means that the change from $x_i$ to $x_{i+1}$ also has to decrease for each step. The convergence criterion is explained in the following theorem [1]:\n",
"\n",
"> **Theorem 1** Assume that $g$ is continuously differentiable, that $g(r) = r$, and that $S = \\left|\\:g'(r)\\:\\right| < 1$. Then fixed-point iteration converges linearly with rate $S$ to the fixed point $r$ for initial guesses sufficiently close to $r$. \n",
"\n",
"Note that in the example above we have $|g'(x)|=\\cos(x)/2 < 1$ which means that the method will converge for all initial guesses.\n",
"\n",
"The convergence can in fact be of higher order. This is explained in the following theorem [2]:\n",
"\n",
"> **Theorem 2** Assume that $g$ is $p$ times continuously differentiable and that the fixed-iteration converges to $r$ for some initial guess $x_0$. If $g'(r)=g''(r)=\\cdots=g^{(p-1)}(r)=0$ and $g^{(p-1)}(r)=0$, then the order of convergence is $p$.\n",
"\n",
"The convergence is said to be of order $p$ when $\\lim_{i\\to\\infty}e_{i+1}/e_{i}^p = \\text{constant}$. \n",
"\n",
"**Excercise:** Prove theorem 1 using the mean value theorem and the theorem 2 using Taylor's theorem."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Example 2: Babylonian Method for Finding Square Roots\n",
"\n",
"As mentioned in the introduction, the choice of $g(x)$ is far from unique. We will now use fixed-point iteration to estimate the square root of a real and positive number $a$. That is, we want to solve $f(x) = x^2-a^2=0$. Two natural choices for $g(x)$ are $g_1(x)=a/x$ and $g_2(x)=x$. However, none of these converges since $|g_1'(\\sqrt a)| = |g_2'(\\sqrt a)| =1$. We can choose the mean as $g(x)$:\n",
"\n",
"$$g(x) = \\frac{1}{2}\\left(\\frac{a}{x} + x\\right).$$\n",
"\n",
"In this case we have \n",
"\n",
"$$\\left|\\:g'(x)\\:\\right|_{x=\\sqrt a}=\\left|\\:\\frac{1}{2}(1-\\frac{a}{x^2})\\:\\right|_{x=\\sqrt a} = 0 < 1$$\n",
"\n",
"The theorem above thus implies that method converges for some initial guess $x_0$. The method is in fact globally convergent due to the convexity of $f(x)$."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"x1:\t0.53500\n",
"x2:\t0.33292\n",
"x3:\t0.27159\n",
"x4:\t0.26467\n",
"x5:\t0.26458\n",
"x6:\t0.26458\n",
"x7:\t0.26458\n",
"x8:\t0.26458\n",
"x9:\t0.26458\n"
]
}
],
"source": [
"a = 0.07 # square of root\n",
"x = 1 # initial guess\n",
"N = 10 # iterations\n",
"\n",
"for i in range(1, N):\n",
" x = (a/x + x)/2\n",
" print(\"x%i:\\t%.5f\"%(i, x))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This method can also be derived from Newton's method (see the section on Newton's method below). The method was however first used by the Babylonian people long before Newton [2]."
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"## Example 3: Adding Stopping Conditions\n",
"\n",
"Underneath follows an implementation of the Fixed-Point Iteration with a testing condition $\\left|\\:x_{i+1} - x_i\\:\\right| < \\alpha$. There are other stopping criteria that may be relevant such as the backward error $ \\left|\\:f(x_a)- f(0)\\right|$. We need an additional stopping criteria in case convergence fails. We therefore stop the computation if the a given number of iterations $N$ is reached."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def FPI(g, x, maxit=100, alpha=1e-5):\n",
" \"\"\" A function using fixed-point iteration to compute \n",
" an approximate solution to the fixed point problem g(r) = r.\n",
" \n",
" Arguments:\n",
" g callable function, g(x)\n",
" x float. The initial guess\n",
" alpha float. Error tolerance\n",
" maxit int. Maximum number of iterations\n",
" \n",
" Returns:\n",
" array of each fixed-point approximation\n",
" \"\"\"\n",
" \n",
" result = [x]\n",
" x_next = g(x)\n",
" i = 0 # Counter\n",
" e = abs(x_next - x)\n",
" \n",
" while e > alpha :\n",
" if i > maxit:\n",
" print(\"No convergence.\")\n",
" break\n",
" x = x_next\n",
" x_next = g(x_next)\n",
" result.append(x)\n",
" if x - x_next > e:\n",
" print(\"Divergence.\")\n",
" break\n",
" e = abs(x_next - x)\n",
" i += 1\n",
" return result"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We want to solve the equation,\n",
"\n",
"\\begin{equation*}\n",
"f(x) = x^4 + x - 1 = 0\n",
"%\\label{eq:example}\n",
"\\end{equation*}\n",
"\n",
"using fixed point iteration. First, we rewrite the equation to the form, $x = g(x)$. This particular equation may be rewritten in three ways:\n",
"\n",
"$$g_1(x) = x = 1 - x^4,$$\n",
"\n",
"$$g_2(x) = x = \\sqrt[4]{1 - x},$$\n",
"\n",
"$$g_3(x) = x = \\frac{1 + x^4}{1 + 2x^3}.$$\n",
"\n",
"Each of these has its own convergence rate. One of the solutions is clearly somewhere between 0 and 1, so we will use $x_0=0.5$ as initial guess."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"x_0 = 0.5\n",
"N = 100\n",
"alpha = 1e-5\n",
"\n",
"def g1(x): return 1 - x**4\n",
"def g2(x): return (1 - x)**(1/4)\n",
"def g3(x): return (1 + 3*x**4)/(1 + 4*x**3)\n",
"\n",
"def print_result(result):\n",
" for i in range(len(result)):\n",
" print(\"x%i:\\t%.5f\"%(i, result[i]))"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Divergence.\n",
"x0:\t0.50000\n",
"x1:\t0.93750\n"
]
}
],
"source": [
"result_g1 = FPI(g1, x_0, N, alpha)\n",
"print_result(result_g1)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"x0:\t0.50000\n",
"x1:\t0.84090\n",
"x2:\t0.63157\n",
"x3:\t0.77909\n",
"x4:\t0.68557\n",
"x5:\t0.74883\n",
"x6:\t0.70794\n",
"x7:\t0.73514\n",
"x8:\t0.71739\n",
"x9:\t0.72912\n",
"x10:\t0.72143\n",
"x11:\t0.72650\n",
"x12:\t0.72317\n",
"x13:\t0.72536\n",
"x14:\t0.72392\n",
"x15:\t0.72487\n",
"x16:\t0.72425\n",
"x17:\t0.72465\n",
"x18:\t0.72439\n",
"x19:\t0.72456\n",
"x20:\t0.72445\n",
"x21:\t0.72452\n",
"x22:\t0.72447\n",
"x23:\t0.72451\n",
"x24:\t0.72448\n",
"x25:\t0.72450\n"
]
}
],
"source": [
"result_g2 = FPI(g2, x_0, N, alpha)\n",
"print_result(result_g2)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"x0:\t0.50000\n",
"x1:\t0.79167\n",
"x2:\t0.72986\n",
"x3:\t0.72453\n",
"x4:\t0.72449\n"
]
}
],
"source": [
"result_g3 = FPI(g3, x_0, N, alpha)\n",
"print_result(result_g3)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The method does not converge for $g_1$, but both $g_2$ and $g_3$ converges towards the fixed number $0.72449$, which is a solution of the equation $f(x) = 0$. The method converges much faster for $g_3$ than for $g_2$.\n",
"\n",
"This result can be verified by calculating the convergence rates $S = \\left|\\: g'(r) \\:\\right|$:\n",
"\n",
"\\begin{equation*}\n",
"\\left|\\:g_1'(0.72449)\\:\\right| \\approx 1.521 > 1, \\\\\n",
"\\left|\\:g_2'(0.72449)\\:\\right| \\approx 0.657 < 1, \\\\\n",
"\\left|\\:g_3'(0.72449)\\:\\right| \\approx 0.000 < 1. \\\\\n",
"\\end{equation*}\n",
"\n",
"Theorem 1 tells us that the method will converge for $g_2$ and $g_3$ for some inital guess. Moreover, from theorem 2 we know that the method converges linearly for $g_2$ and the order of convergence for $g_3$ is two.\n",
"\n",
"**Exercise:** Verify that the convergence for $g_2$ is linear and that the convergence of $g_3$ is of second order by using the numerical results above."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## In comparison with the Newton-Rhapson method\n",
"\n",
"The Newton-Rhapson method is a special case of fixed-point iteration where the convergence rate is zero - the fastest possible. The method estimates the root of a differentiable function $f(x)$ by iteratively calculating the expression\n",
"\n",
"\\begin{equation*}\n",
"x_{n+1} = x_n -\\frac{f(x_n)}{f'(x_n)}.\n",
"\\end{equation*}\n",
"\n",
"By using theorem 2, we can show that Newton's method is of second order (if $f'(r)\\neq 0$). If $f''(r)=0$, then the method has a third order convergence.\n",
"\n",
"This Newton iteration may be rewritten as a fixed point iteration\n",
"\n",
"\\begin{equation*}\n",
"x_{n+1} = g(x_n), \\: n = 1, 2 ,3, ...\n",
"\\end{equation*}\n",
"\n",
"where\n",
"\n",
"\\begin{equation*}\n",
"r = g(r) = r - \\frac{f(r)}{f'(r)}\n",
"%\\label{eq:newton}\n",
"\\end{equation*}\n",
"\n",
"Provided that the iteration converges to a fixed point $r$ of $g$. From this we obtain\n",
"\n",
"$$ \\frac{f(r)}{f'(r)} = 0 \\: \\Rightarrow \\: f(r) = 0$$\n",
"\n",
"and thus the fixed point $r$ is a root of $f$.\n",
"\n",
"The Newton's method is further discussed in the module [Root Finding - Newton-Rhapson Method](https://nbviewer.jupyter.org/urls/www.numfys.net/media/notebooks/newton_raphson_method.ipynb)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Final comments\n",
"\n",
"We have now used the fixed-point iteration method to solve the equation $f(x) = 0$ by rewriting it as a fixed-point problem. This is a simple method and which is easy to implement, but unlike the [Bisection Method](https://nbviewer.jupyter.org/urls/www.numfys.net/media/notebooks/bisection_method.ipynb) it only converges if the initial guess is sufficiently close to the root $r$. It is in general only locally convergent. However, the convergence rate may, or may not, be faster than that of the Bisection Method, which is 1/2.\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"## References\n",
"\n",
"[1] Sauer, T.: Numerical Analysis international edition, second edition, Pearson 2014 \n",
"[2] Gautschi, W.: Numerical Analysis, second edition, Birkhäuser 2012 (https://doi.org/10.1007/978-0-8176-8259-0)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"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.5.4"
}
},
"nbformat": 4,
"nbformat_minor": 2
}