Material delivered by both zoom and online.

Fridays are lab assignments, homework also due. May have midterms

Final project exists, will be described later, due last week of class

Found in Chapter 1 of text book

Numerical analysis – study of algorithms used for numerical approx. for problems of mathematical analysis

Mathematical analysis – limits, differentiation, integration, infinite series, etc

Approximation is particularly important, and numbers represented as 32(64) bits, a sign, exponent, mantissa

• This approximation is innaccurate for many numbers

Representations of numbers can be inaccurate – there may be round-off error

• Numbers between two representations may only be approximated!

• Round – use the nearest representable value

• Chop – use the representable part of the value

Function approximation by series

• Some functions are representable by series, $e^{x}$, $\sin x$, etc.

• These may be represented as Taylor polynomials, which are rather generally useful.

Basic method in Numerical Analysis

Builds on Intermediate value theorem

• Given function continuous from $[a, b]$, and $K$ in $[a, b]$, then there exists $c$ such that $f(c) = K$

• When $f(a)$ and $f(b)$ have opposite signs, then there exists $p$ between $(a,b)$, $f(p) = 0$, that is, it has a root in the interval

Repeated Application of above

• Assume continuity, and $f(a)$, $f(b)$ have opposite signs

• f, a, b

• $p = (a + b)/2$

• if stopping criteria: return $p$

• If $f(a)$ and $f(p)$ have opposite signs, then recurse with $f, a, p$

• Else, recurse with $f, p, b$

Tolerance Criteria

• $TOL$ is error tolerance value

• Stop when one of the following:

• $|p_{n} - p_{n - 1}| < TOL$ (absolutely tiny difference)

• $|p_{n} - p_{n - 1}|/|p_{n}| < TOL$ (normalized difference)

• $|f(p_{n})| < TOL$ (very close to zero)

Characteristics of the Method

• Disadvantages

  • Slow to converge, as N may become large before $p_{n} - p_{n - 1}$ becomes small

  • Possible that intermediate approx may be inadvertently discarded!

Can also find fixed points

• Let $g(x)$ be a function, $p$ is a fixed point if $g(p) = p$

• Transformation requires finding root $f(x) = g(x) - x$

Interpolation, used for estimation and prediction

• Census data, pretty well known, overall population of US

• But 10-year cycle, what about 1975, 2005, 2015?

• Requires interpolation

• This estimates/predicts, by finding a function that fits the data

Can be done in higher dimensions!

• Interpolating from, for instance, a set of elevation points, provides a smoother, more realistic curve

Other interpolation methods

Matlab!

Equations easily defined

Each script is separate function (gross)

High-degree polynomials fail

• Sometimes data doesn't fit polynomials!

• Runge's Phenomenon – Edges of interval oscilated at polynomial interpolation with high degree

Splines used – spline function consists of polynomial pieces joined with smoothness conditions

Piecewise, works for $t_{0}$ to $t_{1}$, $t_{1}$ to $t_{2}$, \ldots, $t_{{n - 1}}$ to $t_{n}$

• $t$ are called "knots", ordered

• $\mathcal{S}$ is spline of degree $k$ if

  • each polynomial is degree $k$

  • and is $k - 1$ times continuously differentiable

  • Note, end points must match

  • k = 1, linear

  • k = 2, quadratic

  • k = 3, cubic (most used)

Quadratic:

• Must find a, b, c

• Can be tough, given matching condition

• Boundary equation can also be added, $Q_{0}^{{\prime}}(t_{0}) = 0$

• Recurrence equation:

  • $z_{0} = 0$

  • $z_{{i + 1}} = -z_{i} + 2\left(\frac{y_{i + 1} - y_{i}}{t_{i+1} - t_{i}}$

  • $Q_{i}(x) = \frac{z_{{i + 1}} - z_{i}}{2(t_{i + 1} - z_{i})} (x - t_{i})^{2} + z_{i}(x - t_{i}) + y_{i}$

Most complex yet, but frequently used

Given $n$ data points $[(x_{i}, a_{i})]$, find $b_{j}$, $c_{j}$, $d_{j}$ coefficients for $S_{j}(x)$ for $0 \leq j \leq n - 1$

• Find intervals matching at knots, and assuming normal boundary conditions

• Then, $h_{j}$ is the difference between two knots, $x_{j + 1 }- x_{j}$

  • Then $h_{j-1 }c_{{j - 1}} + 2(h_{{j-1}} + h_{j})c_{j} + h_{j} c_{j+1} = \frac{3}{h_{j}}(a_{j+1} - a_{j}) - \frac{3}{h_{j-1}}(a_{j} - a_{j-1})$

  • $b_{j} = \frac{1}{h_{j}}(a_{{j+1}} - a_{j}) - \frac{h_{j}}{3}(2c_{j} + c_{j+1})$

  • $d_{j} = \frac{1}{3h_{j}}(c_{{j + 1}} - c_{j})$

  • $S_{j}(x) = a_{j} + b_{j}(x - x_{j}) + c_{j} (x - x_{j})^{2} + d_{j} (x - x_{j})^{3}\ldots$

  • Can be represented as matrix, $A$ of $h$ values, $b$ as vector, and $c$ as vector

  • $\mathbf{b}$ is vector of constants (from $b_{j}$), $\mathbf{x}$ is $c$

  • Solve $A\mathbf{x} = \mathbf{b}$

Not particularly difficult – matrix solution possible, trivial even

Cubic Spline most frequently used for single-variable functions

Solved using diagonal/triangular matrix-vector equation

$r$ cubic spline interpolation based on recurrence equations

$f(t)$ from $\mathbb{R} \rightarrow \mathbb{R}$

$S_i(t) = a_i + b_i(t - t_i) + c_i(t - t_i)^2 + d_i(t - t^i)^3$

• $a$ are knots

• Smoothness conditions, require shared knots

• first derivative sharing

• second derivative sharing

Equations should show this

$h_t = t_{{i + 1}}_{} - t_{i}$

• $a_{i} + b_{i} h_{i} + c_{i }h_{i}^{2} + d_{i} h_{i}^{3} = a_{{i + 1}}$

• Finally, we can solve for d, b, a

• Can follow on and continue to solve

• Can eventually get equations only in a, h, i (15)

Then, Can be found using matrix equations, and recurrence equation form can be used

Note, assumes equal spacing

Consider fairly flexible

• Normally natural cubic spline interpolation

We'll create triangles for our sample points, using known algorithm

Gets a triangulated irregular network

Will use a number of triangles in other triangles

Coefficients used are shape functions, half the determinant of a matrix

Project all points to plane

Triangulate all points, i.e., connect all points to create triangles

Then, use triangulation, extend to triangulated irregular network, add back dimension in Z

Elements of TIN are a set of functions w, which divide area of triangle

• $w(x, y) = N_{1}(x,y)w_{1} + N_{2}(x, y)w_{2} + N_{3}(x, y)w_{3}$

• N are based on various things, including $\mathcal{A} = \frac{1}{2} \mathrm{det}\begin{bmatrix}1 & x_{1} & y_{1} \\1 & x_{2} & y_{2} \\1 & x_{3} & y_{3}\end{bmatrix}$

Interpolation function must be calculated carefully

This presents the finite element method

• Triangulations and shape functions are significant part of method

• Interpolation provides a method of integrating approximations of complex functions that can't be directly integrated

Interpolation in Matlab

• scatteredInterpolant uses the same kind of methods, specify as 'linear'

• Takes x, y, z

• Can generate a grid, and then visualise using mesh

Test Overview

• Root finding

  • Bisection, Newton's, Convergence

• Single-variable interpolation

  • Lagrange

  • Divided Differences

  • Piecewise Interpolation: Linear, Quadratic, Cubic Spline

• No spatial interpolation, no Matlab

Use textbook for quadratic/cubic spline interpolation

Approximating derivatives using polynomial approximating the function

Limit Theorem ($f^{{\prime}}(x_{0}) = \lim_{{h \to 0}}\frac{f(x_{0} + h) - f(x_{0})}{h}$

• $h$ should be as small as possible for achiving good accuracy

• Can have issues like subtraction of numbers with nearly equal values

• Errors can be amplified (division with tiny denominator)

• Just use it by dividing by $h$

Taylor's theorem gives us a definition of $f(x + h)$

• Involves an error term, $\zeta$ ($f^{\prime}(x) = \frac{1}{h}(f(x + h) - f(x)) - \frac{1}{2!}hf^{\prime\prime}(\zeta)$, $\zeta \in [x, x + h]$)

• Looking backward adds error term, rather than subtract, linear function of $h$

Interpolation error theorem

• Generally, for an interpolated function, a zeta function is given.

Three-point formulas used instead:

• With an error term, long, complex equations

• Uniform sampling, means $h$ is known

• These have simple error terms, making defining derivatives simple

Indefinite integral is class of functions $F$

Definite integral over fixed integral is number

F can be hard to get symbolically

Lower & Upper sums

Lower sum based on lower end of rectangles

• Uses largest lower bound, smallest width

Upper sum based on upper end of rectangles

• Uses smallest upper bound, small widths

Are both incorrect estimates

Goal is for something to be Reimann-integrable

Richardson's Extrapolation gives a formula for approximating M for any $h$

• Even powers of $h$ give better extrapolation results

Estimating polynomial written in lagrange form, using 3 points (see above)

Finds best-fit quadratic polynomial

Composite form of rule:

• Requires even integrals

• $S(f; P) = \frac{h}{3}(f(a) + f(b) 4\sum_{{k = 1}}^{{n/2}} f(a + (2k - 1)h) + 2\sum_{k = 1}^{(n - 2)/2} f^{{(4)}}(a + 2kh))$

• Error: $S - \frac{1}{180}(b - a)h^{4}f^{(4)}(\zeta)$

3/8 rule, from cubic spline interpolation, 4 point closed rule (slide 5)

• $(3/8)h(f(x0) + 3f(x1) + 3f(x2) + f(x3)) - (3/80)h^{5} f^{(4)}(\zeta)$

Boole's Rule

• Much more complex, $h^7$ error

• Similarly, Newton-Cotes rule has about the same error

Romberg's Integration

Given a new $T$ function, trapezoid rule across the whole

Random sampling

Monte Carlo, Monaco is the source of the name

Computational methods require pseudo-random number generators

Monte Carlo method for finding value of Pi

• shoot a thousand darts at a board with a circle, of unit size, and divide number in over total, and multiply by 4

Can be used easily for integration

Current provided by battery

• Can be solved with simple set of linear equations…

• Can be solved using, for instance, Gaussian elimination

Can be over or under-determined

Can also be ill conditioned

Only going to consider linear equations

Na\"ive Gaussian Elimination is possible, but potentially difficult

$Ax = B$ also gives $x = A^{{-1}} b$

• Inverse must be square, and $det(A) \neq 0$ ($A$ is full ranked)

• Least Squares method gives $\hat{x} = (A^{T} A)^{{-1}} A^{T}b$

• Must not have redundant rows

Scaled Partial Pivoting

• Involves scales of $|a_{n1}|/|l_{n}|$

• Row with largest scale value used for first elimination

• Then from second column against all, with same maximum value

• Residual Vector $R = A\hat{x} - b$

Indicates how well-conditioned/ill-conditioned a matrix is

Ill conditioned matrices are harder to handle and give less precise results

Find chacateristic polynomial, and then from this eigenvalues.

Condition number $\kappa$ is $\frac{|\lambda_{{max}}|}{|\lambda_{{min}}|}$ (assuming matrix is normal)

• Normally, however, $\|A\|_{2} \cdot \|A^{-1}\|_{2}$

Norm is the $\sqrt{\rho(A^{T} A)}$

• $\rho$ function returns lambda max

Polynomial interpolation can have some problems

Piecewise linear approximations can have many points, difficult to evaluate

Can instead remove some points such that upper/lower bounds are satisfied

Fairly simple with time series and error threshold, start, min, max given.

• For loop finding max and min of slope.

• Add lines, given slopes, etc.

Hooke's Law

Used when there's a linear relationship

Hooke's law derived from least-squares approximation

Columns for $x$ and $y$

• Appears to be almost a line, but not quite

• Best fit approximating all points has several definitions, minimum occurs when partial derivatives are 0

• Square is used for error

• Linear Approximation is fairly simple

• Matlab calculation

Three Point Formulas

Upper/Lower bounds

Trapezoidal Rule

Simpson's Rule

Composite Trapezoidal, Simpson's

Error Estimation – intervals bounding the error

Systems of Linear Equations

Gaussian Elimination, Na\"ive & Partial Pivoting

Residual Vector

Condition Number

Include error terms for composite rules