MA514 - Numerical Analysis - Spring 2017
  • ascher and greif, a first course in numerical methods

Floating Point Arithmetic

Denoted as \(fl(x)\) \(fl(x) = -1 (1 d_1 d_2 d_3 ...)_2 \cdot 2^{exp}\) sign bit - binary fraction - exponent

When a number is normalized, 1 is added to the binary fraction: \(fl(x) = (-1)^{d_s} (1 + d_1 d_2 d_3 ... d_n)_2 \cdot 2^{exp}\)

Let \(t = 2\) be the number of digits, \([L, U] = [-2, 1]\) be the range of possible exponents.

What are the possible values that can be stored in our floating point system?

Binary fractions

  • 1.00 = 1
  • 1.01 = 1 + 1/4 = 5/4
  • 1.10 = 1 + 1/2 = 3/2
  • 1.11 = 1 + 1/2 + 1/4 = 7/4

*When exponentiated*

\(2^{-2}\) \(2^{-1}\) \(2^{0}\) \(2^{1}\)
1/4 1/2 1 2
5/16 5/8 5/4 5/2
3/8 3/4 3/2 3
7/16 7/8 7/4 7/2
  • Floating point systems can represent a finite number of values.
  • Floating point systems have a minimum and maximum, denoted as \(realmin\) and \(realmax\).
  • The spacing between consecutive numbers increases as the magnitude of the number increases.
  • The spacing between 0 and \(realmin\) is larger than the space between \(realmin\) and the next number.

Double Precision Floating Point

binary64 - specified by IEEE 754

  • 1 sign bit
  • 11 exponent bits
  • 52 fraction bits

    64 bits total

    \(fl(x) = (-1)^{d_s}(1 + d_1 d_2 ... d_52) \cdot 2^{exp - 1023}\)

    This shifted exponent is known as a biased exponent and gives us a range of exponents \([-1023, 1024]\). However, we reserve -1023 and 1024 to represent NaN and Inf, respectively.

    \(realmin = (1.0000...0) \cdot 2^{-1022}\) \(realmax = (1.1111...1) \cdot 2^{1023}\)

Special Numbers

  • \(0(0...0)_2 \cdot 2^{1...1} = +Inf\)
  • \(1(0...0)_2 \cdot 2^{1...1} = -Inf\)
  • \(0(\neq 0)_2 \cdot 2^{1...1} = +NaN\)
  • \(1(\neq 0)_2 \cdot 2^{1...1} = -NaN\)
  • \(0(0...0)_2 \cdot 2^{0...0} = +0\)
  • \(1(0...0)_2 \cdot 2^{0...0} = -0\)

#+endexamples

  1. Find the value of \((7ff0\ 0000\ 0000\ 0000)_{16}\) 0111 1111 1111 0… = +Inf
  2. Find the value of \((4014\ 0000\ 0000\ 0000)_{16} = 0100\ 0000\ 0001\ 0100\ 0\dots\) \(\underbrace{0}_\text{sign bit}\underbrace{100\ 0000\ 0001}_\text{exp}\ \underbrace{0100\ 0...}_\text{frac} = (1+2^{-2}) \cdot 2^{1025 - 1023} = 5\)

Subnormal numbers

\(fl(x) = (-1)^{d_s} (0 + d_1 d_2 d_3 ... d_n)_2 \cdot 2^{exp}\)

  • \(0(\neq = 0)_2 \cdot 2^{0...0} = subnormal\)

Smallest exp - -1074 Largest exp - 1023

Error

Relative error

\(\frac{|x - f(x)|}{|x|}\)

In floating point, the relative error is bounded by a constant \(\eta\), which is given by

\[\eta = \begin{cases} 2^{-t} & \text{if truncated} \\ 2^{-(t+1)} & \text{if rounded} \end{cases}\]

where \(t\) is the number of bits truncated or rounded to.

Derivation

In Octave, exp() gives the distance between a number and the next smallest floating point number.

  1. x = 1 exp(1) = (1.0 ⋯ 1)2 ⋅ 20 = 2-52

Arithmetic Operations

Let \(x,y \in \mathbb{R}\)

Multiplication

\(f(x)\) can be expressed as \(x(1 + \delta)\) where \(|\delta| \leq \eta\)

\(f(x \cdot y) = (x(1 + \delta_x) \cdot y(1 + \delta_y)){(1 + \delta_m)}{rounding after multiplication}\) = xy(1 + δx + δy + δm) + {O(η2)}{≤ c ⋅ η2}

\(\frac{|x \cdot y - f(x \cdot y)|}{|x \cdot y|} = |\delta_x + \delta_y + \delta_m| + O(\eta^2) \leq 3 \eta + O(\eta^2)\)

Addition

\(f(x + y) = \[x(1 + \delta_x) + y(1 + \delta_y)\]\underbrace{(1 + \delta_a)}_{\text{error after rounding addition result}} &= x(1 + \delta_x)(1 + \delta_a) + y(1 + \delta_y)(1 + \delta_a) = x(1 + \delta_x + \delta_a) + y(1 + \delta_y + \delta_a) + O(\eta^2) \\ &= (x + y) + x(\delta_x + \delta_a) + y(\delta_y + \delta_a) + O(\eta^2)\)

\(|\frac{f(x + y) - (x + y)}{x + y}| \leq |\frac{x( \delta_x + \delta_a)}{x + y}| + |\frac{x( \delta_x + \delta_a)}{x + y}|\)

If \(x\) and \(y\) have same sign, \(|\frac{x}{x+y}| \leq 1\) and \(|\frac{y}{x+y}| \leq 1\), then $ ≤|δx + δa| + |δy + δa| ≤ 4 η$

But if \(x\) and \(y\) have different sign, then \(|x - y|\), can be very small and potentialy causes the result to blow up.

We conclude that the result of any floating point arithmetic operation must be equal to the result using infinite precision then rounding to \(t\) binary digits.

  1. Let \(b >> 4ac\), \(x_1 = \frac{-b - \sqrt{b^2 - 4ac}}{2a} \approx \frac{0}{2a}\) \(x_2 = \frac{-b + \sqrt{b^2 - 4ac}}{2a}\)

    To avoid cancellation, we can calculate the roots a different way:

    \(x_1 x_2 = \frac{c}{a} \Rightarrow x_1 = \frac{1}{x_2} \frac{c}{a}\)

  2. Let \(x = 1.1103 \cdot 10^{-1}\), \(y = 9.963 \cdot 10^{-3}\)

    \(x - y = 1.337 \cdot 10^{-3}\)

Subtraction

\(100.0 - 99.99 \Rightarrow 1.000 \cdot 10^2 - 0.9999 \cdot 10^2 = 0.001 \cdot 10^2\)

\(fl(x-y) = .001 \cdot 10^2 = 0.1\) \(\text{relative error} = \frac{0.1 - 0.01}{0.01} = 9\)

with guard digit

\(100.0 - 99.99 \Rightarrow 1.0000 \cdot 10^2 - 0.9999 \cdot 10^2 = 0.01 \cdot 10^2\)

Condition of a problem

condition number \(k_f(x) = |f'(x)| \frac{|x|}{|f(x)|}\)

Proof
  1. \(f(x) = \tan{x}\)

    \(k_f(x) = \frac{|\sec^2 (x)| |x|}{|\tan x|} = \frac{|{x}|}{|\sin (x)||\cos (x)|}\)

    As \(x \rightarrow \pi/2\), \(k_f(x) \rightarrow \infty\).

  2. Let \(y_n = \int_0^1 \frac{x^n}{x + 10} dx\)

    \(y_n\) should monotically decrease as \(n\) increases. However, in computation we'll see that \(y_n\) increases, then become negative.

    \(y_n = g(n) + (-10)^n y_0\)

    \(\frac{dy_n}{dy_0} = (-10)^n\)

    \(k_f(x) = \frac{10^n \cdot |x|}{|f(x)|} = \frac{10^n \cdot |y_0|}{|y_n|}\)

    \(k_f(x) > 10^n\)

Stability of an Algorithm

So far, we've defined two types of error. These errors require that we have the exact answer to the problem, \(y\).

forward error \(\frac{|y - \hat{y}}{|y|}\) absolute forward error \(|y - \hat{y}\)

For some problems, we don't have access to the exact answer, so we instead compute a different kind of error.

backward error \(\frac{|x - \bar{x}|}{|x|}\)

  1. \(fl(x_1 + x_2) = \[x_1(1 + \delta_1) + x_2(1 + \delta_2)\](1 + \delta_3) \\ &= x_1(1 + \delta_1)(1 + \delta_3) + x_2(1 + \delta_2)(1 + \delta_3) = \bar{x_1} + \bar{x_2}\) \(\bar{x_1} = x_1(1 + \delta_1)(1 + \delta_3) = x_1 + x_1(\delta_1 + \delta_3) + O(\delta^2) \rightarrow \frac{|\bar{x_1} - x_1|}{x_1} \leq |\delta_1 + \delta_3| = 2 \eta\)
  2. Let \(y = fl(x_1, x_2) = x_1^2 - x_2^2\) using decimal arith. to 3 digits. Find the backwards error

backwards stability An algorithm is backwards stable if the computed \(\hat{y}\) satisfies

\(\hat{y} = f(x + \delta)\) where \(|\delta| \leq \varepsilon |x|\)

for any \(x\) and a sufficiently small \(\varepsilon\)

\(y_n = \int_0^1 \frac{x^n}{x + 10} dx\) - the thing were trying to calculate

\(y_n + 10y_{n-1} = \frac{1}{n}, n \geq 1\) - the problem \(y_0 = \ln(11) - \ln(10)\)

when we calculate \(y_n\) using floating point arithmetic, this is the algorithm

Finding Roots

\(y_{i+1} = -f(y_i)/f'(y_i)\)

Intermediate Value Theorem

Intermediate value theorem

Let \(f\) be a continuous function in the domain \((a,\ b)\). For any \(y\) s.t. \(f(a) \leq y \leq f(b)\), there exists \(x\) such that \(a \leq x \leq b\)

Secant Method

Convergence

Newton's Method

\(x_{k+1} = x_k - \frac{f(x_k)}{f'(x_k)}\), where \(x_1\) is given. This find the point \(x^*\) such that \(f(x^*) = 0\)

Convergence

  • Newton's method converges quadratically.
Proof
  • Newton's method exhibits local convergence

Bisection Method

Convergence

  • the Bisection method converges globally

Fixed point iteration

This can be interpreted \(x_{k+1} = g(x_k)\) When the solution \(x^*\) is found, \(x^* = g(x^*)\)

There are many functions \(g\) such that \(x^* = g(x^*)\), like \(g(x) = x + \alpha f(x)\)

Let \(g\) be a function such that.

  1. G is continuous on \([a,b]\)
  2. \(g(a) \geq a\), \(g(b) \leq b\)
  3. differentiable on \([a,b]\)
  4. \(|g'(x)| < 1\) on \([a,b]\)

Then

  • there exists a point \(x^* \in [a,b]\) such that \(g(x^*) = x^*\)
  • the stationary point \(x^*\) is unique
  • \(x_{k+1} = g(x_k)\) will converge to \(x^*\) given \(x_1 \in [a,b]\)

Assume \(g(a) - a > 0\) and \(g(b) - b < 0\)

From IVT, there exists some \(x^* \in [a,b]\) such that \(g(x^*) - x^* = 0\)

Uniqueness

Proof

Convergence

Proof

Implementation

Root finding scripts combine multiple methods to achieve speed and reliability in finding roots of \(f\)

  1. check \(f\) at an initial guess \(x_1\) and finds two points \(a,b\) around \(x_1\) such that the signs of \(f(a)\) and \(f(b)\) are different
  2. iterate newtons method a few times until error stops decreasing significantly
    • (for example, iterate until \(\text{curr. err} \leq \frac{\text{prev. err.}}{2} = |x_{k+1} - x_k | \leq \frac{|x_k - x_{k-1}|}{2}\)

Condition Number

We can't use the relative condition number, \(k_r(x) = \frac{|x| \cdot |f'(x)|}{|f(x)|}\), since \(k_r \rightarrow \infty\) as \(f(x) \rightarrow 0\).

Instead we use the absolute condition number to quantify the problem.

absolute condition number

\(k_a(x) = |f'(x)|\)

Proof

\[\log_{10} | \thing{\frac{y - \hat{y}}{y}}{10^{-p}} | = \log_{10} k_r + \log_{10} | \thing{\frac{x - \hat{x}}{x}}{10^{-16}} |\]

\(p = 16 - \log_{10} k_r\)

When finding roots, we find \(x\) given \(y\), so \(x = f^{-1}(y) \rightarrow k_r = \frac{|y| \cdot |f^{-1}'(y)|}{ |f^{-1}(y)| }\).

\(f^{-1}'(y) = \frac{d}{dy} f^{-1} (y) = \frac{dx}{dy} = \frac{1}{dy/dx} = \frac{1}{f'(x)}\)

\(k_{r,f^{-1}} = \frac{|f(x)|}{|x| \cdot |f'(x)|}\)

\(\hat{y} - y = (\hat{x} - x) f'(\psi)\) \(|\hat{y} - y| = |(\hat{x} - x)| |f'(\psi)|\) where \(\psi \in [x,\hat{x}]\)

Approximate \(f'(\psi)\) with \(f'(x)\)

\(|\hat{y} - y| \approx |f'(x)| \cdot |\hat{x} - x|\)

\(|\hat{x} - x| = \frac{1}{|f'(x)|} |\hat{y} - y|\)

\(k_a = \frac{1}{|f'(x)|}\)

  1. Find the roots of \(\sin (x)\) \(y = 0 = x - \sin (x) = f(x)\) \(f'(x) = x - \sin (x)\)

    \(k_a = \frac{1}{|f'(x)|} = \frac{1}{|1 - \cos x|}\)

  2. Find the roots of \((x - 1)(x - 2) \dots (x - 20)\) (Wilkinson's polynomial)

Matrices

\[Ax = b\]

\(x_1 + x_2 = 3\) \(x_1 - x_2 = 1\)

\[\begin{matrix} 1 & 1 \\ 1 - 1 \end{matrix} \cdot \begin{matrix} x_1 \\ x_2 \end{matrix} = \begin{matrix} 3 \\ 1\]

Relative condition: A(x + δ x) = b + δ b

\[\frac{|| \delta x ||}{|| x ||} = k_r \cdot \frac{|| \delta b ||}{|| b ||}\]

\(k_r = || A || \cdot || A^{-1} ||\)

\[|| A ||_1 = max_{1 \leq j \leq n} \sum_{i=1}^n |a_{ij}|\] - max column sum norm \[|| A ||_\infty = max_{1 \leq i \leq n} \sum_{j=1}^n |a_{ij}|\] - max row sum norm

  1. H = hilb(5)
    b = sum(H, 2)
    
    H\b
    
  2. P = pascal(5)
    b = sum(P, 2)
    

Finding Minimums

f = @(x) x.^2 - x - 2; # parabola with minimum at x = .5

x = fminbnd(f, -2, -3)

Assume we are trying to minimize a function \(\phi (x)\), where \(x^*\) is the location of the minimum (minimizer).

Assume \(\phi (x)\) is twice differentiable.

Writing \(\phi\) as a taylor series: \[\phi (x) = \phi (x^*) + (x - x^*) \phi '(x^*) + \frac{(x - x^*)^2}{2} \phi ''(x^*) + O((x - x^*)^3)\]

Since \(x^*\) is a minimizer, \(\phi '(x^*) = 0\), so \[\phi (x) = \phi (x^*) + \frac{(x - x^*)^2}{2} \phi ''(x^*) + O((x - x^*)^3)\]

If \(\phi ''(x^*) > 0\), then \(x^*\) is a minimizer, and if \(\phi ''(x^*) < 0\), then \(x^*\) is a maximizer. If \(\phi ''(x^*) = 0\), there is not enough information to know.

Implementation

A common minimum finder employs two methods:

  1. Golden section search
  2. Parabolic interpolation

Assume \(f\) is a unimodal function (only one minimum).

\(\tau = \frac{\sqrt{5} - 1}{2} \approx 0.618\)

  1. Begin with two points, \(x_1 = a + (1 - \tau)(b - a)\), \(x_2 = a + \tau (b - a)\)

    case 1: if \(f(x_1) < f(x_2)\), then there must be a minimum in the range \([a, x_2]\)

    case 2: if \(f(x_1) > f(x_2)\), then there must be a minimum in the range \([x_1, b]\)

  2. Find two new points \(x_1, x_2\) in the new range, noticing that the unused existing point is already in the right place, ie:

    case 1: \(x_2' = x_1\)

    case 2: \(x_1' = x_2\)

Proof

The new length of the interval after each iteration is \(\tau\) of its previous length.

Parabolic Interpolation

  1. Start with 3 points, \(u, v, w\), such that v is the minimizer: $f(v) < f(u), f(w)$$

  2. Draw a parabola \(p(x)\) that goes through the points \((u,\ f(u))\ ,(v,\ f(v))\ ,(w,\ f(w)))\)

  3. Find the minimum of the parabola, \(m\)

  4. case 1: if \(m < v\), then \(u'=u,\ v'=m\ ,w' = v\)

    case 2: if \(m > v\), then $u

Solving for minimum of parabola

\(L(x) = \frac{f(u)(x - v)(x - w)}{(u-v)(u-w)} + \frac{f(v)(x - u)(x - w)}{v - i)(v - w)} + \frac{f(w)(x - u)(x - v)}{w - u)(w - v)}\)

\(L'(x) = 0\)

\(p = (v - u)^2(f(w) - f(v) - (v - w)^2(f(u) - f(v))\)

Interpolation

Interpolation allows us to:

  • data fitting - if we only have a limited number of sample points, we can find a function that fits the data
  • approximate functions that are difficult to compute - we can compute a function at a handful of points then interpolate to fill in the gaps

Power basis

Given a set of coordinates \((x_i, y_i), \cdots, (x_n, y_n)\), find a function that interpolates the values. We limit the the function to a polynomial of \(n\) degrees to keep the function simple.

Let \(P_n(x) = a_1 x^n + \cdots a_n x + a_{n+1}\) be an nth degree polynomial.

\[\underbrace{\left[ \begin{matrix} x_1^n & \cdots & x_1 & 1 \\ x_2^n & \cdots & x_2 & 1 \\ \vdots & \vdots & \vdots & \vdots \\ x_{n+1}^n & \cdots & x_{n+1} & 1 \end{matrix} \right]}_{V} \left[ \begin{matrix} a_1 \\ a_2 \\ a_3 \\ a_3 \end{matrix} \right] = \left[ \begin{matrix} y_1 \\ y_2 \\ \vdots \\ y_{n+1}\end{matrix} \right]\]

If \(x_i \neq x_ x_j\), then \(Va = y\) has a unique solution. \(a\) contains the coefficients of the interpolation polynomial.

The basis here is \({x^n, x^{n-1}, \cdots, x, 1}\) and is known as a power basis.

#+beginexamples

  1. (1,16), (3,21), (5,15), (15,12). Fit to a 3rd degree polynomial: \(P_3(x) = a_1 x^3 + a_2 x^2 + a_3 x + a_4\)

    \[\left[ \begin{matrix} 1 & 1 & 1 & 1 \\ 27 & 9 & 3 & 1 \\ 125 & 25 & 5 & 1 \\ 216 & 36 & 6 & 1 \end{matrix} \right] \left[ \begin{matrix} a_1 \\ a_2 \\ a_3 \\ a_3 \end{matrix} \right] = \left[ \begin{matrix} 16 \\ 21 \\ 15 \\ 12\end{matrix} \right]\]

    \(V\) is the Vander Monde matrix, and if \(x_i \neq x_j\), then \(V\) is non-singular

#+endexapmles

However Vander Monde matrices are known for being ill-conditioned, so this is not a good method for computing interpolating polynomials.

General matrix form

Let the basis functions be \(\phi_1 (x), \phi_2 (x), \cdots, \phi_{n+1} (x)\)

\(y = c_1 \phi_1 (x) \cdot + \phi_{n+1} (x)\)

\[\left[ \begin{matrix} \phi_1 (x_1) & \cdots & \phi_n (x_1) & \phi_{n+1} (x_1) \\ \phi_1 (x_2) & \cdots & \phi_n (x_2) & \phi_{n+1} (x_2) \\ \vdots & \vdots & \vdots & \vdots \\ \phi_1 (x_{n+1}) & \cdots & \phi_n (x_{n+1) & \phi_{n+1} (x_{n+1}) \end{matrix} \right] \left[ \begin{matrix} a_1 \\ a_2 \\ a_3 \\ a_3 \end{matrix} \right] = \left[ \begin{matrix} y_1 \\ y_2 \\ \vdots \\ y_{n+1}\end{matrix} \right]\]

Lagrange basis

\[L_i(x) = \frac{(x - x_1)(x - x_2) \cdot (x - x_{i-1}) (x - x_{i+1}) \cdot (x - x_{n+1})}{(x_i - x_1) (x_i - x_2) \cdots (x_i - x_{i-1}) (x_i - x_{i+1}) \cdots (x_1 - x_{n+1})}\]

\(L_i\) has the property that \[L_i(x_j) = \begin{cases} 1 & i = j \\ 0 & i \neq j \end{cases}\]

\(L(x) = y_1 L_1 (x) + y_2 L_2 (x) + \cdots + y_{n+1} L_{n+1} (x)\)

Each \(L_i (x)\) is an nth degree polynomial.

\[L(x) = \sum_1^{n+1} y_i L_i (x)\]

Where \[L_i(x) = \frac{\prod_1^{n+1} (x - x_i)}{\prod_1^{n+1} (x - x_j)}\]

Let \(\psi(x) = \prod_{j = 1}^{n + 1} \frac{f(x_i)}{ x - x_j} w_i\)

\[w_i = \frac{1}{\prod_{j = 1,\ j \neq i}^{n+1} (x_i - x_j)}\]

$Li(x) =

\[L = \sum_i=1^{n+1} y_i \frac{\phi (x)}{x - x} \cdot w_ifdA \]L9x) ϕ(x) \

\(L(x) = \frac{\sum_{i=1}^{n+1} \frac{z_i x - x_i}}\)

Newton Form Interpolating Polynomial

\(P_n(x) = f[x_1] + \underbrace{f[x_1, x_2]}_{\text{divided difference}}(x - x_1) + \dots + f[x_1, x_2, \dots, x_n](x - x_1) \dots (x - x_n)\)

\[P_n(x) = \sum_{i=1}^{n+1} \left[ f[x_1, \cdots, x_i] \prod_{j=1}^{i=1} (x - x_j) \right]\]

The divided difference approaches the derivative if all \(x_i\) begin to approach the same point: \[\lim_{(x_0, \cdots, x_n) \rightarrow (x, \cdots, x)} f[x_0, \cdots, x_n] = f^{(n)}(x)\]

So the Newton form polynomial is equivalent to a Taylor polynomial if all \(x_i\) begin to approach the same point: \[\lim_{(x_0, \cdots, x_n) \rightarrow (x, \cdots, x)} f[x_1] + f[x_1, x_2](x - x_1) + \dots + f[x_1, x_2, \dots, x_n](x - x_1) \dots (x - x_n) = f(x) + (\varepsilon - x) f'(x) + (\varepsilon - x)^2 \frac{f(x)}{2!} + \cdots + (\varepsilon - x)^n \frac{f(x)^(n)}{(n!)}\]

Divided difference

\[f[x_0, \cdots, x_n] = \frac{f[x_1, \cdots, x_n] - f[x_0, \cdots, x_{n-1}]}{x_n - x_0}\]

where

\[f[x_0, x_1] = \frac{f[x_1] - f[x_0]}{x_1 - x_0}\]

\[f[x_1, x_2] = \frac{f[x_2] - f[x_1]}{x_2 - x_1}\]

\[f[x_0, x_1, x_2] = \frac{f[x_1, x_2] - f[x_0, x_1]}{x_2 - x_0}\]

\[\cdots\]

Interpolation with Derivatives

\[f[x_1, x_2] = \frac{f[x_2] - f[x_1]}{x_2 - x_1}\]

\[\lim_{x_2 \rightarrow x_1} f'(x_2)\]

so define \(f[x_1, x_1] = f'(x_1)\),\(\underbrace{f[x_1, x_1,\dots , x_1]}_{\text{k+1 times}} = \frac{f^{(k)}}(x_1){k!}\)

Given the points \(t_1, t_2, \cdots, t_q\) and multiplicities \(m_1, m_2, \cdots, m_q\)

we can find an inteprolating polynomial such that \(p^{(k)}(x_i) = f^{(k)}(x_i)\), for \(k \in [0,m_i]\)

We have \(\sum_{i=1}^q (m_i + 1) = \sum_{i= 1}^q (m_i + q)\) constraints.

  1. Assume we want to find a polynomial that interpolates the points, \(f(0) = 1\), \(f'(x) = -1\), \(f(1) = 1\)

    \(x_i\) \(f[x]\) \(f[x,x]\) \(f[x,x,x]\)
    \(x_0 = 0\) 1
    \(x_0 = 0\) 1 \(f[x_0, x_0] = f'(x_0) = -1\)
    \(x_1 = 1\) 1 \(f[x_0, x_1] = \frac{f[x_1] - f[x_0]}{x_1 - x_0} = 0\) \(f[x_0, x_0, x_1 = \frac{f[x_0, x_1] - f[x_0, x_0]}{x_1 - x_0} = 0\)

    Interpolating polynomial: \(1 - 1 \cdot(x-0) + 1 \cdot(x - 0)^2\) = x2 - x + 1$

  2. Find a polynomial such that \(f(0) = -1, f'(0) = 1, f(1) = 0, f'(1) = 2\)

    \(x_i\) \(f[x]\) \(f[x, x]\) \(f[x, x, x]\) \(f[x, x, x, x]\)
    \(t_1 = 0\) -1
    \(t_1 = 0\) -1 f[t1, t2] = f'(t1) = -1$
    \(t_2 = 1\) 0 \(f[t_1, t_2] = \frac{f[1] - f[0]}{1 - 0} = 1\) \(f[t_1, t_1, t_2] = \frac{1 - 1}{1 - 0} = 0\)
    \(f[t_2, t_2] = f'(t_2) = 2\) \(f[t_1, t_2, t_2] = \frac{f[t_2, t_2] - f[t_1, t_2]}{t_2 - t_1} = 1\)
    \(f[t_1, t_1, t_2, t_2] = \frac{f[t_1, t_2, t_2] - f[t_1, t_1, t_2]}{t_2 - t_1} = 1\)

    \(P_3(x) = - 1+ 1 \cdot(x - 0) + (x - 0)^2 + 1 \cdot (x - 0)^2 (x - 1) = x^3 - x^2 + x - 1\) \(P_3'(x) = 3x^2 - 2x + 1\)

    We can verify this polynomial at various values

    x \(P_3(x)\) \(P_3'(x)\)
    0 -1 1
    1 0 2

    So this meets the original constraints.

Error

\[e_n(x) = |f(x) - P_n(x)| = \frac{f^{(n+1)}(\zeta) }{(n+1)!} \phi(x)\], where \(\phi(x) = (x - x_1) \cdots (x - x_{n+1})\) and \(\frac{f^{(n+1)}(\zeta)}{(n+1)!} = f[x_1, x_2, \cdots, x_{n+1}, x]\)

for \(n=1\), \(P_1(x) = f[x_1] + f[x_1, x_2](x - x_1)\).

  1. \[e_1(x) = |f(x) - P_1(x)| = |\frac{f''(\zeta)}{2}| \cdot |(x - x_1)(x - x_2)| \leq \frac{1}{2} \text{max} |f''(\zeta)| \cdot \text{max}_{x \in[x_1, x_2]} | (x - x_1)(x - x_2)|\]

    where the max of \((x - x_1)(x - x_2)\) is \(x = \frac{x_1 + x_2}{2}\) which evaluates to \(\frac{(x_2 - x_1)^2}{4}\)

    so \(e_1(x) \leq \frac{1}{2} \cdot \frac{(x_2 - x_1)^2}{4} \cdot \text{max}_{\zeta \in [x_2, x_2]} f''(\zeta)\)

Error in the polynomial is driven by two terms:

  • the derivative of f
  • placement of interpolated points - equally spaced points usually drive this term high and lead to a large error. See chebyshev polynomial for information about ideal point placement.

Horner's Method

\(P = c_{n+1} (x - x_n)\)

\(P = \left[ P + c_n \right] (x - x_{n-1})\)

$P = [P+cn − 1] (x - xn-2)

jth iteration

\(P = \left[ P + c_j \right] (x - x_{j-1})\)

Pseudo code

Chebyshev Interpolating points.

\(n \theta = (2k - 1) \frac{\pi}{2}\), for \(k \in [1, \cdots, n]\)

\(\theta = \frac{2k - 1}{2n} \pi\)

Piecewise Interpolation

Hermite Interpolation

Assume we are interpolating a function \(f(x)\). We can perform nth degree Hermite interpolation at \((x_1, x_2, \cdots, x_n)\) if we know \(f(x_i), f'(x_i), \cdots, f^{(n)}(x_i)\).

Numerical Integration

Suppose we want to integrate

\[I = \int_a^b f(x) dx\]

We can replace \(f(x)\) by an interpolating polynomial and solve numerically.

\[Q = \int_a^b P_n(x) dx\]

where the error is \[E = \int_a^b f(x) - P_n(x) dx\]

Midpoint Method

Let \([a, b]\) be a small interval.

Approximate \(f\) over \([a,b]\) with \(f(\frac{a+b}{2})\).

\(Q = \int_a^b f\left( \frac{a+b}{2} \right) dx = f \left(\frac{a+b}{2} \right) (b - a)\)

This is equivalent to using an interpolating polynomial of degree 0.

midpoint method drawing

Error

\(f(x) = f(m) + f'(\varepsilon) \cdot (x - m) + \frac{f''(\varepsilon)}{2!}(x - m)^2\), where \(\varepsilon \in [a,b]\)

The first part of the error term is

\(\int_a^b f'(m) (x - m) dx = 0\)

The second part of the error term is

\(\text{error} = \frac{1}{2} \int_a^b f''(\varepsilon) (x - m)^2 dx\)

Since \((x-m)^2\) doesn't change sign, we can use the mean value theorem for integrals to conclude,

\(\text{error} = f''(\eta) \frac{1}{2} \int_a^b (x-m)^2 dx = \frac{f''(\eta) (b-a)^2}{24}\), for \(\eta \in [a,b]\)

MVT for Integrals

Assume \(g(x) \in C[a,b]\). \(\phi(x)\) is integrable that is nonnegative or nonpositive on \([a,b]\). Then there is a point \(\eta \in [a,b]\) such that \(\int_a^b g(x) \phi(x) dx = g(\eta) \inta^b \phi(x) dx\)

Trapezoidal Method

Let \([a, b]\) be a small interval.

Approximate \(f\) over \([a,b]\) with \(f(a) + f[a,b](x - a)\).

\[\begin{aligned} Q &= \int_a^b \left[ f(a) + f[a,b] (x - a) \right] dx \\ &= f(a)(b - a) + f[a,b] \frac{(x - a)^2}{2}_\eval{a,b} \\ &= \frac{b-a}{2} \cdot \left[ f(a) + f(b) \right] \end{aligned}\]

trapezoidal method drawing

Error

\(f(x) = f(a) + f[a,b](x - a) + \frac{f''(\varepsilon)}{2} (x-a)(x -b)\)

\(\text{error} = \int_a^b \frac{f''(\varepsilon)}{2} \underbrace{(x - a)}_{\geq 0}\underbrace{((x - b)}_{\leq 0}\), where \((x - a)(x - b)\)

Since \((x - a)(x - b)\) is nonpositive,

\[\text{error} = \frac{f''(\eta)}{2} \int_a^b (x - a)(x - b) dx = -f''(\eta)}\frac{(b - a)^3}{12}\]

Simpson's Method

Replace \(f(x)\) by a 2nd order polynomial in Lagrange form.

\[Q = \int_a^b P_2(x) dx = \frac{h}{6} \left[ f(a) + 4 f(m) + f(b) \right]\]

Error

\(\text{error} = \frac{f^{(4)}(\eta) (b - a)^5}{90 \cdot 2^5}\)

Degree of Precision

The degree of precision

Composite Midpoint Rule

Sum

\[M = \frac{b - a}{n} \sum_{i=1}^n f \left( \frac{ x_i + x_{i+1}}{2} \right)\]

Proof

Error

\[\text{error} = \frac{h^2}{24}(b - a) \cdot f''(\eta)\]

for \(a \leq \eta \leq b\)

Proof

Composite Trapezoidal Rule

Sum

\[Q = \frac{h}{2} \left[ f(a) + 2 \sum_{i=1}^n f(x_i) + f(b) \right]\]

Proof

Error

\(\text{error} = \frac{-h^2}{12}(b-a) f''(n)}\)

where \(a \leq \eta \leq b\)

Proof

Composite Simpson's Rule

Sum

Assume equal subintervals with endpoints \([x_1, \cdots, x_{n+1}]\)

\[S = \frac{h}{3} \left[ f(x_1) + 4 \sum_{i=2,4, ...}^n f(x_i) + 2 \sum_{i=3,5, ...}^n f(x_i) + f(x_{n+1}) \right]\]

Proof

Error

\[\text{error} = \frac{h^4}{180}(b - a) f^{(4)}(\psi)\]

where \(a \leq \psi \leq b\)

Proof

Gaussian Quadrature

Assume we are computing \(I = \int_{-1}^1 f(x) dx\). With Gaussian quadrature method, we use carefully chosen points \([x_1, \cdots, x_{n+1}] \in [-1, 1]\) to make the integral precision as high as possible.

We begin by defining a family of orthogonal polynomials known as Legendre polynomials

Legendre polynomial

\(P_0(x) = 1\)

\(P_1(x) = x\)

\(P_2(x) = \frac{1}{2}(3x^2 - 1)\)

\(P_3(x) = \frac{1}{2} (5x^3 - 3x)\)

\(P_4(x) = \frac{1}{8}(35x^4 - 30x^2 + 3)\)

\(P_5(x) = \frac{1}{8}(63x^5 - 70x^3 + 15x)\)

\[P_{j+1}(x) = \frac{2j + 1}{j+1} x P_j(x) - \frac{j}{j+1} P_{j-1}(x)\]

These polynomials have the following properties:

  • j is odd -> \(P_j\) is even
  • j is even -> \(P_j\) is odd
  • \(P_j(1) = 1\)
  • \(\int_{-1}^1 P_k(x) x^j dx = 0\) for \(j = 0, \cdots, k-1\)

We use the roots of \(P_k(x)\) as the interpolation points for quadratic interpolation.

Proof
Proof

Approximating Functions

Given a function \(f\) on an interval \([a, b]\) and a set of basis functions, \(\phi_j\) for \(j \in 0 \cdots n\),

we try to find an approximation \[v(x) = \sum_{j=0}^n c_j \phi_j (x)\]

such that the following quantity is minimized

\[\int_1^b [f(x) - v(x) ]^2 dx = || f - v ||_2^2\]

We denote the minimizing \(v\) as \(v^* = \sum_0^n c_j^* \phi_j(x)\)

\[min ||f - v||_2^2 = \int_a^b \left( f(x - \sum_{j=0}^{n+1} c_j \phi_j(x) \right)^2 dx\]

Differentiate with respect to \(c_k\)

\[\frac{d}{dk} \int_a^b \left( f(x - \sum_{j=0}^{n+1} c_j \phi_j(x) \right)^2 dx = 2 \int_a^b \left( f(x - \sum_{j=0}^{n} c_j \phi_j(x) \right) ( - \phi_k(x) ) dx = 0\]

\[\int_a^b f(x)(- \phi_k(x)) dx = \int_a^b \sum_{j=0}^{n} \ c_j phi_j(x) \phi_k(x) dx = \sum_{j=0}^n c_j \int_a^b \phi_j(x) \phi_k(x) dx\]

so

\[\underbrace{\sum_{j=0}^n c_j \int_a^b \phi_j(x) \phi_k(x) dx}_{\sum_{j=0}^n c_j B_{kj}} = \underbrace{\int_a^b f(x)(- \phi_k(x)) dx}_{b_k}\]

If we choose the monomial basis,

\(c_i = \frac{b_i}{B_{ii}}\)

where \(b_i = \int_a^b w(x) f (x) \phi_i(x) dx\)

Let the basis functions be the Legendre polynomials

\(\phi_0 = 1\), \(\phi_1 = x\), \(\phi_2 = \frac{1}{2}(3x^2 - 1)\)

If we are minimizing \(||f - v(x) ||_2^2\)

then \(v^*(x) = \frac{e - e^{-1}}{2} + 3e^{-1}x + \frac{5}{2}(e - 7e^{-1}) \frac{1}{2}(3x^2 - 1)\)

\(c_1 = \frac{b1}{B_{11}}\), \(b_1 = \int_{-1}^1 e^x \phi_1(x) d = 2e^{-1}\), \(B_{11} = \frac{2}{3}\)

\(c_2 = \frac{b2}{B_{22}}\), \(b_2 = \int_{-1}^1 e^x \phi_2(x) d = e - 7e^{-1}\), \(B_{22} = \int_{-1}^1 \left(\frac{1}{2}(3x^2 - 1) \right)^2 dx = \frac{2}{5}\)

Properties of Norms

  1. \(|| g || \geq 0\), \(||g|| = 0 \Rightarrow g(x) = 0\) on \([a,b]\)

  2. \(||\alpha g|| = |\alpha|||g||\)

  3. \(||f + g||_2^2 \leq \left( ||f||_2 + ||g||_2 \right)^2\)

Least Squares

\[Va \approx y\]

We can get the residual error by $r = y - Va\[ Since $Va^*$ is a projection of $y$ to a lower dimension, $r$ is orthogonal to $Va^*$ and therefore $V^T$. $V^T(y - Va^*) = 0 \Longrightarrow V^Ty = V^TVa$ ** Orthogonal Polynomials $\phi_{-1} = 0$ $\phi_0 = 1$ \]ϕj(x) - x ϕj-1(x) = b ϕj-1 - c ϕj-1} + d ∑k=0j-3\[ three term recurrence relation: $\phi_j(x) = (x - b_j) \phi_{j-1}(x) - c_j \phi_{j-2}(x)$ $(\phi_j, \phi_{j-1}) = 0 = (x \phi_{j-1}, \phi_{j-1}) - \b_j (\phi_{j-1}, \phi_{j-1}) - c_j (\phi_{j-2}, \phi_{j-1})$ Multiply by $\phi_{j-1}$, integrate from a to b with weight for $w$ to find $b_j$ \]bj = = \[ Repeat with $\phi_{j-2}$ to find $c_j$ \]cj = \[ ** Legendre Polynomial $w(x) = $ $\phi_1(x) = (x - b_1) \phi_0 - c_1 \phi_{-1} = (x - b_1) 1 - c_1 0 = x - b_1$ From our earlier results, \]bj = = \[ \]c2 = = = = 1/3\[ \]ϕ2 = x2 - 1/3\[ *** Orthogonality \]Pn(x) = (1 - x2)n\[ $P_0(x) = 0$, $P_1(x) = 1$ We will show $\int_{-1}^1 P_j(x) P_k(x) dx = 0$ where $j < k$ \]-11 (1 - x2)j ⋅ (1 - x2)k dx

Chebyshev Polynomial

Family of orthogonal polynomials on the interval \([-1, 1]\) with weights \(w(x) = \frac{1}{\sqrt{1 - x^2}}\).

i.e. the weights approach infinity when \(x \rightarrow \pm 1\), which puts more weight at the ends of the interval.

\(\phi_0 = 1\), \(\phi_1 = x\), \(\phi_j = (x - b_j) \phi_{j-1} - c_j \phi_{j - 2}\)

\(\phi_2 = x^2 - 1/2\), \(\phi_3 = x^3 - 3/4\), \(\phi_4 = x^4 - x^2 + 1/8\)

Alternative Representation

\(T_j = \cos(j \arccos(x)) = \cos(j \theta)\), where \(\theta = \arccos(x)\)\(\cos(\theta) = x\)

which has roots \(j \theta = (2k - 1) \frac{\pi}{2} \Rightarrow \theta_k = \frac{2k - 1}{j} \frac{\pi}{2}\) for \(k = 1,2, \cdots, j\)

Notice that \(T_j\) has the same roots as \(\phi_j\), but scaled.

\(T_2 = 2x^2 - 1\) \(\phi_2 = x^2 - 1/2\)

The function has extreme points at \(\cos(j \theta) = \pm1 \Rightarrow \theta = \frac{k \pi}{j}\)

\(T_j = \cos(j \theta)\)

\(T_{j+1} = \cos((j + 1) \theta) = \cos(j\theta + \theta) = \cos(j \theta) \cos( \theta) - sin(\j \theta) \sin(\theta)\)

\(T_{j-1} = \cos((j - 1) \theta = \cos(j \theta - \theta) = \cos(j \theta) \cos( \theta) + sin(j \theta)\sin(\theta)\)

\(T_{j+1} + T_{j - 1} = 2 \cos(j \theta) \cos(\theta) = 2 T_j (\theta) \cdot x = 2x t_j(x) - T_{j-1}(x)\)

\(T_0 = 1\)

\(T_1 = x\)

\(T_2 = 2x^2 - 1\)

\(T_3 = 2x(2x^2 - 1) - x = 4x^3 - 3x\)

\(\phi_j = \prod_{k=1}^j (x - \cos(\theta_k))\) for \(\theta_k = \frac{2k - 1}{j} \frac{\pi}{2}\)

\(\phi_j = \frac{1}{2^{j-1}} T_j(x)\)

The minimum of \(max_{-1 \leq x \leq 1} |P_n(x)\) is \(max_{-1 \leq x \leq 1} |\phi_n(x)| = \frac{1}{2^{n-1}}\)

Proof

todo

  • algo vs problem
  • problem conditioning
  • backwards stability
  • find k of integration prob
  • types of convergence
    • linear
    • quadratic
    • forbinius
    • 3.5
  • hessian matrix
  • read chapter 10
  • review parabolic interpolation
  • piecewise polynomial interpolation

  • chebyshev polynomial ch12

  • everything up to piecewise polynomial interpolation

  • read chebyshev polynomial 10.6