{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"# Preliminaries \n",
"\n",
" \n",
"**Anne Kværnø**\n",
"\n",
"Date: **Oct 19, 2018**\n",
"\n",
"# Introduction\n",
"In your mathematical courses so far, you have learned how to solve different\n",
"mathematical problems, like linear and nonlinear equations and differential\n",
"equations. You have learned how to differentiate and to integrate functions. \n",
"Unfortunately, only simplified models from \"real life\" applications can be\n",
"treated by these techniques, for more complex and realistic problems the best we\n",
"can aim for is some kind of an approximate solution, found on a computer by some\n",
"clever numerical algorithms. In this part of the course, our concern is how to\n",
"develope, analyse, implement and test a selection of such algorithms. \n",
"\n",
"In this note, we will present some mathematical results that will be \n",
"used frequently, as well as some definitions and concepts. Most of it should be known from\n",
"previous courses, in particular from Mathematics 1 and 3. \n",
"\n",
"Whenever some theoretical error analysis has been established, it should be\n",
"verified numerically. This is done based on a test problem with a known exact\n",
"solution, so the error in our numerical experiment can be evaluated, and the\n",
"theory numerically verified. Examples of such verifications will be demonstrated. \n",
"\n",
"# Vector spaces and norms\n",
"**Real vector space.**\n",
"\n",
"A *real vector space* is a set $V$ together with the operations $+$ (addition) and $\\cdot$ (multiplication with a scalar) which for all $x,y,z \\in V$ and $\\alpha, \\beta \\in \\mathbb{R}$ satisfy the following conditions:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{array}{rlrl}\n",
"1) & x+y \\in V \\\\ \n",
"2) & x+y=y+x &\n",
"3) & x+(y+z)=(x+y)+z \\\\ \n",
"4) & \\text{There exist a $0\\in V$ such that } x+0=x \\\\ \n",
"5) & \\text{For all $x\\in V$ there is a $-x \\in V$ such that $x+(-x)=0$} \\\\ \n",
"6) & \\alpha x \\in V &\n",
"7) & \\alpha (\\beta x) = (\\alpha \\beta) x \\\\ \n",
"8) & 1 x = x &\n",
"9) & \\alpha (x+y) = \\alpha x + \\alpha y \\\\ \n",
"10) & (\\alpha + \\beta)x = \\alpha x + \\beta x\n",
"\\end{array}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Examples:** \n",
"The following vector spaces will be used thoughout the course:\n",
"\n",
"* $\\mathbb{R}^m$ is the set of all real vectors with $m$ components. \n",
"\n",
"* $\\mathbb{R}^{m\\times n}$ is the set of all $m\\times n$ real matrices.\n",
"\n",
"* $C^m[a,b]$ is the set of all functions with continuous first $m$ derivatives on the interval $[a,b]$. It is common to use $C[a,b]$ rather than $C^0[a,b]$ for all continuous functions. \n",
"\n",
"* $\\mathbb{P}_n$ is the set of all polynomials of degree $n$ or less. \n",
"\n",
"Notice that $C^n[a,b] \\subset C^m[a,b]$ if $n > m$. \n",
"Further, $\\mathbb{P}_n \\subset C^{\\infty}[\\mathbb{R}]$.\n",
"\n",
"### Norms\n",
"\n",
"The norm $\\|\\cdot\\|$ of an element $x$ in a vector space $V$ is essentially a measure of the size of the element. \n",
"The norm obeys the following rules:\n",
"\n",
"**Norm $\\|\\cdot\\|$.**\n",
"\n",
"For all $x,y \\in V$ and $\\alpha \\in \\mathbb{R}$ the following holds"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"\\|x\\| & \\geq 0 && \\text{and} \\qquad \n",
" \\|x\\|=0 \\quad \\Leftrightarrow \\quad x \\equiv 0 \\\\ \n",
"\\|\\alpha x\\| &= |\\alpha |\\, \\|x\\|, \\\\ \n",
"\\|x+y\\| &\\leq \\|x\\|+\\|y\\| && (\\text{Triangle inequality}). \n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Norms for $C[a,b]$ and $\\mathbb{R}^m$:**\n",
"* For $f\\in C[a,b]$: \n",
"\n",
" * either $\\|f\\|_{\\infty} = \\max_{x\\in [a,b]}|f(x)|$.\n",
"\n",
" * or $\\|f\\|_2 = \\sqrt{\\int_a^b f(x)^2dx}$. \n",
"\n",
"\n",
"* For $\\mathbf{x} \\in \\mathbb{R}^m$, we will use: \n",
"\n",
" * either $\\|\\mathbf{x}\\|_{\\infty} = \\max_{i=1}^m |x_i|$, \n",
"\n",
" * or $\\|\\mathbf{x}\\|_2 = \\sqrt{\\sum_{i=1}^m x_i^2}$.\n",
"\n",
"\n",
"Note that we will use bold symbols for vectors $\\mathbf{x} \\in \\mathbb{R}^m$. \n",
"\n",
"As is demonstrated in these examples, the norm to which we\n",
"refer is usually marked with a subscript. \n",
"We will always use the absolute value as the norm of a real number, thus\n",
"$\\|x\\|=|x|$ whenever $x\\in \\mathbb{R}$. \n",
"\n",
"**Example 1:**\n",
"Let $\\mathbf{x}=[1,-6,3,-1,5]^T\\in \\mathbb{R}^5$. Then"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
" \\|\\mathbf{x}\\|_2 &= \\sqrt{1+36+9+1+25} = 8.4853, \\\\ \n",
" \\|\\mathbf{x}\\|_{\\infty} &= \\max\\{1,6,3,1,5\\}= 6.\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Norms of vectors in $\\mathbb{R}^m$ are implemented in Python."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"\n",
"from numpy import *\n",
"from numpy.linalg import norm\n",
"from matplotlib.pyplot import *\n",
"newparams = {'figure.figsize': (8.0, 4.0), 'axes.grid': True,\n",
" 'lines.markersize': 8, 'lines.linewidth': 2,\n",
" 'font.size': 14}\n",
"rcParams.update(newparams)"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# The norm of a vector in R^n\n",
"x = np.array([1, -6, 3, -1, 5])\n",
"nx_2 = norm(x) # The 2-norm is the default\n",
"nx_inf = norm(x,ord=inf) # The max-norm\n",
"print('x = \\n', x)\n",
"print('The 2-norm of x: {:8.4f}'.format(nx_2))\n",
"print('The max-norm of x: {:8.4f}'.format(nx_inf))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Example 2:**\n",
"Let $f(x)=\\sin(x)$ on $[0,2\\pi]$, so $f\\in C^{\\infty}[0,2\\pi]$. In this case"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
" \\|f\\|_2 &= \\displaystyle \\sqrt{\\int_0^{2\\pi} \\sin^2(x) dx} = \\sqrt{\\pi} = 1.7725 \\\\ \n",
" \\|f\\|_{\\infty} &= \\displaystyle \\max_{x\\in[0,2\\pi]}|\\sin(x)| = 1.\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Norms on function spaces require some sort of numerical approximations. In this case, the interval $[0,2\\pi]$ has been divided into $N=1000$ uniform subintervals, so $x_i=(2\\pi i)/N$, $i=0,1,\\dots,N$. The norm $\\|f\\|_{\\infty}$ is approximated by $\\max_i|f(x_i)|$, and the integral required for $\\|f\\|_2$ is computed by the trapezoidal rule."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Approximations to the norm of a function in C[a,b]\n",
"def f(x): # The function\n",
" return sin(x)\n",
"a, b = 0, 2*pi # The interval\n",
"\n",
"x = linspace(a, b, 1001) # The values of x on [a,b]\n",
"norm_of_f_2 = sqrt(trapz(f(x)**2, x=x)) # Approximate the integral by the trapezoidal rule\n",
"norm_of_f_inf = max(abs(f(x))) # Approximate the max-norm by max|f(x_i)| \n",
"\n",
"print('The 2-norm of f: {:8.4f}'.format(norm_of_f_2))\n",
"print('The max-norm of f: {:8.4f}'.format(norm_of_f_inf))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Convergence and errors\n",
"\n",
"**Convergence of a sequence.**\n",
"\n",
"Let $\\{x_k\\}_{k=0}^{\\infty}$ be an infinite sequence of real numbers. The sequence converges to $x$, if, for any $\\varepsilon >0$ there exist a positive integer $N(\\varepsilon)$ such that $|x_k - x|<\\varepsilon$ whenever $k > N(\\varepsilon)$. \n",
"\n",
"Common notations:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\lim_{k\\rightarrow\\infty} x_k = x \\qquad \\text{or} \\qquad x_k \\rightarrow x \\text{ as } k\\rightarrow \\infty.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Example 3:**\n",
"It is known that the sequence $x_k=(1+1/k)^k \\rightarrow e = 2.7182\\dotsc$ as $k\\rightarrow\n",
"\\infty$. The sequence is monotone, so $|x_{k+1}-e|<|x_k-e|$. The following\n",
"program demonstrates the concept of convergence for this sequence: Given an $\\varepsilon$, the\n",
"positive integer $N(\\epsilon)$ is returned."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Demonstrate the convergence of a sequence \n",
"epsilon = 1.e-3\n",
"x_exact = exp(1) # Exact solution\n",
"Nmax = 10**7 # Maximum N, to avoid infinite loops.\n",
"step = 1 # increments of k. Make it large for epsilon small\n",
"for k in range(1, Nmax, step): # \n",
" xk = (1+1/k)**k\n",
" if abs(xk-x_exact) < epsilon:\n",
" break\n",
"if k+1==Nmax:\n",
" print('Maximum N reached')\n",
"print('epsilon = {:.2e}, N = {:d}'.format(epsilon, k))\n",
"print('xk = {:f}, |xk-e| = {:e}'.format(xk, abs(xk-x_exact)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Clearly, this code can only be used for monotone series. \n",
"\n",
"Warning: The sequence converges very slowly, so given a very small `epsilon`\n",
"will result in $N$ too large. This can be compensated by making `Nmax` larger,\n",
"and increment $k$ in larger steps.\n",
"\n",
"## Convergence of an iterative process\n",
"Let $X$ be the exact solution of a problem, and $X_k$ a numerical approximation\n",
"achieved by some iterative process $X_{k+1} = G(X_k)$. In this case the\n",
"iterations converge towards $X$ if"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\lim_{k\\rightarrow \\infty} \\|X-X_k\\| = 0\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let $e_k = \\|X-X_k\\|$ measure the error. In practice, you have to choose an\n",
"appropriate norm, which depends on the problem and what you might be interesting\n",
"in measuring. The order of convergence is $p$ if there exist a positive constant $M$ such that"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"e_{k+1} \\leq M e_k^p\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Notation:**\n",
"The case $p=1$ is called *linear* convergence, $p=2$ is called *quadratic* convergence and $p=3$ *cubic* convergence. \n",
"\n",
"### Numerical verification of the order\n",
"\n",
"To verify the order, we make the assumptions that $e_{k+1} = C_k e_k^p$, and\n",
"that the $C_k$ do not change much from one iteration to the next one. \n",
"This assumption is usually reasonable when the error becomes small.\n",
"The order $p$ can then be computed numerically by the following procedure:\n",
"Take the expressions for the error for two subsequent iterations, assuming that\n",
"$C_{k+1} \\approx C_k \\approx C$. Then divide the expression \n",
"by each other to get rid of the unknown constant $C$, take the logarithm on both sides and\n",
"solve for the order $p$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{array}{rcl}\n",
" e_{k+2} &\\approx& C e_{k+1}^p \\\\ \n",
" e_{k+1} &\\approx& C e_{k}^p \n",
" \\end{array}\n",
" \\quad \\Rightarrow \\quad \n",
" \\frac{e_{k+2}}{e_{k+1}} \\approx \\left( \\frac{e_{k+1}}{e_{k}} \\right)^p \n",
" \\quad \\Rightarrow \\quad \n",
" \\log \\left( \n",
" \\frac{e_{k+2}}{e_{k+1}}\\right) \\approx p \\log \\left( \\frac{e_{k+1}}{e_{k}} \\right)\n",
" \\quad \\Rightarrow \\quad \n",
" p \\approx \\frac{\\log{(e_{k+2}/e_{k+1})}}{\\log{(e_{k+1}/e_{k})}}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We are usually not so interested in the constant $C$, but given $p$ and the\n",
"error for two iterations, this can easily be approximated. \n",
"\n",
"**Example 4:**\n",
"Newton's method applied to the equation $f(x)=0$ is given by"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"x_{k+1} = x_k - \\frac{f(x_k)}{f'(x_k)}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let $r$ be a solution of the equation. \n",
"It can be proved that the error $r-x_k$ satisfies"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"r-x_{k+1} = -\\frac{f''(\\xi_k)}{2 f'(x_k)}(r-x_k)^2, \\qquad \\text{where $\\xi_k$ a\n",
"real number between $x_k$ and $r$.}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Since the error in this case is a real number, its norm is $e_k = |r-x_k|$\n",
"and"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"e_{k+1} = C_k e_{k}^2 \\quad \\text{where} \\quad C_k = \\frac{|f''(\\xi_k)|}{2 |f'(x_k)|}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Notice that $C_k \\rightarrow |f''(r)|/(2|f'(r)|$ as $x_k \\rightarrow r$. \n",
"\n",
"The constant $M$ is given by the upper\n",
"bound of $C_k$. Or more precisely: let $I_{\\delta}=[x-\\delta,x+\\delta]$ be some interval \n",
"around the solution $r$. Assume there exists \n",
"constants $L$ and $K$ such that \n",
"$ |f'(x)|\\leq L$ and $|f''(x)|\\geq K$ for all $x\\in I_{\\delta}$. Then $M=K/(2L)$\n",
"and"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"e_{k+1} \\leq M e_k^2.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The convergence is quadratic, and the iterations converge for all starting\n",
"values $x_0$ chosen such that"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"M e_0 < 1,\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"often described as \"sufficiently close to the solution\". \n",
"\n",
"Let us now verify the theoretical result by applying Newtons method \n",
" to the problem $x^2-a=0$ for some $a>0$. The iterations become"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"x_{k+1} = x_k - \\frac{x_k^2-a}{2x_k} = \\frac{x_{k}^2+a}{2x_k}, \\qquad k=0,1,2,\\dotsc.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"with the exact solution $r=\\sqrt{a}$. \n",
"From the discussion above, we expect $C\\approx |f''(r)/(2f'(r)|$ with\n",
"$r=\\sqrt{a}$, which in our case becomes $C \\approx 1/(2\\sqrt{a})$. Use the\n",
"following code to see if the theoretical considerations hold in practice:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Order of convergence for iterations\n",
"# Test problem: Newtons method for x^2-a=0\n",
"a = 4\n",
"def g(x):\n",
" return (x**2+a)/(2*x) # g(x) = x-f(x)/f'(x)\n",
"x_exact = sqrt(a) # Exact solution\n",
"x = 1 # Starting value\n",
"errors = [abs(x_exact-x)] # Array to store errors\n",
"Nit = 10 # Number of iterations\n",
"\n",
"# Start the iterations\n",
"print('The Newton iterations:')\n",
"for k in range(Nit):\n",
" x = g(x) # One iteration\n",
" ek = abs(x_exact-x) # Find the error\n",
" print('k = {:2d}, x_k = {:10.8f}, e_k = {:8.2e}'.format(k, x, ek))\n",
" if ek < 1.e-15: # If the error is small, terminate.\n",
" Nit = k+1\n",
" break \n",
" errors.append(ek) # Append the new error to the array of errors\n",
"\n",
"# Find the order and the error constant C\n",
"print('\\nThe order p and the error constant C')\n",
"for k in range(Nit-2):\n",
" p = log(errors[k+2]/errors[k+1])/log(errors[k+1]/errors[k])\n",
" C = errors[k+2]/errors[k+1]**p\n",
" print('k = {:2d}, p = {:4.2f}, C = {:6.4f}'.format(k, p, C))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Convergence of $h$-dependent approximations\n",
"Let $X$ be the exact solution, and $X(h)$ some numerical solution depending on a\n",
"parameter $h$, and let $e(h)$ be the norm of the error, so $e(h)=\\|X-X(h)\\|$. The numerical approximation $X(h)$ converges to $X$ if $e(h) \\rightarrow 0$ as $h\\rightarrow 0$. \n",
"The order of the approximation is $p$ if there exists a positive constant $M$ such that"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"e(h) \\leq M h^p\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**The Big $\\mathcal{O}$-notation:**\n",
"A function $f(x) = \\mathcal{O}(g(x))$ as $x\\rightarrow a$ if and only there exist positive numbers $\\delta$ and $M$ such that"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"|f(x)| \\leq M|g(x)| \\qquad \\text{when} \\qquad 0 < |x-a|<\\delta.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let $a=0$ and $f(h)=e(h)$, thus the error of an approximation of order $p$ can be written as"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"E(h) = \\mathcal{O}(h^p).\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This is often used when we are not directly interested in any expression for the constant $M$, we only need to know it exists. \n",
"\n",
"### Numerical verification\n",
"\n",
"The following is based on the assumption that $e(h)\\approx C h^p$ for some\n",
"unknown constant $C$. This assumption is usually reasonable for sufficiently\n",
"small $h$. \n",
"\n",
"Choose a test problem for which the exact solution is known and compute the\n",
"error for a sequence of smaller $h$'s, for instance $h_k=H/2^k$,\n",
"$k=0,1,2,\\dots$. The procedure is then quite similar to what was done for\n",
"iterative processes."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{array}{rcl} e(h_{k+1}) &\\approx& C h_{k+1}^p \\\\ e(h_k) &\\approx& C h_k^p \\end{array}\n",
" \\qquad \\Rightarrow \\qquad \n",
" \\frac{e(h_{k+1})}{e(h_k)} \\approx \\left( \\frac{h_{k+1}}{h_k} \\right)^p \n",
" \\qquad \\Rightarrow \\qquad \n",
" p \\approx \\frac{\\log{(e(h_{k+1})/e(h_k))}}{\\log{(h_{k+1}/h_k)}}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Since"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"e(h) \\approx Ch^p \\qquad \\Rightarrow \\qquad \\log{e(h)} \\approx \\log{C} + p \\log{h}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"a plot of $e(h)$ as a function of $h$ using a logarithmic scale on both axes\n",
"will be a straight line with slope $p$. Such a plot is referred to as\n",
"an *error plot* or a *convergence plot*. \n",
"\n",
"### Some terminology\n",
"\n",
"Let $X$ be the exact solution of some problem, and $X(h)$ the numerical\n",
"approximation of that problem. The following concepts are of interest: \n",
"\n",
"* The error $E(h)$: $E(h) = X - X(h)$. This is something which obviously is only known if the exact solution is known (which it will be in our test problems, but not in real life problems). Still, most error analysis will start trying to find an expression for this error, but it will typically contain some higher order derivatives evaluated in some unknown point.\n",
"\n",
"* The error bound: Typically of the form $\\|X-X(h)\\|\\leq M h^p$. If $M$ is known, this can be used to decide how small $h$ has to be to guarantee that the error is below some tolerance. \n",
"\n",
"* Error estimate $\\mathcal{E} \\approx E$. This is an approximation to the error, and something that can be computed, and included in practical codes. How to compute those will be described for each problem we will discuss later in the course. \n",
"\n",
"**Example 5:**\n",
"Consider the trapzoidal rule for numerical integration. It is known that"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\int_a^b f(x)dx = T(h) + E(h)\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where $T(h)$ is the numerical approximation given by"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"T(h) = h\\left(\\frac12 f(x_0)+\\sum_{i=0}^n f(x_i)+\\frac12 f(x_n) \\right), \\qquad x_n=a+ih, \\quad h = \\frac{b-a}{n},\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"and the error $E(h)$ is known to be"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"E(h) = -\\frac{b-a}{12} f''(\\xi) h^2, \\qquad \\xi \\in (a,b).\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Assume there exist an $M$ such that $|f''(x)|\\leq M$ for all $x\\in (a,b)$. Let\n",
"$e(h)=|E(h)|$ (notice that $E(h)$ is a scalar) so"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"e(h) \\leq M h^2.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"So the error of the trapezoidal rule is of order 2, and $E(h) = \\mathcal{O}(h^2)$.\n",
"\n",
"Use this to verify the order of the trapezoidal rule, as given above. As test example, choose"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\int_0^{\\pi}\\sin(x)dx = 2.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this case, we expect the order to be $p$. The constant $C = |f''(\\xi)|\\pi/12$ for some unknown $\\xi\\in[0,\\pi]$. Thus $0 < C <\n",
"\\pi/12=0.2617\\dotsc$, but we can not be more precise. The upper bound for the error is $e(h)\\leq \\pi/12\\; h^2$. \n",
"\n",
"\n",
"The following code can be used to confirm the result. It also returns an\n",
"approximation to $C$, so we can at least check if it is within the expected\n",
"bound."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Find the order and the error plot for the trapezoidal rule\n",
"\n",
"def trapezoidal(f, a, b, n):\n",
" # The trapezoidal rule\n",
" h = (b-a)/n # The stepsize\n",
" x = linspace(a, b, n+1) # Create a uniform grid x = [a=x_0,x_1,...,x_{n-1},x_n=b]\n",
" res = h*sum(f(x[1:n])) # h*(f(x_1)+f(x_2)+...+f(x_{n-1})\n",
" res = res + 0.5*h*(f(a)+f(b)) # add the function values of the endpoints a and b\n",
" return res\n",
"\n",
"def f(x): # Define the function\n",
" return sin(x)\n",
"\n",
"a, b = 0, pi # integration interval\n",
"exact = 2\n",
"\n",
"# Find an numerical approximation for different values of h. \n",
"# Store the stepsize h and the error\n",
"n = 1 # initial stepsize, h=(b-a) \n",
"h = (b-a)/n\n",
"steps = [] # arrays to store stepsizes and errors\n",
"errors = []\n",
"Nmax = 10\n",
"for k in range(Nmax):\n",
" numres = trapezoidal(f, a, b, n) # Numerical approximation\n",
" eh = abs(exact - numres) # Error e(h)\n",
" print('h = {:8.2e}, T(h) = {:10.8f}, e(h) = {:8.2e}'.format(h, numres, eh))\n",
" steps.append(h) # Append the step to the array\n",
" errors.append(eh) # Append the error to the array\n",
" n = 2*n # Reduce the stepsize with a factor 2\n",
" h = (b-a)/n\n",
"\n",
"# Find the order and the error constant\n",
"print('\\nThe order p and the error constant C')\n",
"for k in range(1, Nmax-1):\n",
" p = log(errors[k+1]/errors[k])/log(steps[k+1]/steps[k])\n",
" C = errors[k+1]/steps[k+1]**p\n",
" print('h = {:8.2e}, p = {:4.2f}, C = {:6.4f}'.format(steps[k], p, C))\n",
"\n",
"# Make an error plot\n",
"clf()\n",
"loglog(steps, errors, 'o-')\n",
"xlabel('h')\n",
"ylabel('e(h)')\n",
"title('Error plot for the trapezoidal rule')\n",
"grid('True')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Taylor-expansions\n",
"Given a function $f\\in C^{\\infty}[a,b]$. Choose a point $x$ and an increment $h$ such that $x, x+h \\in [a,b]$. The Taylor expansion of $f$ around $x$ is then given by"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"f(x+h) = \\sum_{k=0}^{\\infty} \\frac{h^k}{k!}f^{(k)}(x).\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The function is called analytic in $x$ if the series converges for sufficiently small values of $|h|$.\n",
"In numerics, we will usually work with the trunctated series, also called the Taylor polynomial:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"f(x+h) = \\sum_{k=0}^{m} \\frac{h^k}{k!}f^{(k)}(x) + R_{m+1}(x)\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where the remainder term is given by"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"R_{m+1}(x) = \\frac{h^{m+1}}{(m+1)!} f^{(m+1)}(\\xi)\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where $\\xi$ is some unknown point between $x$ and $x+h$. The truncated series only require $f\\in C^{m+1}$ near $x$. The truncated expansion will often simply be written as"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"f(x+h) = \\sum_{k=0}^{m} \\frac{h^k}{k!}f^{(k)}(x) + \\mathcal{O}(h^{m+1}).\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Some other useful results\n",
"\n",
"**Result 1:**\n",
"Let $f\\in C[a,b]$ and let $u$ be some number between $f(a)$ and $f(b)$, then there exist at least one $\\xi \\in (a,b)$ such that $f(\\xi)=u$.\n",
"\n",
"**Result 2:**\n",
" Let $f\\in C[a,b]$. Given $k$ nodes $x_i\\in [a,b]$ and $k$ positive weights\n",
"$w_i>0$, $i=1,\\dotsc,k$. Then there exists at least one $\\xi \\in (a,b)$ such that"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"$$\n",
"\\begin{equation} \\sum_{i=1}^{k} w_i f(x_i) = f(\\eta) \\sum_{i=1}^k w_i. \n",
"\\label{eq:gen_mid} \\tag{1}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Result 3:**\n",
" Let $f\\in C^1[a,b]$. Then there exists at least one $\\xi\\in (a,b)$ such that"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"f'(\\xi) = \\frac{f(b)-f(a)}{b-a}\n",
"\\label{_auto1} \\tag{2}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Result 4:**\n",
"(Rolle's theorem) Let $f\\in C^1[a,b]$ and $f(a)=f(b)=0$. Then there exists at least one $\\xi \\in (a,b)$ such that $f'(\\xi)=0$."
]
}
],
"metadata": {},
"nbformat": 4,
"nbformat_minor": 2
}