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
*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 |
binary64 - specified by IEEE 754
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}\)
#+endexamples
\(fl(x) = (-1)^{d_s} (0 + d_1 d_2 d_3 ... d_n)_2 \cdot 2^{exp}\)
Smallest exp - -1074 Largest exp - 1023
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.
In Octave, exp() gives the distance between a number and the next smallest floating point number.
Let \(x,y \in \mathbb{R}\)
\(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)\)
\(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.
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}\)
Let \(x = 1.1103 \cdot 10^{-1}\), \(y = 9.963 \cdot 10^{-3}\)
\(x - y = 1.337 \cdot 10^{-3}\)
\(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 number \(k_f(x) = |f'(x)| \frac{|x|}{|f(x)|}\)
\(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\).
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\)
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|}\)
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
\(y_{i+1} = -f(y_i)/f'(y_i)\)
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\)
\(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\)
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.
Then
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\)
Root finding scripts combine multiple methods to achieve speed and reliability in finding roots of \(f\)
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)|\)
\[\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)|}\)
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|}\)
Find the roots of \((x - 1)(x - 2) \dots (x - 20)\) (Wilkinson's polynomial)
\[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
H = hilb(5)
b = sum(H, 2)
H\b
P = pascal(5)
b = sum(P, 2)
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.
A common minimum finder employs two methods:
f = @(x) x.^2 + 4 * cos(x);
fplot(f, [0 3]);
options = optimset('Display', 'iter');
[xmin, fmin] = fminbnd(f, 0,3, options);
octave:4> [xmin, fmin] = fminbnd(f, 0,3, options);
Func-count x f(x) Procedure
1 1.85410 2.319570 initial
2 2.29180 2.611785 golden
3 1.83004 2.323648 parabolic
4 1.88882 2.316881 parabolic
5 1.89619 2.316809 parabolic
6 1.89554 2.316808 parabolic
7 1.89549 2.316808 parabolic
8 1.89549 2.316808 parabolic
9 1.89549 2.316808 parabolic
10 1.89549 2.316808 parabolic
11 1.89549 2.316808 golden
Optimization terminated:
the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-08
Assume \(f\) is a unimodal function (only one minimum).
\(\tau = \frac{\sqrt{5} - 1}{2} \approx 0.618\)
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]\)
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\)
The new length of the interval after each iteration is \(\tau\) of its previous length.
Start with 3 points, \(u, v, w\), such that v is the minimizer: $f(v) < f(u), f(w)$$
Draw a parabola \(p(x)\) that goes through the points \((u,\ f(u))\ ,(v,\ f(v))\ ,(w,\ f(w)))\)
Find the minimum of the parabola, \(m\)
case 1: if \(m < v\), then \(u'=u,\ v'=m\ ,w' = v\)
case 2: if \(m > v\), then $u
\(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 allows us to:
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,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.
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]\]
\[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}}\)
\(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!)}\]
\[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\]
\[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.
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$
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.
\[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)\).
\[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:
\(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
\(n \theta = (2k - 1) \frac{\pi}{2}\), for \(k \in [1, \cdots, n]\)
\(\theta = \frac{2k - 1}{2n} \pi\)
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)\).
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\]
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
\(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]\)
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\)
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
\(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}\]
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]\]
\(\text{error} = \frac{f^{(4)}(\eta) (b - a)^5}{90 \cdot 2^5}\)
The degree of precision
Sum
\[M = \frac{b - a}{n} \sum_{i=1}^n f \left( \frac{ x_i + x_{i+1}}{2} \right)\]
Error
\[\text{error} = \frac{h^2}{24}(b - a) \cdot f''(\eta)\]
for \(a \leq \eta \leq b\)
Sum
\[Q = \frac{h}{2} \left[ f(a) + 2 \sum_{i=1}^n f(x_i) + f(b) \right]\]
Error
\(\text{error} = \frac{-h^2}{12}(b-a) f''(n)}\)
where \(a \leq \eta \leq b\)
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]\]
Error
\[\text{error} = \frac{h^4}{180}(b - a) f^{(4)}(\psi)\]
where \(a \leq \psi \leq b\)
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:
We use the roots of \(P_k(x)\) as the interpolation points for quadratic interpolation.
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}\)
\(|| g || \geq 0\), \(||g|| = 0 \Rightarrow g(x) = 0\) on \([a,b]\)
\(||\alpha g|| = |\alpha|||g||\)
\(||f + g||_2^2 \leq \left( ||f||_2 + ||g||_2 \right)^2\)
\[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
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\)
\(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}}\)
piecewise polynomial interpolation
chebyshev polynomial ch12
everything up to piecewise polynomial interpolation
read chebyshev polynomial 10.6