# CS 370: Introduction to Scientific Computing https://www.student.cs.uwaterloo.ca/~cs370/ - also piazza group (Mark Prosser) # 2013-06-21 ## TODO: Missing notes # 2013-06-19 ## TODO: Transcribe hand-written notes # 2013-06-17 # Complex Numbers (Review) \[ \begin{align} z_1 &= cos(\theta_1) + i sin(\theta_1) &= e^{i \theta_1} \\ z_2 &= cos(\theta_2) + i sin(\theta_2) &= e^{i \theta_2} \\ z_1 z_2 &= cos(\theta_1 \theta_2) + i sin(\theta_1 \theta_2) &= e^{i(\theta_1 \theta_2)} \\ \end{align} \] Recall the unit circle on the complex plane. \[ \begin{align} e^{\pi i} &= -1 & e^{i \frac{\pi}{2}} &= i \\ e^{2\pi i} &= 1 & e^{-i \frac{\pi}{2}} &= -i \end{align} \] Given N: \[ \begin{align} \omega &= e^{i \frac{2 \pi}{N}} \\ \omega^{N} &= (e^{i \frac{2 \pi}{N}})^N &= e^{i 2 \pi} &= 1 & \Rightarrow N^{th} \text{complex roots of unity} \end{align} \] Given N observations: \[ \begin{align} 1, & \omega, & \omega^2, & ... & \omega^{N-1} & \text{repeats once in } N \text{ observations} \\ 1, & \omega^2, & \omega^4, & ... & \omega^{2(N-1)} & \text{repeats twice in } N \text{ observations} \\ \vdots \\ 1, & \omega^n, & \omega^2n, & ... & \omega^{n(N-1)} & \text{repeats n times in } N \text{ observations} \end{align} \] Going in the other direction: \[ \begin{align} 1, & \omega^{-1}, & \omega^{-2}, & ... & \omega^{-(N-1)} & \text{repeats once in } N \text{ observations} \\ 1, & \omega^{-2}, & \omega^{-4}, & ... & \omega^{-2(N-1)} & \text{repeats twice in } N \text{ observations} \\ \vdots \\ 1, & \omega^{-n}, & \omega^{-2n}, & ... & \omega^{-n(N-1)} & \text{repeats n times in } N \text{ observations} \end{align} \] # DFT Definitions Input: $f_0, f_1, ..., f_{N-1}$ Output: $f_0, f_1, ..., f_{N-1}$ \[ \begin{align} F_k &= \frac{1}{N}\sum_{n=0}^{N-1} f_n \omega^{-kn} & \omega &= e^{i\frac{2 \pi}{N}} \\ k &= 0, 1, ..., N-1 \\ F_k &= \frac{1}{N}\sum_{n=0}^{N-1} f_k \overline{\omega}^{kn} \\ F_0 &= \frac{1}{N}(f_0 + f_1 + ... + f_{N-1}) \\ F_1 &= \frac{1}{N}(f_0 + f_1 \omega^{-1} + f_2 \omega^{-2} + ... + f_{N-1} \omega^{-(N-1)} \end{align} \] ## Properties of $\omega = e^{i \frac{2 \pi}{N}}$ 1. $\omega ^ N = 1$ 2. $\omega ^ k = 1$ if $k = lN$ where $l$ is an integer 3. $1 + \omega + \omega^2 + ... + \omega^{N-1} = \left\{ \begin{array}{1 1} N & \text{ if } \omega = 1 \\ 0 & \text{ if } \omega \neq 1 \end{array}\right.$ Proof of 3: \[ \begin{align} & (1-\omega)(1 + \omega + \omega^2 + ... + \omega^{N-1}) &= 1-\omega^N = 0 \\ \text{if } \omega \neq 1 \text{, then} & (1 + \omega + \omega^2 + ... + \omega^{N-1}) &= 0 \\ \text{else } \omega = 1 \text{, then} & (1 + \omega + \omega^2 + ... + \omega^{N-1}) &= N \end{align} \] $\square$ Similarly, \[ \begin{align} 1 + \omega^2 + \omega^4 + ... + \omega^{2(N-1)} &= 0 \\ 1 + \omega^{(N-1)} + \omega^{2(N-1)} + ... + \omega^{(N-1)(N-1)} &= 0 \\ 1 + \omega^N + \omega^{2N} + ... + \omega^{(N-1)N} &= N \end{align} \] ## DFT \[ \left[ \begin{array}{c} F_0 \\ F1 \\ \vdots \\ F_{N-1} \end{array} \right] = \frac{1}{N} \left[ \begin{array}{cccc} \overline{\omega}^{0 \cdot 0} & \overline{\omega}^{0 \cdot 1} & \dots & \overline{\omega}^{0 \cdot (N-1)} \\ \overline{\omega}^{1 \cdot 0} & \overline{\omega}^{1 \cdot 1} & \dots & \overline{\omega}^{1 \cdot (N-1)} \\ \vdots & \vdots & & \vdots \\ \overline{\omega}^{(N-1) \cdot 0} & \overline{\omega}^{(N-1) \cdot 1} & \dots & \overline{\omega}^{(N-1) \cdot (N-1)} \end{array} \right] \left[ \begin{array}{c} f_0 \\ f1 \\ \vdots \\ f_{N-1} \end{array} \right] \] ## Example: Computing DFT \[ \left[\begin{array}{c}f_0\\f_1\\f_2\\f_3\end{array}\right] = \left[\begin{array}{c}1\\3\\3\\2\end{array}\right] \\ \left[\begin{array}{c}F_0\\F_1\\F_2\\F_3\end{array}\right] = \frac{1}{N} \left[\begin{array}{cccc} 1 & 1 & 1 & 1 \\ 1 & \overline{\omega}^1 & \overline{\omega}^2 & \overline{\omega}^3 \\ 1 & \overline{\omega}^2 & \overline{\omega}^4 & \overline{\omega}^6 \\ 1 & \overline{\omega}^3 & \overline{\omega}^6 & \overline{\omega}^9 \end{array}\right] \left[\begin{array}{c}1\\3\\3\\2\end{array}\right] \] \[ \begin{align} \omega &= e^{i\frac{2 \pi}{4}} = i \\ \overline{\omega} &= -i \\ \overline{\omega}^2 &= -1 \\ \overline{\omega}^3 &= i \\ \overline{\omega}^4 &= 1 \end{align} \] ### Remarks 1. Each column repeats a pattern. 2. First column sums to 1. The other columns sum up to 0. # 2013-06-14 Assignment 2: Please write down section # on the cover page. # Overview of Discrete Fourier Transform (DFT) ## Discrete Fourier Transform Given discrete samples $ f_0, f_1, ..., f_{N-1} $ $(f_0, f_1, ..., f_{N-1}) \overrightarrow{DFT} (F_0, F_1, F_{N-1})$ DFT has the following properties: 1. transformation is reversible 2. transformation can be done efficiently (FFT) 3. $ F_0, F_1, ... F_{N-1} $ contains info about how fast a pattern repeats in the original data. (This is hard to understand.) ## Wide Applications - data compression (MP3, JPEG) - signal processing - image processing/comparison > ## Example 1 > > spot prices, e.g. orange price, electricity price > > objective: determine trends (repetition patterns in the price) > > data: $ 20 $ years of monthly spot price of oranges > > 1. seasonal patterns: high in spring, low in fall > > 2. el nino (every $ 7 $ years) > > 3. sunspot activity > > 4. business cycles > ## Example 2 > sound recording which is a superimposition of the sound of birds chirping and train whistling ## Four Parts 1. Definition of DFT 2. Reversible (IDFT) 3. How to compute fast: FFT 4. What does it mean? # Introduction to Discrete Fourier Transform Assume data comes from sampling. Represent unknown function in a form which helps understanding how fast repeat. In general: $cos t, sin t$ repeats once in $[0, 2\pi)$, ..., $cos kt, sin kt$ repeat $ k $ times in $ [0, 2\pi) $. If $ f(t) \approx C_0, a_1 cos(t) + b_1 sin(t) + a_2 cos(2t) + b_2 sin(2t) + ... + a_k cos(kt) + b_k sin(kt) $ , then: $ a_1^2 + b_1^2 $ measure strength of repeating once in $ [0, 2 \pi) $ , ... , $ a_k^2 + b_k^2 $ measure strength of repeating k times in $ [0, 2 \pi) $ # Discrete Case Input: $ N $ observations $ f_0, f_1, ... f_N $ Want to know: How many times a pattern is repeated in $N$ observations. # Review: Complex Numbers $z = a + bi$ $a$ and $b$ are real numbers. $i = \sqrt{-1}$ : imaginary unit $\overline{z} = a - b i$ $|z| = \sqrt{a^2 + b^2}$ $e^{i \theta} = cos \theta + i sin \theta$ (unit circle) # TODO: Transcribe hand-written notes # 2013-06-03 # Improved Euler \[ \begin{cases} \frac{dy}{dt} &= f(t, y) \\ y(t_0) &= y_0 \end{cases} \] \[ \begin{align} Euler: & y_{n+1} = y_n + h \cdot f(t_n, y_n) & Explicit \\ Trapezoidal: & y_{n+1} = y_n + h \cdot f(t_{n+1}, y_{n+1}) & Implicit \end {align} \] To make Trapezoidal explicit. 1. Use forward Euler as the first approximation, denote it as $y_{n+1}^{*}$ 2. Use Trapezoidal formula but replace $f(t_{n+1}, y_{n+1})$ by $f(t_{n+1}, y_{n+1}^{*})$ > # Improved Euler > > \[ \begin{align} y_{n+1}^{*} &= y_n + h \cdot f(t_n, y_n) & Explicit \\ y_{n+1} &= y_n + h \cdot f(t_{n+1}, y_{n+1}^{*}) \end {align} \] # Local Truncation Error (LTE) Analysis LTE: One Step Accuracy How well the formula approximates Taylor expansion for y(t_{n+1}). 1. Write method as an equation $ y_{n+1} - y_n - \frac{h}{2}(f(t_n, y_n) + f(t_{n+1}, y_{n+1})) = 0 $ 2. Replace y_n # 2013-05-31 # Recap ## Standard Form IVP \[ \begin{cases} y'(t) &= f(t, y(t)) \\ y(t_0) &= y_0 \end{cases} \] Determine: \[ \begin{align} t_0 < t_1 < ... < t_N \\ n_0 < n_1 < ... < n_n \end{align} \] $y_n \approx y(t_n), n=1,2,...,N$ ## Forward Euler Given $\alpha, t_0 < t_1 < ... < t_N, f(t, y),$ let $h_n = t_{n+1} + t_n, y_0 = \alpha$ 1. for $n = 0, 1, ..., N-1$ a. $y_{n+1} = y_n + h_n f(t_n, y_n)$ 1. end **How good is this method?** # Local Truncation Error (LTE) Main Tool: Taylor Expansion $y(t) = y(t_n) + y'(t_n)(t-t_n) + \frac{y''(t_n)}{2!}(t-t_n)^2 + ... + \frac{y^{(p)}(t_n)}{p!}(t-t_n)^p + Rest$ $Rest = \frac{y^{(p+1)}(\eta)}{(p+1)!}(t-t_n)^{(p+1)}$, $\eta$ is between $t_n$ and $t$ Assume that $\left|\frac{y^{(p+1)}(\eta)}{(p+1)!}\right| \leq C$ Therefore, $|Rest| \leq C|t-tn|^{(p+1)}, t->t_n$ Therefore $Rest = O((t-t_n)^{p_n})$ Forward Euler: $y_{n+1} = y_n + h_n f(t_n, y_n)$ Assume at %t_n, y_n = y(t_n), t_{n+1} - t_n = h% \[ \begin{align} y(t_{n+1}) &= y(t_n) + y'(tn) \cdot h + O(h^2) \\ &= y_n + f(t_n, y_n) \cdot h + O(h^2) \\ &= y_{n+1} + O(h^2) \end{align} \] ## Trapezoid Method Following Taylor \[ y(t_{n+1}) = y(t_n) + y'(t_n) \cdot h + \frac{y''(t_n)}{2}h^2 + O(h^3)$ [1] \] Apply Taylor to $y'(t)$ \[ \begin{align} y'(t_{n+1}) &= y'(t_n) + y''(t_n) \cdot h + O(h^2) & [2] \\ y''(t_n) &= \frac{y(t_{n+1}) - y'(t_n)}{h} + O(h) \\ \end{align} \] Substitute [2] into [1] y(t_{n+1}) = y(t_n) + y'(t_n) \cdot h + \frac{1}{2} \cdot \left[\frac{y(t_{n+1}) - y'(t_n)}{h} + O(h)\right\] \cdot h^2 + O(h^3) just a little bit of manipulation y(t_{n+1}) = y(t_n) + \frac{h}{2} \cdot \left[f(t_n, y(t_n)) + f(t_{n+1}, y(t_{n+1})) + O(h^3) This suggests: y_{n+1} = y_n + \frac{f(t_n, y_n) + f(t_{n+1}, y_{n+1})}{2} > ## Trapezoidal Method > Given \alpha, t_n = t_0 + n \cdot h, f(t, y) \\ > y_0 = \alpha \\ > for n=0, 1, ..., N-1 \\ > solve: y_{n+1} = y_n + > \frac{h}{2}(f(t_n, y_n)+f(t_{n+1}, y_{n+1})) > end > > Note that y_{n+1} is on both sides of the equation. > This method is more accurate, but requires solving an > equality ("implicit method", <=> > forward euler / "explicit method") > ## Example > \[ > \begin{cases} > y'(t) &= t - y(t)^2 \\ > y(0) &= 2 > \end{cases} \\ > t_N = 10, N= 50, h = 0.2 \\ > y_0 = 2 > ... blah blah blah > \] # 2013-05-29 Assignment 1: data for handwriting can be saved in a .mat file and submitted to the dropbox. # System of 1st order ODE IVP \frac{dy}[dt} = f(t,y) y(t_0) = y_0 f(t,y): dynamics function > ## Example > \frac{dy}{dy} = -10 y > y(0) = 1 > > f(t,y) = -10 y > > y=100 > t=10 > > => f(10, 100) = -1000 = \frac{dy}{dt} \frac{dy_1}{dt} = f(t, y_1, ..., y_n) ... \frac{dy_n}{dt} = f(t, y_1, ..., y_n) y_1(t_0) = \alpha_1 ... y_n(t_0) = \alpha_n y = [y_1(t), ..., y_n(t)] f = [f_1(t, y_1, ..., y_n), ..., f_n(t_, y_1, ..., y_n)] y_0 = [\alpha_1, ..., \alpha_n] # Higher Order ODE IVP An n^{th} order IVP. \frac{d^ny}{dt^n} = g(t, y, \frac{dy}{dt}, ... \frac{d^{n-1}y}{dt^{n-1}}) y(t_0)=\alpha_1, \frac{dy}{dt}(t_0) = \alpha_2, ... \frac{d^{n-1}y}{dt^{n-1}}(t_0) = \alpha_n > ## Example > z_1 = y > z_2 = \frac{dy}{dt} > ... > z_n = \frac{d^{n-1}y}{dt^{n-1}} > > \frac{dz_1}{dt} = \frac{dy}{dt} = z_2 > \frac{dz_2}{dt} = \frac{d}{dt} \frac{dy}{dt} = \frac{d^2y}{dt^2} = z_3 > \frac{dz_{n-1}}{dt} = ... = z_n > \frac{dz_n}{dt} = ... = \frac{d^ny}{dt^n} = g(t, z_1, ..., z_{n-1}) > NOTE: We need Taylor expansion for the next week or two. # Euler's Method (Forward Euler) Goal: Compute y(t_0), y(t_1), ... , y(t_n) for t_0 < t_1 < ... < t_n. We want to obtain y_0, y_1, ..., y_n with y_k \approx y(t_k), k=0, ..., n Let h_{k-1} = t_k - t_{k-1} (step size) ## Idea At t_{k-1}, y_{k-1} \approx y(t_{k-1}). Use y_{k-1} and f(t_{k-1}, y_{k-1}) to compute y_k. At t = t_0, y(t_0) = y_0 is known. \frac{dy}{dt}(t_0) = f(t_0, y(t_0)) => slope of y(t) at t_0 is known. Let \frac{y_1-y_0}{t_1-t_0} = f(t_0, y_0) => y_1 = y_0 + h_1 f(t_0, y_0) In general, assume y_{k-1} = y(t_{k-1}). Forward Euler: y_k = y_{k-1} + h_{k-1} \cdot f(t_{k-1}, y_{k-1}). y_0 is given. > ## Algorithm > > y_0 = \alpha > for k = 1:N > y_k = y_{k-1} + h_{k-1} \cdot f(t_{k-1} > end > ## Example > y'(t) = t - y(t)^2 > y(0) = 2 > > t_0 = 0, t_n = 10, n=50, h=\frac{10}{50} = 0.2 > f(t, y) = t - y^2 > y_0 = 2 > y_1 = 2 + 0.2(0 - 2^2) = 2 - 0.2(4) = 1.2 # 2013-05-27 # Ordinary Differential Equations - <=> Partial Differential Equations - We will not be dealing with partial differential eqns in this course. ## Problem Description Given a function f, system dynamics function, and initial value y, solve: y'(t) = f(t, y(t)), t_0 \lte t \lte t_{final} y(t_0) = y_0 - y(t), y_0 and f can be *scalars* or *vectors*. - We know the location of an object at t_0, y(t_0) = y_0 - We known how it's going to move at any time t; at time t, location is y(t); speed/rate of change is f(t, y(t)) - Goal: Determine when the object is at *any* time t; i.e. determine y(t). > ## Example > Bank account has balance y(t_0) at t=t_0. > Assume a continuously compounding interest rate, r. > \frac{dy}{dt} = y \dot r > y(t_0} = y_0 > Answer: y(t) = y_0 * e^{r(t-t_0)} > ## Example > Population growth model. > r = birth rate - death rate > ## Example > Multiple objects, population of species. > x'(t) = a x(t) - \alpha x(t) y(t) > y'(t) = -b y(t) + \Beta x(t) y(t) > a, \alpha, b, \Beta constants > x(0) = x_0 > y(0) = y_0 > How do we express this? > Let z(t) = [z_1(t), z_2(t)] = [x(t), y(t)] > z(0) = [z_1(0), z_2(0)] = [x_0, y_0] > f(t, z_t) = [a z_1(t) - \alpha z_1(t) z_2(t), > -b z_2(t) + \Beta z_1(t) z_2(t)] > z'(t) = f(t, z_t). > This does not have an explicit (closed form) solution. > Compute solution *numerically*. ## Initial Values are Important Solution to y'(t) = f(t, y(t)) is not unique. > ## Example > \frac{dy}{dt} = t-y > f(t,y) = t-y > For any c, y(t) = c e^{-t} + t - 1 satisfies > \frac{dy}{dt} = t-y . > With an initial condition y(0) = 1, we get c=2, > therefore y(t) = 2 e^{-t} + t - 1 . ## Higher Order ODEs > ## Example > Pendulum movment > angle theta, arm length L, force of gravity m*g > > \frac{d^2\theta}{dt^2} + \frac{c}{mL} \frac{d\theta}{dt} + \frac{g}{L}sin(\theta) = 0 (**) > > damping force: -c \frac{d\theta}{dt} > > \theta(0) = \theta_0 : initial position > \frac{d\theta}{dt}(0) = v_0 : initial velocity > > This is a *second order* ODE. > \frac{d^2u}{dt^2} = f(t,u,\frac{du}{dt}) > u(t_0) = u_0 > \frac{du}{dt}(t_0) = v_0 > > We can transform this into a system of first order ODEs > Consider (**). Let > y_1 = \theta, > y_2 = \frac{d\theta}{dt} > Rewrite (**) in terms of y_1 and y_2 only: > Note \frac{d y_1}{dt} = \frac{d \theta}{dt} = y_2 > \frac{d y_2}{dt} = \frac{d^2\theta}{dt^2} = - \frac{c}{mL}y_2 - \frac{g}{mL}sin(y_1) > ... etc. # 2013-05-24 > *What is a "piecewise cubic spline"?* > A polynomial of degree d \leq 3, that goes through a > given set of points, and has continuous value, first > derivative, and second derivative. # Compute Cubic Spline Interpolant s(x) Interpolating (x_1, y_1) ... (x_n, y_n) => 2(n-1) eqn. Continuity of s'(x) and s''(x) => 2(n-2) eqn. Boundary conditions => 2 eqn. s(x) has 4(n-1) unknowns. Naive solution: solve the system of 4(n-1) equations. There's a better way! Solve n equations in a special structure. Idea today, actual computation will be in the final lecture (!) for solving systems of equations. # More efficient computation of s(x). - Consider Hermite piecewise cubic interpolation. (x_1, y_1, s_1), ..., (x_1, y_1, s_n) => s(x) I don't know the slope, but I will choose the slopes s_i so that s''(x) is continuous and the boundary conditions are satisfied. => Solve of linear system of equations in n unknowns. > csape(x, y, cmd) > (MATLAB) computes a cubic spline interpolant _Once you determine s(x), you usually want to compute values of s(x) so that you can e.g. plot s(x), or interpolate a a value._ To compute spline: x_ref = linspace(0, 10, 100) > ppval(s, xref) > (MATLAB) evaluate a _Piecewise Polynomial_ y_ref = ppval(s, x_ref) x_1 < x_2 < ... < x_n => y_i = f(x_i) for some unknown f ## Parametric Curves in 2D Sometimes, our data points (x_i, y_i) cannot be represented in the form y = f(x). (Example: unit circle.) We can represent such 2D data by a pair of continuous functions x(t), y(t) with share a common parameter (t). > Example: unit circle > x(t) = cos(\pi t), t \elt [0, 2] > y(t) = sin(\pi t), t \elt [0, 2] ## Computing parametric curves Given (x_1, y_1), ..., (x_n, y_n), we want to compute (x(t), y(t)) to represent this data smoothly (using splines). 1. _Keep track of points by an auxillary parameter t._ Represent points as (t_1, x_1), ..., (t_n, x_n) and (t_1, y_1), ..., (t_n, y_n). 2. _Compute two interpolating splines using csape._ s_x(t) : interpolates (t_1, x_1), ..., (t_n, x_n) s_y(t) : interpolates (t_1, y_1), ..., (t_n, y_n) Question: So... how do we choose t? 1. Position/order of points as t. (For (x_i, y_i), let t_i = i.) Problem: This order does not capture the density of points. 2. Use arc length as t. t_1 = 0 t_2 = t_1 + distance (x_1, y_1) to (x_2, y_2) ... t_n = t_{n+1} + distance (x_{n-1}, y_{n-1}) to (x_n, y_n) 3. To draw s_x and s_y, determine t_ref (see course notes to see how to generate t_ref). For handwriting drawing question, you need to break up letters into multiple segments and generate a smooth representation for each segment. Note that cursive letters actually have discontinuities, you will need to decompose the letters into splines split where the derivative is not continuous. # 2013-05-22 ## Review: Piecewise (low degree) polynomial Given x_1 < x_2 < ... < x_n p(x) = p_1(x), x \elt [x_1, x_2] = p_2(x), x \elt [x_1, x_2] ... p_{n-1}(x), x \elt [x_{n-1}, x_{n}] Interpolate (x_1, y_1), ..., (x_n, y_n). => p(x_i) = y_i, i.e. p_i(x_{i+1}) = y_{i+1} = p_{i+1}(x_{i+1}), for i = 1, ..., n-2 Problem? Derivative may not be defined at (x_i, y_i). ## Piece-wise Hermite Interpolation Given (x_i, y_i, s_i) for i = 1..n, x_1 < ... < x_n There exists a unique piecewise cubic polynomial: p(x) = p_1(x), x \elt [x_1, x_2] = p_2(x), x \elt [x_2, x_3] ... = p_{n-1}(x), x \elt [x_{n-1}, x_n] p_i(x) is cubic p_i(x) = a_i + b_i(x-x_i) + c_i(x - x_i)^2 + d_i(x-x_i)^3 such that p(x_i) = y_i, p'(x_i) = s_i, i = 1, ..., n In addition, a_i = y_i, b_i = s_i, c_i = \frac{3 y'_i - 2 s_i - s_{i+1}}{\delta x_i} d_i = \frac{s_{i+1} + s_i + 2 y'_i}{(\delta x_i}^2} y'_i = \frac{y_{i+1} - y_i}{\delta x_i} p(x) and p'(x) are continuous at x_2, ..., x_{n-1}. ## Cubic Splines Given x_1 < x_2 < ... < x_n, a cubic spline s(x) is a piecewise cubic polynomial with s(x), s'(x), s''(x) are continuous in [x_1, x_n]. A cubic spline has the additional property s(x_i) = y_i, i=1,...,n. Can it be done? s(x) = s_1(x), x \elt [x_1, x_2] ... s(x) = s_{n-1}(x), x \elt [x_{n-1}, x_n] Let s_i(x) = a_i + b_i(x - x_i) + c_i(x - x_i)^2 + d_i(x - x_i)^3. Defining s(x) has 4(n-1) unknowns (four unknowns for each s_i(x), i = 1, ..., n_1). s(x) interpolates data, so: s_1(x_1) = y_1, s_{n_1}(x_n) = y_n. s(x) must be continuous, so: s_i(x_{i+1}) = y_{i+1} s_{i+1}(x_{i+1}) = y_{i+1} # of equations: 2 + 2(n-2) = 2(n-1). Still need more equations to make system of equations for s(x) defined. s'(x) must be continuous at x = x_2, ... x_n-1, so: s'_i(x_{i+1}) = s'_{i+1}(x_{i+1}), i=1,...,n-2 s''(x) must be continuous, so: s''_i(x_{i+1}) = s''_{i+1}(x_{i+1}), i=1,...,n-2. This yields another 2(n-2) equations, for a total of 2(n-1) + 2(n-2) equations for a total of 4n - 6 equations. We need 4n - 8 equations to make s(x) defined. Therefore, we need two more conditions to define s(x). ## Boundary Conditions at x_i and x_n ### Natural Spline s''_1(x_1) = 0, s''_{n-1}(x_n) = 0 This results in no curvature at x_1 and x_n. ### Clamped Spline s'_1(x_1) = s_1 s'_{n-1}(x_n) = s_n ### Not-a-Knot Spline s'''(x) at x_2 and x_n-1 are continuous. Look at x_1, x_2, x_3; x_{n-2}, x_{n-1}, x_n. ## Computation of Spline ### Example Can we find a_1, a_2, ..., a_5 such that s(x) is a natural spline? s(x) = 20 + 25x + 18x^2 + a_1 x^3, x \elt [-3, -1] s(x) = 26 + a_2 x + a_3 x^2 + a_4 x^3, x \elt [-1, 0] s(x) = 26 + 19x + a_5 x^2 + x^3, x \elt [0, 3] ### Solution First, compute piecewise s'(x), s''(x). s'_1(x) = 25 + 36x + 3 a_1 x^2 s'_2(x) = a_2 + 2 a_3 x + 3 a_4 x^2 s'_3(x) = 19 + 2 a_5 x + 3x^2 s''_1(x) = 36 + 6 a_1 x s''_2(x) = 2 a_3 + 6 a_4 x s''_3(x) = 2 a_5 + 6x Tip: Look at higher order derivatives first -- they're usually linear, and thus are simpler. s''(-3) = 0 ? s''(3) = 0? s''_1(-1) = s''_2(-1)? s''_2(0) = s''_3(0)? s'_1(-1) = s'_2(-1)? s'_2(0) = s'_3(0)? s_1(-1) = s_2(-1)? s_2(0) = s_3(0)? x_1 = -3, x_2 = -1, x_3 = 0, x_4 = 3 # 2013-05-17 Mg 01-1, next Tues 2-3pm ## Hermic Polynomial Interpolation Given (x_1, y_1), (x_2, y_2), ..., (x_n, y_n) with x_1, x_2, ..., x_n distinct, a unique polynomial p(x) of degree not exceeding n-1 can be determined to satisfy p(x_i) = y_i, i=1, 2, .. n - monomial representation - lagrange representation We can interpolate derivatives as well ### Example #### Problem Suppose we have x_L, y_L, s_L, x_R, y_R, s_R. Find a cubic polynomial p(x) such that p(x_L)=y_L p(x_R)=y_R p'(x_L)=s_L p'(x_R)=s_R #### Solve p(x) = a + b(x-x_L) + c(x-x_L)^2 + d(x-x_L)^3 > p(x_L) = y_L, so a = y_L p'(x) = b + 2c(x-x_L) + 3d(x-x_L)^2 > p'(x_L) = s_L, so b = s_L > p(x_R) = y_R, so a + b(x_R-x_L) + c(x_R-x_L)^2 + d(x_R-x_L)^3 = y_R > p'(x_R) = s_R, so b + 2c(x-x_L) + 3d(x-x_L)^2 = s_R #### Solution a = y_L b = s_L c = \frac{3y_L' - 2s_L - s_R}{\delta x_L} d = \frac{s_L+s_R-2y_L'}{(\delta x_L)^2} where \delta x_L = x_R - x_L y_L' = \frac{y_R - y_L}{x_R - x_L} Hermite interpolation: interpolating + derivatives. ## Problems with high degree polynomial interpolation For large n, polynomial interpolation is problematic. Range Function f(x) = \frac{1}{1+25x^2} > - High degree polynomial tends to wiggle. > - Small change to data can lead to large change in solution. - conclusion: avoid high degree polynomials - in practice, polynomial of degree \gte 5 is rarely used! ## Piecewise Polynomial Interpolation Given (x_1, y_1), ..., (x_n, y_n); x_1, x_2, ... x_n distinct. Assume x_1 < ... x_n. Find: 1. p(x) interpolates data > p(x_i) = y_i for x=1,...,n 2. p(x) is a polynomial in each [x_i, x_{i+1}] p(x) = p_1(x), x \elt [x_1, x_2] = p_2(x), x \elt [x_2, x_3] ... p(x) = p_{n-1}(x), x \elt [x_{n-1}, x_n] 3. p(x) has desired smoothness in [x_i, x_n] ## Piecewise Linear Polynomial Interpolation - p(x) connects the data points by straight lines - p(x) is continuous (because of the (x_i, y_i) points connecting the p_i(x)) - p'(x) at x_i, i=2, ..., n-1 is not generally continuous ## Piecewise Hermite Interpolation p(x) = p_1(x), x \elt [x_1, x_2] ... p(x) = p_{n-1}(x), x \elt [x_{n-1}, x_n] Find p(x) such that p(x_i) = y_i, i=1, ..., n, and p'(x_i) = s_i, i=1, ..., n - since p'(x) is defined for all x_i, i=1,...,n - p(x) is continuous - p'(x) is continuous at x_i, i=2...n Next week: Splines! # 2013-05-15 ## Polynomial Interpolation (Section 2) Given n points (x_1, y_1), (x_2, y_2)... (x_n, y_n) with x_1, x_2, ... x_n distinct Find a polynomial p(x) of degree \leq n-1 satisfying p(x_1) = y_1, ... P(x_n) = y_n p(x) \approx f(x) p(x) interpolates f(x) as x_i, i=1,2,...n ### Example #### n = 1 (x_1, y_1), p(x) = y_1 #### n = 2 (x_1, y_1), (x_2, y_2) p(x) = y_1 + \frac{y_2 - y_1}{x_2 - x_1}(x - x_i) What about p(x) for arbitrary n? ## Monomial Representation p(x) = C_1 + C_2 x + C_3 x^2 + ... + C_n x^{n-1} p(x_1) = y_1 C_1 + C_2 x_1 + C_3 x_1^1 + ... + C_n x_1^{n-1} = y_1 p(x_2) = y_2 C_1 + C_2 x_2 + C_3 x_2^1 + ... + C_n x_2^{n-1} = y_2 ... p(x_n) = y_n C_1 + C_2 x_n + C_3 x_n^1 + ... + C_n x_n^{n-1} = y_2 (Oh, hey, a system of n equations in n unknowns!) V * \vec{C} = \vec{y} \[ \left[ \begin{array}{5}{ccccc} 1 & x_1 & x_1^2 & ... & x_1^{n-1} \\ 1 & x_2 & x_2^2 & ... & x_2^{n-1} \\ & & ... & & \\ 1 & x_n & x_n^2 & ... & x_n^{n-1} \end{array} \right] \] \[ \left[ \begin{array}{1}{c} c_1 \\ c_2 \\ ... \\ c_n \end{array} \right] ] = \[ \left[ \begin{array}{1}{c} y_1 \\ y_2 \\ ... \\ y_n \end{array} \right] ] This has a unique solution when the matrix is invertable (not singular). > ## V: Vandermonde matrix > FACT: > V is nonsingular when x_1, x_2, ..., x_n are distinct. > > det(V) = \Pi_{1 \leq i < j \leq n} (x_j - x_i) ### Theorem Given n points (x_1, y_1), ..., (x_n, y_n) with x_1, x_2, ..., x_n distinct, there exists a unique polynomial p(x) of degree not exceeding n-1 such that p(x_i) = y_i, i=1, 2, ..., n ### Question If degree of P(x) is allowed to be > n-1, how many p(x) interpolations are there for n data points? Answer: infinitely many ### Question If the degree of p(x) < n-1, how many? Answer: Either 0 or 1. ## Lagrange Form Given x_1 = 1, x_2 = 3, x_3 = 4 Consider L_1(x) = \frac{(x-3)(x-4)}{(1-3)\dot(1-4)} = \frac{(x-3)(x-4)}{6} L_2(x) = \frac{(x-1)(x-4)}{(3-1)\dot(3-4)} = \frac{(x-1)(x-4)}{(-2)} L_3(x) = \frac{(x-1)(x-3)}{(4-1)\dot(4-3)} = \frac{(x-1)(x-3)}{3} Hmm. L_1(x_1) = 1, L_1(x_2) = 0, L_1(x_3) = 0 L_2(x_1) = 0, L_2(x_2) = 1, L_2(x_3) = 0 L_3(x_1) = 0, L_3(x_2) = 0, L_3(x_3) = 1 Now find the interpolating polynomial for: (1, 17), (3, -2), (4, 11) p(x) = 17 \dot L_1(x) - 2 \dot L_2(x) + 11 \dot L_3(x) Because of the earlier theorem, we note that this p(x) has to be the same as the one in monomial form! For a general x_1, x_2, ..., x_n distinct, we have L_1(x), L_2(x), ..., L_n(x) of degree n-1 L_i(x) = \frac{(x-x_1)(x-x_2)...(x-x_{i-1})(x-x_{i+1})...(x-x_n)}{(x_i-x_1)...(x_i-x_{i-1})(x_i-x_{i+1})...(x_i-x_n)} L_i(x_j) = 1 if j=1, 0 if j \neq i Determine polynomial interpolant function (x_1, y_1)(x_2, y_2)...(x_n)(y_n) P(x) = y_1 \dot L_1(x) + ... + y_n \dot L_n(x) L_i(x) <- lagrange bases ## Closing Remarks Vandermonde and Lagrange yield the same polynomial. How to verify? - evaluate both at n different points (degree of n-1) - rewrite lagrange in standard form and compare c_i's to vandermonde # 2013-05-13 ## Review FPNS: F(\Beta, t, L, U) fl(x) + (1 + \delta)x, |\delta| < \Epsilon | \frac{x (+) y - (x + y)}{x+y} | \leq (2E + E^2) * \frac{|x| + |y|}{|x+y|} Watch out for catastrophic error when the addition is really subtraction! E = \Beta ^ {1-t} if truncation = \frac{\Beta^{1-t}}{2} if rounding Machine epsilon E is the smallest x such that f(1+x) > 1. NOTE: E does NOT have to be representable in the FPNS. (Not in the notes!) ## Roundoff Error Compute I_n = \int_0^1 \! \frac{x^n}{x+\alpha} \mathrm{d}x \alpha: a given parameter n: integer I_0 = \int_0^1 \! \frac{1}{x+\alpha} \mathrm{d}x = log(x + \alpha) \bigg|_0^1 I_0 = log(1+x)-log(\alpha)-log(\frac{1+\alpha}{\alpha}) I_{n+1} = \int_0^1 \! \frac{x^{n+1}}{x+\alpha} \mathrm{d}x = \int_0^1 \! \frac{x^n(x+\alpha-\alpha)}{x+\alpha} \mathrm{d}x = \int_0^1 \! x^n \mathrm{d}x - \alpha \int_0^1 \! \frac{x^n}{x+\alpha} \mathrm{d}x = \frac{x^{n+1}}{n+1}\bigg|_0^1 - \alpha I_n = \frac{1}{n+1} - \alpha I_n I_{n+1} = \frac{1}{n+1} - \alpha I_n ## Example ### Algorithm I_0 = log(\frac{1+\alpha}{\alpha}) for K = 1:n I_K = \frac{1}{K} - \alpha * I_{K-1} end ### Results Using a single precision computer, \alpha = 0.5, I_100 = 0.00664 \alpha = 2, I_100 = 2.1*10^22 (Hmm...) But \alpha \geq 0 \int_0^1 \frac{x^n}{x+\alpha} \mathrm{d}x = \int_0^1 {x^{n-1}}\frac{x}{x+\alpha} \mathrm{d}x \leq \int_0^1 x^{n-1} \mathrm{d}x = \frac{x^n}{n} \bigg|_0^1 = \frac{1}{n} ## Stability Analysis of Algorithm Assume error only occurs at n=0. What is the effect of this initial error at the end of the computation? Consider: (1) I_0 = Exact. I_n = \frac{1}{n} - \alphaI_{n-1} (2) \hat{I_0} = Approx. \hat{I_n} = \frac{1}{n} - \alpha\hat{I_{n-1}} Let e_n = I_n - \hat{I_n} (1) - (2) I_n - \hat{I_n} = -\alpha(I_{n-1} - \hat{I_{n-1}}) e_n = -alpha e_{n-1} = (-\alpha)^2 e_{n-2} ... = (-\alpha)^n e_0 |\alpha| < 1, e_n \rightarrow 0 as n \rightarrow +\inf |\alpha| > 1, e_n = |\alpha|^n|e_0| \rightarrow +\inf as n \rightarrow +\inf Assume |\alpha| > 1. Is there a stable algorithm to compute I_n? > ## Stability An algorithm is stable if error does not become exponentially amplified. ## Example Compute the intersection of two straight lines. Seems straightforward enough. Where might there be catastrophic error? If the lines are almost parallel. The problem is very sensitive to a little bit of perturbation. > "If you get a bad answer, is it because your problem is not formulated well, or is it because your algorithm is bad?" > Accuracy of computed solution also depends on the _condition_ of the problem. - stability <=> algorithm - condition <=> problem instance > ## Condition > A problem is well conditioned if the change in solution is small when the (input data of the) problem is changed slightly. ## Interpolation (Section 2) Input: (x_1, y_1), (x_2, y_2), ..., (x_n, y_n) Output: A smooth representation of data Assume: x_1 < x_2 < ... < x_n Imagine y_k = f(x_k) for some unknown function Task: Find a smooth P(x) such that P(x_i) = y_i for i = 1..n P(x) \approx f(x) # 2013-05-10 ## Review FPNS: F(\Beta, t, L, U) = - \Beta = base - t = number of digits - L = lower bound on exponent - U = upper bound on exponent fl(x) = x(1+\delta) |\delta| \leq E E = \Beta^{1-t} if truncation = \frac{1}{2} \Beta^{1-t} if rounding ## Machine epsilon / unit round-off E ### Example F(10, 3, -10, 10) x = \pi = 3.1415926... fl(\pi) = 3.14 \delta = \frac{0.314*10 - 0.31415936...*10}{3.1415926...} \delta = \frac{-0.000159...*10}{3.1415926...} = \frac{-0.159...*10^{-2}}{3.14159126...} |delta| < {1/2}*10^{-2} E is the smallest number in F such that 1+E > 1 1 = 0.10... * \Beta E = 0.0....01 * \Beta 1 + E = 0.10....01 * \Beta > 1 ## Arithmetic Operations e.g. addition Let x, y be real numbers. x (+) y = fl( fl(x) + fl(y) ) Example compute x (+) y using 4 binary digits for x=0.1111 and y=0.1110 fl(x) = x fl(y) = y x(+)y = fl( fl(x) + fl(y) ) = fl( 0.1111 + 0.1110 ) = fl( 1.1101 ) = 0.1110 * 2 ## FACT a(+)b = b(+)a ## FACT (a (+) b) (+) c \neq a (+) (b (+) c) ## FACT | \frac{x(+)y - (x+y)}{x+y} | \leq (2E + E^2)\frac{|x| + |y|}{|x + y|} - avoid/be aware of the instance where your result is a small number determined by subtracting two large numbers ## Catastrophic Cancellation When X and Y have different signs, and x+y \approx 0 => large relative error in fl(x+y) ### Example Using four-digit arithmetic to compute the roots of x^2+62.1x+1 = 0 x_1 = \frac{-b + \sqrt{b^2 - 4ac}}{2a} and x_2 = \frac{-b - \sqrt{b^2 - 4ac}}{2a} Answer: computed x_1 = -0.0200 computed x_2 = -62.1 ## FACT x_1 = \frac{-b + \sqrt{b^2 - 4ac}}{2a} = \frac{-2c}{b+\sqrt{b^2 - 4ac}} ### Example (cont'd) Using new x_1, computed x_1 = -0.0161 Correct answer: x_1 = -0.0161072374... x_2 = -62.08389... ### Example Use 5 digit arithmetic to compute e^{-5.5} We know: e^x = 1 + x + \frac{x^2}{2!} + \frac{x^3}{3!} + ... e^x = \frac{1}{e^-x} #### Direct application of taylor expansion e^{-5.5} = 1 - 5.5 + \frac{(5.5)^2}{2} - \frac{(5.5)^3}{3!} + ... + (25 terms) = 0.0026363 [computed answer 1] e^{-5.5} = \frac{1}{1 + 5.5 + \frac{5.5}{2!} + ... + (25 terms)} = 0.0040868 [computed answer 2] correct answer: 0.00408771 # 2013-05-08 - we need to understand the implication of error when doing computation on computers (which have a fixed precision) - Real numbers have infinite # of digits - Computer numbers have finite precision => Computer arithmetic is approximate. The errors can be large, causing major problems. ## Normalized Floating Point Numbers - Any real number can be represented as $ 0 $ or $+/- * d_1 d_2 ... d_t d_t+1 * \Beta^P$ where $ d_1 \neq 0 $ $\Beta$ : base, often $\Beta = 10, 2$ For example, in base 10, $\pi = 0.31415.... * 10^1$ - recall that, for any two non-equal real numbers, there is always at least one real number between them ## Floating Point Number System (FPNS) $F(\Beta, t, L, U) =$ ${ 0 } U { +/- * 0.d_1 d_2 ... d_t * \Beta^P, d_1 \neq 0, L \lte P \lte U }$ FPNS has a limited extent and limited density Example: $F(10, 2, -1, 2)$ - largest positive number is $0.99 * 10^2 = 99$ - smallest positive number is $0.10 * 10^{-1} = 0.0010$ - next smallest positive number is $0.11 * 10^{-1} = 0.0011$ - spacing between numbers is $0.01 * 10^{-1} = 0.0001$ | | How many | spacing -|-|-|- $p = -1$ | $[0.1 * 10^{-1}, 0.1 * 10^0)$ | $9 * 10 = 90$ | $0.01*10^{-1} = 0.001$ $p = 0$ | $[0.1 * 10^0, 0.1 * 10^1)$ | $9 * 10 = 90$ | $0.01*10^{0} = 0.01$ $p = 1$ | $[0.1 * 10^1, 0.1 * 10^2)$ | $9 * 10 = 90$ | $0.01*10^{1} = 0.1$ $p = 2$ | $[0.1 * 10^2, 0.1 * 10^3)$ | $9 * 10 = 90$ | $0.01*10^{2} = 1$ Total Numbers in $F$ $1 + 2 * 9 * 10 * 4 = 721$ $(zero, +/-, d_1, d_2, P)$ Single Precision: $F(2, 24, -126, 127)$ Double Precision: $F(2, 53, -1022, 1023)$ ## Floating Point Error $x$ : real number $fl(x) \elt F$ : floating point number representation Let $x > 0$ $x = d_1 d_2 ... d_t d_t+1 ... * \Beta^P, d_1 \neq 0, L \lte P \lte U$ truncation: $fl(x) = 0.d_1 d_2 ... d_t * \Beta^P$ rounding: $fl(x) = $ $d_1 ... d_t * \Beta ^ P if d_t+1 < \Beta/2$ $d_1 ... d_t * \Beta ^ P + \Beta^{P-t}$ otherwise Error: $| fl(x) - x |$ Relative error: $| fl(x) - x | / |x|$ Let $\delta = (fl(x) - x) / x$ note this is the same as $x * \delta = fl(x) - x$ or $\delta*x + x = fl(x)$ or $fl(x) = (1+\delta) * x$ can we bound $\delta$? Truncation: $x>0$ $\delta = 0.00000d_t+1 * \Beta^P / d_1...d_1.... * \Beta^P$ $= 0.d_t+1 * \Beta^-t / d_1.d_2... * \Beta^-1$ $|\delta| <= | 0.d_{t+1} ... / d_1.d_2... | * \Beta^(1-t)$ $<= \Beta^(1_t)$ <--- machine epsilon (~= upper bound on error in floating point) # 2013-05-06 ## Grades - 4 assignments (32%, 8% each) - midterm (28%) - June 26, DC 1350, 7pm-9pm - final (40%) - must pass final to pass course - cheating (-100%) ## Logistics - assignment/midterm solutions will be discussed in tutorials - open book exam (lecture notes, assignments available) - Programming will be done in MatLab - course notes at the copy centre (BUY) (DC? MC?) - ~$14 ? - optional texts available ## Why take this course? - big data - Google - signal processing - image processing - analysis & visualization in science and engineering applications ## What will be covered? - solving continuous mathematical problems using computers ## Review of Calculus - Taylor Expansion - approximate a function by using a polynomial and the derivatives of the function - Taylor expansion of $ f(x) $ at $ x_c $ - $ f(x_c + h) = f(x_c) + f'(x_c)*h + \frac{f''(x_c)*h^2}{2!} + ... + \frac{f^{(n)}(x_c)*h^n}{n!} $ - Why do we use the Taylor Expansion? - Computers have finite precision - Taylor expansions allow us to approximate arbitrary functions \w a computer ## Floating Point Arithmetic - computers have finite precision (!) - error can accumulate if you don't handle floating point correctly - Disasters due to floating point error - Gulf War, '91 - Patriot defense system - error in how time is tracked - 28 died, 100 injured because missile target missed - Vancouver Stock Exchange Index - started with a value of 1000 - 22 months of trading; value dropped to ~500 - after each transaction, the value of the index was truncated (!) ## Smooth representation of data - improves our ability to predict - leads to simpler models - graphics - spline interpolation - audio (interpolation) - A1 (?) ## (Ordinary Differential Eqns) - $\frac{dy}{dt} = f(t,y); y(t_0) = y_0$ - simulation (angry birds, olympic jumpers, etc.) - A2 ## Fourier Transform - "not an easy topic" - we observe signals (e.g. audio) - try to understand the repetition patterns observed in the signal - $f(t) => \hat{f}(xi)$ (time domain to frequency domain) - also will cover FFT - A3 (usually image compression) ## Solving Linear Systems of Equations Ax=b - "fast", "accurate" - e.g. pagerank - A4