ODEs, Polynomials approx., Linear Least Squares, and Errors

Machine epsilon

$fl(x)=x(1+ϵ)\text{where }|ϵ|\le u$ $|\frac{fl(x)-x}{x}|=|ϵ|\le u\text{is called relative error.}$ $\text{Cancellations occur when subtracting nearby number containing roundoff.}$

Taylor series

$$\begin{array}{rl}f(x)& =\sum _{k=0}^{\inf }\frac{f^{(k)}(c)}{k!}(x-c{)}^k\\ E_{n+1}& =\frac{f^{(n+1)}(\xi )}{(n+1)!}(h:=x-c{)}^{n+1}\\ |E_{n+1}|\le c{h}^{n+1}& \\ & \end{array}$$

Polynomial Interpolation

$$\begin{array}{rl}v(x)=& \sum _{j=0}^n{c}_j{\varphi }_j(x)\to \text{linearly independent iff}v(x)=0\forall x\to c_j=0\forall j)\\ & \\ \text{Linear system: }& \left[\begin{array}{cccc}{\varphi }_0(x_0)& {\varphi }_1(x_0)& \cdots & {\varphi }_n(x_0)\\ {\varphi }_0(x_1)& {\varphi }_1(x_1)& \cdots & {\varphi }_n(x_1)\\ ⋮& ⋮& \ddots & ⋮\\ {\varphi }_0(x_n)& {\varphi }_1(x_n)& \cdots & {\varphi }_n(x_n)\end{array}\right]\left[\begin{array}{c}{c}_0\\ c_1\\ ⋮\\ c_n\end{array}\right]=\left[\begin{array}{c}{y}_0\\ y_1\\ ⋮\\ y_n\end{array}\right]\end{array}$$ $$\begin{array}{rl}\text{Monomial basis: }& {\varphi }_j(x)=x^j,j=0,1,...,n\to v(x)=\sum _{j=0}^n{c}_j{x}^j\\ & p_n(x_i)=c_0+c_1{x}_i+c_2{x}_i^2+\cdots +c_n{x}_i^n=y_i\\ & \\ X:& \text{Vandermonde matrix}\to \text{det}(X)=\prod _{i=0}^{n-1}\left[\prod _{j=i+1}^n(x_j-x_i)\right]\\ \text{if }& x_i\text{are distinct:}\\ & \bullet \text{det}(X)\ne 0\\ & \bullet X\text{is nonsingular}\\ & \bullet \text{system has unique solution}\\ & \bullet \text{unique polynomial of degree}\le n\text{that interpolates the data}\\ & \bullet \text{can be poorly conditioned, work is }O(n^3)\\ & \end{array}$$ $$\begin{array}{rl}\text{Lagrange basis: }& L_j(x_i)=\{\begin{array}{ll}0& \text{if }i\ne j\\ 1& \text{if }i=j\end{array}\\ & L_j(x)=\prod _{i=0,i\ne j}^n\frac{x-x_i}{x_j-x_i}\\ & p_n(x_i)=\sum _{j=0}^n{y}_j{L}_j(x_i)=\sum _{j=0}^{i-1}{y}_j{L}_j(x_i)+y_i{L}_i(x_i)+\sum _{j=i+1}^n{y}_j{L}_j(x_i)=y_i\\ & \end{array}$$ $$\begin{array}{rl}\text{Newton’s basis: }& {\varphi }_j(x)=\prod _{i=0}^{j-1}(x-x_i),j=0:n\\ & p_n(x_i)=c_0+c_1(x_i-x_0)+\cdots +c_n(x_i-x_0)(x_i-x_1)\cdots (x_i-x_{n-1})=f(x_i)\\ & \end{array}$$ $$\begin{array}{rl}& \text{Divided differences: }f[x_i,\cdots ,x_j]=\frac{f[x_{i+1},\cdots ,x_j]-f[x_i,\cdots ,x_{j-1}]}{x_j-x_i}\\ & \bullet \text{at }x=x_0\text{ then }{c}_0=f(x_0)=f[x_0]\\ & \bullet \text{at }x=x_1\text{ then }{c}_1=\frac{f(x_1)-f(x_0)}{x_1-x_0}=f[x_0,x_1]\\ & \bullet \text{at }x=x_2\text{ then }{c}_2=\frac{f(x_2)-c_0-c_1(x_2-x_0)}{(x_2-x_0)(x_2-x_1)}=\frac{\frac{f(x_2)-f(x_1)}{x_2-x_1}-\frac{f(x_1)-f(x_0)}{x_1-x_0}}{x_2-x_0}=f[x_0,x_1,x_2]\\ & \\ & \therefore \forall x\in [a,b]\exists \xi =\xi (x)\in (a,b):f(x)-p_n(x)=\frac{f^{n+1}(\xi )}{(n+1)!}\prod _{i=0}^n(x-x_i)\\ & \therefore \text{Error: }|f(x)-p_n(x)|\le \frac{M}{4(n+1)}{h}^{n+1}\\ & \text{where: }\\ & \bullet M=ma{x}_{a\le t\le b}|f^{n+1}(t)|\\ & \bullet h=\frac{b-a}{n}\\ & \bullet x_i=a+ih\text{ for }i=0,1,\cdots ,n\end{array}$$

Basic Numeric Integration

$$\begin{array}{rl}& I_f={\int }_a^{b}f(x)dx\approx \sum _{j=0}^n{a}_{j}f(x_j)\text{(quadrature rule)}\\ & \bullet x_0,\cdots ,x_n\text{be distinct points in }[a,b]\\ & \bullet p_n(x)\text{be interpolating polynomial of }f\to {\int }_a^{b}f(x)dx\approx {\int }_a^b{p}_n(x)dx\\ & \bullet \text{Uses Lagrange form: }{\int }_a^{b}f(x)dx\approx \sum _{j=0}^{n}f(x_j){\int }_a^b{L}_j(x)dx=\sum _{j=0}^{n}f(x_j)a_j\\ & \end{array}$$ $$\begin{array}{rl}\text{Trapezoidal rule: }& f(x)\approx p_1(x)=f(x_0)L_0(x)+f(x_1)L_1(x)(n=1,x_0=a,x_1=b)\\ \therefore & I_f={\int }_a^{b}f(x)dx\approx f(a){\int }_a^b\frac{x-b}{a-b}dx+f(b){\int }_a^b\frac{x-a}{b-a}dx\\ & =\frac{b-a}{2}[f(a)+f(b)]\\ \text{Error: }& f(x)-p_1(x)=\frac{1}{2}{f}^{{}^{\prime \prime }}(\xi (x))(x-a)(x-b)\\ \text{then: }& {\int }_a^b(f(x)-p_1(x))dx=\frac{1}{2}{\int }_a^b{f}^{{}^{\prime \prime }}(\xi (x))(x-a)(x-b)dx\\ \text{From MVT: }& \exists \eta \in (a,b):{\int }_a^b{f}^{{}^{\prime \prime }}(\xi (x))(x-a)(x-b)dx=f^{{}^{\prime \prime }}(\eta ){\int }_a^b(x-a)(x-b)dx\\ \therefore & \text{Error of Trapezoidal rule: }{I}_f-I_{trap}=-\frac{f^{{}^{\prime \prime }}(\eta )}{12}(b-a{)}^3\\ & \end{array}$$ $$\begin{array}{rl}\text{Midpoint rule: }& I_f\approx I_{mid}=(b-a)f(\frac{a+b}{2})\\ & \text{Let }m=\frac{a+b}{2}\to f(x)=f(m)+f^{{}^{\prime }}(m)(x-m)+\frac{1}{2}{f}^{{}^{\prime \prime }}(\xi (x))(x-m{)}^2\\ \therefore & I_f={\int }_a^{b}f(x)=(b-a)f(m)+\frac{1}{2}{\int }_a^b{f}^{\prime \prime }(\xi (x))(x-m{)}^{2}dx\\ & \exists \eta \in (a,b):\frac{1}{2}{\int }_a^b{f}^{\prime \prime }(\xi (x))(x-m{)}^{2}dx=\frac{f^{\prime \prime }(\eta )}{24}(b-a{)}^3\\ \therefore & \text{Error of Midpoint rule: }{I}_f-I_{mid}=\frac{f^{{}^{\prime \prime }}(\eta )}{24}(b-a{)}^3\\ & \end{array}$$ $$\begin{array}{rl}\text{Simpson’s rule: }& I_f\approx I_{simp}=\frac{b-a}{6}[f(a)+4f(\frac{a+b}{2})+f(b)]\\ & (p_2(x),n=2,x_0=a,x_1=\frac{a+b}{2},x_2=b)\\ \therefore & \text{Error of Simpson’s rule: }{I}_f-I_{Simpson}=-\frac{f^{(4)}(\eta )}{90}(\frac{b-a}{2}{)}^5,\eta \in (a,b)\\ & \end{array}$$

Composite Numeric Integration

$$\begin{array}{rl}& \bullet \text{subdivide }[a,b]\text{int }r\text{subintervals}\\ & \bullet h=\frac{b-a}{r}\text{length per interval}\\ & \bullet t_i=a+ih\text{for }i=0,1,\cdots ,r\\ & t_0=a,t_r=b\to {\int }_a^{b}f(x)dx=\sum _{i=1}^r{\int }_{t_{i-1}}^{t_i}f(x)dx\\ & \end{array}$$ $$\begin{array}{rl}\text{Composite Trapezoidal rule: }& I_{cf}=\frac{h}{2}[f(a)+f(b)]+h\sum _{i=1}^{r-1}f(t_i)\\ \text{Error: }& I_f-I_{cf}=-\frac{f^{{}^{\prime \prime }}(\mu )}{12}(b-a)h^2\\ \text{Composite Simpson rule: }& I_{cs}=\frac{h}{3}[f(a)+2\sum _{i=1}^{r/2-1}f(t_{2i})+4\sum _{i=1}^{r/2}f(t_{2i-1})+f(b)]\\ \text{Error: }& I_f-I_{cs}=-\frac{f^{(4)}(\zeta )}{180}(b-a)h^4\\ \text{Composite Midpoint rule: }& I_{cm}=h\sum _{i=1}^{r}f(a+(i-1/2)h)\\ \text{Error: }& I_f-I_{cm}=-\frac{f^{{}^{\prime \prime }}(\eta )}{24}(b-a)h^2\\ & \end{array}$$

Linear Least Squares

_Find $c_j$ such that ${\sum }_{k=0}^m(v(x\ast k)-y_k{)}^2=\sum \ast {k=0}^m(\sum \ast {j=0}^n{c}_j{\varphi }_j(x_k)-y_k{)}^2$ is minimised* Conditions: $\frac{\partial \varphi }{\partial a}=0,\frac{\partial \varphi }{\partial b}=0$

$$\begin{array}{rl}\text{Linear fit: }{y}_k& =a{x}_k+b,k=1,\cdots ,m\\ \left[\begin{array}{cc}{\sum }_{k=0}^m{x}_k^2& {\sum }_{k=0}^m{x}_k\\ {\sum }_{k=0}^m{x}_k& m+1\end{array}\right]\left[\begin{array}{c}a\\ b\end{array}\right]& =\left[\begin{array}{c}{\sum }_{k=0}^m{x}_k{y}_k\\ {\sum }_{k=0}^m{y}_k\end{array}\right]\\ p& =\sum _{k=0}^m{x}_k,q=\sum _{k=0}^m{y}_k,r=\sum _{k=0}^m{x}_k{y}_k,s=\sum _{k=0}^m{x}_k^2\\ \to \left[\begin{array}{cc}s& p\\ p& m+1\end{array}\right]\left[\begin{array}{c}a\\ b\end{array}\right]& =\left[\begin{array}{c}r\\ q\end{array}\right]\\ \leftrightarrow A\mathbf{z}& =\left[\begin{array}{cc}{x}_0& 1\\ x_1& 1\\ ⋮& ⋮\\ x_m& 1\end{array}\right]\left[\begin{array}{c}a\\ b\end{array}\right]=\left[\begin{array}{c}{y}_0\\ y_1\\ ⋮\\ y_m\end{array}\right]=\mathbf{f}\text{is overdetermined}\\ & \\ & \end{array}$$ $$\begin{array}{rl}\text{Solving linear system: }r& =b-Ax\\ ||r|{|}_2^2& =\sum _{i=1}^m{r}_i^2=\sum _{i=1}^m(b_i-\sum _{j=1}^n{a}_{ij}{x}_j{)}^2\\ \text{Let }\varphi (x)& =\frac{1}{2}\Vert r{\Vert }_2^2=\frac{1}{2}\sum _{i=1}^m(b_i-\sum _{j=1}^n{a}_{ij}{x}_j{)}^2\\ \text{Conditions}:\frac{\partial \varphi }{\partial x_k}& =0,k=1,\cdots ,n\\ 0& =\sum _{i=1}^m(b_i-\sum _{j=1}^n{a}_{ij}{x}_j)(-a_{ik})\\ \to \sum _{i=1}^m{a}_{ik}\sum _{j=1}^n{a}_{ij}{x}_j& =\sum _{i=1}^m{a}_{ik}{b}_i,k=1,\cdots ,n(\text{equivalent to }{A}^{T}Ax=A^{T}b)\\ & \end{array}$$ $$\begin{array}{rl}{A}^{T}Ax& =A^{T}b\text{is called the normal equations}\\ \text{If }A\text{ has a full-column rank},& \underset{x}{\min }\Vert b-Ax{\Vert }_2\text{has uniq sol:}\\ x& =(A^{T}A{)}^{-1}{A}^{T}b=A^{+}b\\ & \end{array}$$ $$\begin{array}{rl}\text{Adaptive Simpson: find }& Q:|Q-I|\le \text{tol}\\ I& ={\int }_a^{b}f(x)dx=S(a,b)+E(a,b)\\ S_1=S(a,b)& =\frac{h}{6}\left[f(a)+4f\left(\frac{a+b}{2}\right)+f(b)\right]\\ E_1=E(a,b)& =-\frac{1}{90}{\left(\frac{h}{2}\right)}^5{f}^{(4)}(\xi ),\xi \text{ between }a\text{ and }b\\ & \end{array}$$ $$\begin{array}{rl}S=& \text{quadSimpson}(f,a,b,\text{tol})\\ & h=b-a,c=\frac{a+b}{2}\\ & S_1=\frac{h}{6}[f(a)+4f\left(\frac{a+b}{2}\right)+f(b)]\\ & S_2=\frac{h}{12}[f(a)+4f\left(\frac{a+c}{2}\right)+2f(c)+4f\left(\frac{c+b}{2}\right)+f(b)]\\ & {\tilde{E}}_2=\frac{1}{15}(S_2-S_1)\\ & \text{if}|{\tilde{E}}_2|\le \text{tol}\\ & \text{return }Q=S_2+{\tilde{E}}_2\\ & \text{else}\\ & Q_1=\text{quadSimpson}(f,a,c,\text{tol}/2)\\ & Q_2=\text{quadSimpson}(f,c,b,\text{tol}/2)\\ & \text{return }Q=Q_1+Q_2\\ & \end{array}$$

Newton’s Method for Nonlinear equations

$x_{n+1}=x_n-\frac{f(x_n)}{f^{\prime }(x_n)}$

Convergence: if $f,f^{\prime },f^{\prime \prime }$ are continuous in a neighborhood of a root $r$ of $f$ and $f^{\prime }(r)\ne 0$, then $\exists \delta >0:|r-x_0|\le \delta$, then $\forall x_n::|r-x_n|\le \delta ,|r-x_{n+1}|\le c(\delta )|r-x_n{|}^2$ $|e_{n+1}|\le c(\delta )|e_n{|}^2$ (Quadratic convergence, order is 2)

Let $c(\delta )=\frac{1}{2}\ast \frac{{\max }_{|r-x|\le \delta }|f^{\prime \prime }(x)|}{{\min }_{|r-x|\le \delta }|f^{\prime }(x)|}$

For linear system: denote $\mathbf{x}=(x_1,x_2,\cdots ,x_n{)}^T$ and $\mathbf{F}=(f_1,f_2,\cdots ,f_n)$, find ${\mathbf{x}}^{\ast }$ such that $F(x^{\ast })=0$

$$\begin{array}{rl}F(x^{(k)})+F^{\prime }(x^{(k)})(x^{(k+1)}-x^{(k)})& =0\\ F^{\prime }(x^{(k)})& \text{is the Jacobian of }\mathbf{F}\text{at }{x}^{(k)}\\ \text{Let }\mathbf{s}& ={\mathbf{x}}^{(k+1)}-{\mathbf{x}}^{(k)}\\ \therefore F^{\prime }(x^{(k)})s& =-F(x^{(k)})\\ {\mathbf{x}}^{(k+1)}& ={\mathbf{x}}^{(k)}+\mathbf{s}\\ & \end{array}$$

IVP in ODEs.

$$\begin{array}{rl}\text{Given }{y}^{\prime }=f(t,y),y(a)=c,\text{ find }y(t)\text{ for }t\in [a,b]& \\ y^{\prime }& \equiv y^{\prime }(t)\equiv \frac{dy}{dt}\\ \text{System of n first-order: }{y}^{\prime }& =f(t,y),f:\mathbb{R}\times {\mathbb{R}}^n\to {\mathbb{R}}^n\\ & \end{array}$$ $$\begin{array}{rl}\text{Forward Euler’s method (explicit): }{y}_{t_{i+1}}& \approx y(t_i)+hf(t_i,y_{(}{t}_i))\\ \text{where: }h& =\frac{b-a}{N},N>1\\ h& =\text{step size}\\ t_0& =a,t_i=a+ih,i=1,2,\cdots ,N\\ & \end{array}$$

$\text{Backward Euler’s method (implicit): }{y}_{i+1}=y_i+hf(t_{i+1},y_{i+1})$

Non-linear, then apply Newton’s methods

$$\begin{array}{rl}\text{FE Stability: }{y}^{\prime }& =\lambda y,y(0)=y_0\\ \text{Exact sol: }y(t)& =y_0{e}^{\lambda t}\\ \text{FE sol with constant stepsize h: }{y}_{i+1}& =(1+h\lambda )y_i=(1+h\lambda {)}^{i+1}{y}_0\\ \text{To be numerically stable: }h& \le \frac{2}{|\lambda |}\\ & \\ \text{BE Stability: }{y}^{\prime }& =\lambda y,y(0)=y_0\\ |y_{i+1}|& =\frac{1}{|1-h\lambda |}|y_i|\le |y_i|\forall h>0\\ & \end{array}$$

Order, Error, Convergence and Stiffness

$$\begin{array}{rl}\text{Local truncation error of FE: }& d_i=\frac{y(t_{i+1})-y(t_i)}{h}-f(t_i,y(t_i))=\frac{h}{2}{y}^{\prime \prime }({\eta }_i)\text{(q=1)}\\ \text{Local truncation error of BE: }& d_i=-\frac{h}{2}{y}^{\prime \prime }({\xi }_i)\text{(q=1)}\\ & \end{array}$$

$\text{A method of order }q\text{ if}q\text{ is the lowest positive int such that any smooth exact sol of }y(t):{\max }_i|d_i|=O(h^q)$

$$\begin{array}{rl}\text{Global error: }{e}_i& =y(t_i)-y_i,i=0,1,\cdots ,N\\ \text{Consider }{u}^{\prime }& =f(t,u),u(t_{i-1})=y_{i-1},\text{local error: }{l}_i=u(t_i)\\ & \end{array}$$ $$\begin{array}{rl}\text{Convergence: }& \underset{i}{\max }{e}_i=\underset{i}{\max }|y(t_i)-y_i|\to 0\text{ as }h\to 0\\ & \end{array}$$

Stiffness is when the stepsize is restricted by stability rather than accuracy

Runge-Kutta Methods

$$\begin{array}{rl}\text{Implicit trapezoidal: }{y}^{\prime }(t)& =f(t,y),y(t_i)=y_i\\ y_{i+1}& =y_i+\frac{h}{2}[f(t_i,y_i)+f(t_{i+1},y_{i+1})]\\ d_i=O(h^2)& =\frac{y(t_{i+1})-y(t_i)}{h}-\frac{1}{2}[f(t_i,y(t_i))+f(t_{i+1},y(t_{i+1}))]\\ & \\ \text{Explicit trapezoidal: }Y& =y_i+hf(t_i,y_i)\\ y_{i+1}& =y_i+\frac{h}{2}[f(t_i,y_i)+f(t_{i+1},Y)]\\ d_i=O(h^2)& =\frac{y(t_{i+1})-y(t_i)}{h}-\frac{1}{2}[f(t_i,y(t_i))+f(t_{i+1},y(t_i)+hf(t_i,y(t_i)))]\\ & \\ \text{Implicit midpoint: }{y}_{i+1}& =y_i+hf(t_i+h/2,(y_i+y_{i+1})/2)\\ \text{Explicit midpoint: }Y& =y_i+\frac{h}{2}f(t_i,y_i)\\ & \end{array}$$

Classical RK4: based on Simpson’s quadrature rule, $O(h^4)$ accuracy

$$\begin{array}{rl}{Y}_1& =y_i\\ Y_2& =y_i+\frac{h}{2}f(t_i,Y_1)\\ Y_3& =y_i+\frac{h}{2}f(t_i+\frac{h}{2},Y_2)\\ Y_4& =y_i+hf(t_i+\frac{h}{2},Y_3)\\ y_{i+1}& =y_i+\frac{h}{6}[f(t_i,Y_1)+2f(t_i+\frac{h}{2},Y_2)+2f(t_i+\frac{h}{2},Y_3)+f(t_{i+1},Y_4)]\\ & \end{array}$$