Instructor: George Labahn

• Week 1. May 10
• Week 2. May 17
• Week 3. May 25
• Week 4. June 1
• Week 5. June 7
• Week 6. June 15
• Week 7. June 21
• Week 8. June 28
• Week 11. July 19
• Week 12. July 26

Week 1. May 10

prereq:

• calculus (taylor series, derivatives, integrals)
• complex numbers
• big-O
• basic linear algebra
• jupyter

mt (june 24), final (aug 13), 5 assns: 25%, 35%, 40%

Topic 1: floating point number systems

pitfalls of computation

eg. compute $e^{-5.5}$ using 5-digit arithmetic

recall:

$$e^x=1+\frac{x}{1!}+\frac{x^2}{2!}+\frac{x^3}{3!}+...$$

so

$$\begin{aligned} e^{-5.5} &= 1+\frac{-5.5}{1!}+\frac{(-5.5)^2}{2!}+\frac{(-5.5)^3}{3!}+... \\ &\approx 1-5.5+15.125-27.730+...=0.0026363 \end{aligned}$$

$$\begin{aligned} e^{-5.5} &= \frac{1}{e^{5.5}}=(1+\frac{5.5}{1!}+\frac{5.5^2}{2!}+\frac{5.5^3}{3!}+...)^{-1} \\ &\approx (1+5.5+15.125+27.730)^{-1}=0.0040865 \end{aligned}$$

why?

eg. let $I_n=\int_0^1\frac{x^n}{x+\alpha}dx$

note

$$\begin{aligned} I_n&=\int_0^1\frac{x^n+x^{n-1}\alpha-x^{n-1}\alpha}{x+\alpha}dx \\ &=\int_0^1\frac{x^{n-1}(x+\alpha)}{x+\alpha}dx-\int_0^1\frac{\alpha x^{n-1}}{x+\alpha}dx \\ &=\int_0^1x^{n-1}dx-\alpha\int_0^1\frac{x^{n-1}}{x+\alpha}dx \\ &=\frac{1}{n}-\alpha I_{n-1} \end{aligned}$$

where $I_0=\log(\frac{1+\alpha}{\alpha})$.

import numpy as np

def In(alpha, n):
    I = np.log((1 + alpha) / alpha)
    for n in np.arange(1, n+1):
        I = 1./n - alpha*I
    return I

In(0.5, 100)  # 0.006644371123395131
In(2, 100)  # 6058010443603.891 weird

Strange since $I_n=\int_0^1\frac{x^n}{x+\alpha}dx\le\int_0^1\frac{x^n}{1}dx=\frac{1}{n+1}$ if $0\le x\le 1$ and $\alpha>1$. and we should have $I_{100}\le\frac{1}{101}$.

if $I_n^{\mathrm{Exact}}=I_n^\mathrm{Approx}+\mathrm{Err}_n$ then,

$$\begin{aligned} \mathrm{Err}&=I_n^E-I_n^A\\ &=\left(\frac{1}{n}-\alpha I_{n-1}^E\right)-\left(\frac{1}{n}-\alpha I_{n-1}^A\right) \\ &=-\alpha(I_{n-1}^E-I_{n-1}^A) \\ &=-\alpha\mathrm{Err}_{n-1} \\ |\mathrm{Err}_n|&=|\alpha|^n|\mathrm{Err}_0| \end{aligned}$$

• so if $|\alpha|<1$ then $|\mathrm{Err}|\to 0$, the computation is stable.
• if $|\alpha|>1$ then $|\mathrm{Err}|\to \infty$, the computation is unstable.

float point numbers

• computers can represent numbers by fixed point (int) and floating point (real)
• always in normalized form
• issues:
  • density: only keeps a finite number of digits
  • range: only keeps a finite range of exponent

see. CS251

t: significant digits, L: lower bound, U: upper bound

single precision F(b=2, t=24, L=-126, U=127):

• one bit sign, 23 digits for precision, 8 digits for exponent

double precision F(b=2, t=53, L=-1022, U=1023):

• one bit sign, 52 digits for precision, 11 digits for exponent

note implicit leading 1 in precision is used.

special number:

• 0: 00000000 00000000 00000000 00000000
• eg. suppose in F(10, 8, -35, 35)
  • largest number is 0.99999999 * 10^35
  • smallest number: 0.10000000 * 10^-35
  • anything larger: overflow
  • anything smaller: underflow (typically rounded to 0)

defn. if $x=\plusmn0.x_1...x_tx_{t+1}...*\beta^p$, then $\mathrm{fl}(x)=\plusmn0.x_1...x_t*\beta^p$, ie we use truncation. relative error is

$$\delta_x=\frac{\mathrm{fl}(x)-x}{x} \\$$

note $\mathrm{fl}(x)=(1+\delta)x$.

we have

$$\begin{aligned} \delta_x&=\frac{-0.0...0x_{t+1}...}{x_1.x_2...x_tx_{t+1}...}\beta\\ |\delta_x|&\le\frac{0.x_{t+1}...*\beta^{-t}}{1}\beta\\ &\le\beta^{1-t} = E \end{aligned}$$

this quantity is called machine epsilon E.

• it is independent of the numbers involved
• relates to the number of significant digits in the result (t)
• eg. a result is correct to roughly s digits if relative error ≈ 10^-s

it also often defined as the smallest positive number E such that $\mathrm{fl}(1.0+E)>1.0$, and also often called unit round-off error.

if we use rounding we could get $E=\frac{1}{2}\beta^{1-t}$.

arithmetics

to do arithmetics, we do

• put numbers into a computer (fl())
• do normal arithmetics as two reals
• put number back into the computer

eg. suppose x and y are real numbers, then error of addition:

$$\begin{aligned} x\oplus y&=\mathrm{fl}(\mathrm{fl}(x)+\mathrm{fl}(y)) \\ &=(\mathrm{fl}(x)+\mathrm{fl}(y))(1+\delta_z) \end{aligned}$$

$$\begin{aligned} \left|\frac{x+y-x\oplus y}{x+y}\right|&\le \frac{\left|x+y-(x(1+\delta_x)+y(1+\delta_y))(1+\delta_z)\right|}{|x+y|} \\ &=\frac{|x\delta_x+y\delta_y+x\delta_z+y\delta_z+x\delta_x\delta_z+y\delta_y\delta_z|}{|x+y|}\\ &\le \frac{|x||\delta_x|+|y||\delta_y|+|x||\delta_z|+|y||\delta_z|+|x||\delta_x||\delta_z|+|y||\delta_y||\delta_z|}{|x+y|}\\ &=\frac{(|x|+|y|)}{|x+y|}(2E+E^2) \end{aligned}$$

• if x and y are of the same sign then |x+y|=|x|+|y| then the relative error $\le(2E+E^2)$
• if x and y are of the opposite sign and nearly the same size then |x+y| is very small, so the relative error is very large (this is called catastrophic cancellation, why the first way of computing e^-5.5 is wrong)

Week 2. May 17

Polynomial interpolation

interpolation: given n points (x's distinct), find a nice function that go through the points.

• nice means polynomial (of degree at most n-1)
• nice means spline
• points might not be distinct

vandermonde system

with n points $(x_1,y_1),...,(x_n,y_n)$ and distinct $x_i$'s and a polynomial p(x) of degree at most n-1 satisfying $p(x_1)=y_1,...,p(x_n)=y_n$ results in n equations with n unknowns. In matrix form this is

$$\begin{bmatrix} 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{bmatrix} \begin{bmatrix} c_1\\c_2\\...\\c_n \end{bmatrix}= \begin{bmatrix} y_1\\y_2\\...\\y_n \end{bmatrix}$$

the matrix of coefficients is called vandermande matrix V.

we have $|V(x_1,...,x_n)|=\prod_{i<j}(x_i-x_j)$. when x's are distinct then V is nonsigular. the system has a unique solution that always exists.

eg. have four points (1, 1), (-1, 7), (2, 1), (3, -5), find $p(x)=c_1+c_2x+c_3x^2+c_4x^3$.
substitute the points note we have

$$\begin{aligned} p(1)=\quad&c_1&+c_2&+c_3&+c_4&=1\\ p(-1)=\quad&c_1&-c_2&+c_3&-c_4&=7\\ p(2)=\quad&c_1&+2c_2&+4c_3&+8c_4&=1\\ p(3)=\quad&c_1&+3c_2&+9c_3&+27c_4&=-5 \end{aligned}$$

use gaussian elimination to get $p(x)=1-2x+3x^2-x^3$.

lagrange interpolation

lagrange form: given distinct points $(x_1,y_1),...,(x_n,y_n)$ and distinct $x_i$'s define

$$L_i(x)=\frac{(x-x_1)...(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)}=\prod_{j\in[n]-\{i\}}\frac{x-x_j}{x_i-x_j}$$

then the interpolating polynomial is $p(x)=y_1L_1(x)+...+y_nL_n(x)$.

proof.

• for each $i$ degree of $L_i$ is n-1
• $L_i(x_j)=0$ if $i\ne j$
• $L_i(x_i)=1$
• so $p(x_i)=y_i$, p is the unique interpolating polynomial. $\square$

cubic hermite interpolation

given 2 points $(x_1,y_1), (x_2,y_2)$ and two derivative values $s_1,s_2$. find a polynomial $S(x)$ of degree at most 3 with

$$S(x_1)=y_1,S(x_2)=y_2\\ S'(x_1)=s_1,S'(x_2)=s_2$$

then the polynomial is

$$S(x)=a+b(x-x_1)+c(x-x_1)^2+d(x-x_1)^3$$

which is 4 eqs with 4 unknowns. the solution is

$$\begin{aligned} a&=y_1\\ b&=s_1\\ c&=\frac{3y'-2s_1-s_2}{\Delta x}\\ d&=\frac{s_1+s_2-2y'}{\Delta x^2} \end{aligned}$$

where $\Delta x=x_2-x_1,\Delta y=y_2-y_1,y'=\frac{\Delta y}{\Delta x}$.

Week 3. May 25

issues with poly interpolation

• fits point correctly, but might not be so correct elsewhere
• not easy to compute (numerically) for large n

eg.

Spline interpolation

piecewise polynomial interpolation: given n points $(x_1,y_1),...,(x_n,y_n)$ with $x_i$'s distinct. we want to find a piecewise polynomial S(x) s.t.

$$S(x)=\left\{\begin{matrix} S_1(x)&,x_1\le x\le x_2 \\ ...\\ S_{n-1}(x)&,x_{n-1}\le x \le x_n \end{matrix}\right.$$

with each $S_i(x)$ a polynomial (usually same degrees) and $S(x)=y_i$.

linear piecewise polynomial: given n such points, each $S_i(x)$ is line joining $(x_i,y_i)$ to $(x_{i+1},y_{y+1})$ (a poly with degree <= 1).

• problem: not continuous

cubic spline interpolation

given n points $(x_1,y_1),...,(x_n,y_n)$ with $x_i$'s distinct. we want to find a piecewise polynomial S(x) s.t.

$$S(x)=\left\{\begin{matrix} S_1(x)&,x_1\le x\le x_2 \\ ...\\ S_{n-1}(x)&,x_{n-1}\le x \le x_n \end{matrix}\right.$$

with each $S_i(x)$ a polynomial of degree at most 3 satisfying

• $S(x_i)=y_i$
• $S(x)\in C^2$

remark. cubic spline interpolation is not unique.

• n points => n-1 intervals, and each interval has 4 unknowns (deg 3 poly) => 4(n-1) unknowns
• we have 2(n-1)+2(n-2) conditions
  • for every interval, $S(x_i)=y_i$ on both sides => 2(n-1) conditions
  • S' continuous (n-2 joining points, their derivative has to match) => (n-2) conditions
  • S'' continuous => (n-2) conditions
• so 2 more unknowns than conditions $\square$

we can use boundary conditions to let the result be unique.

defn. (natural): $S''(x_1)=S''(x_n)=0$.

defn. (clamped): $S'(x_1)=s_1$ and $S'(x_n)=s_n$ where $s_1,s_n$ are given in advance.

defn. (periodic): $S'(x_1)=S'(x_n)$ and $S''(x_1)=S''(x_n)$ (assuming $y_1=y_n$).

defn. (not-a-knot): $S_1'''(x_2)=S_2'''(x_2)$ and $S_{n-1}'''(x_{n-1})=S_n'''(x_{n-1})$.

• forcing first two polys are same, and last two polys are same

eg. by using cubic hermite for each interval $[x_i,x_{i+1}]$, we have

$$S_i(x)=a_i+b_i(x-x_i)+c_i(x-x_i)^2+d_i(x-x_i)^3$$

and $S_i(x_i)=y_i,S_i(x_{i+1})=y_{i+1},S_i'(x_i)=s_i,S_i'(x_{i+1})=s_{i+1}$, where $s_i,s_{i+1}$ are unknowns. from here we know $S$ interpolates the points and the derivative is continuous. to let second derivative be continuous, we also have $S_i''(x_{i+1})=S_{i+1}''(x_{i+1})$ for $i=1,...,n-2$. note

$$S_i''(x)=2c_i+6d_i(x-x_i),\quad S_{i+1}''(x)=2c_{i+1}+6d_{i+1}(x-x_{i+1})$$

plugging in $x=x_{i+1}$ we get (delta_i is from i to i+1)

$$c_i+3d_i\Delta x_i=c_{i+1}\\ \frac{3 y_{i}^{\prime}-2 s_{i}-s_{i+1}}{\Delta x_{i}}+3 \frac{s_{i}+s_{i+1}-2 y_{i}'}{\Delta x_{i}^{2}}\Delta x_i=\frac{3 y_{i+1}-2 s_{i+1}-s_{i+2}}{\Delta x_{i+1}}\\ \Delta x_{i+1} s_{i}+2\left(\Delta x_{i}+\Delta x_{i+1}\right) s_{i+1}+\Delta x_{i} s_{i+2}=3\left(\Delta x_{i+1} y_{i}^{\prime}+\Delta x_{i} y_{i+1}^{\prime}\right)$$

which are n-2 equations for n unknowns.

suppose we want natural boundary conditions: $S''(x_1)=S''(x_n)=0$.

this is the same as $S_1''(x_1)=S_{n-1}''(x_n)=0$. combine the equations for $\{S_1'',c_1,d_1\}$ and $\{S_{n-1}'',c_{n-1},d_{n-1}\}$ we have additional two eqs.

$$3y_1'=2s_1+s_2,\quad 3y_{n-1}'=2s_n+s_{n-1}$$

note the coefficient matrix looks like diagonal and each row has 2 or 3 non-zero values, and it takes O(n) to determine $s_1,...,s_n$.

parametric curves

given n points $(x_1,y_1),...,(x_n,y_n)$ with $x_i$'s not distinct. we look for a curve $C(t)=(x(t), y(t)),t_1\le t\le t_n$ that goes through the points. steps:

1. get parameter values $t_1,...,t_n$
2. separate x and y into two distinct interpolation problems:
   • find x(t) that interpolates $(t_1,x_1),...,(t_n,x_n)$
   • find y(t) that interpolates $(t_1,y_1),...,(t_n,y_n)$

how to find parameter values?

1. use $t_1=0,t_2=1,...,t_n=n-1$
2. use arc length: $t_1=0,t_{i+1}=t_i+\sqrt{(x_{i+1}-x_i)^2+(y_{i+1}-y_i)^2}$

to plot the curve C(t):

1. determine a set of parameter values $\tau_0,...,\tau_{N},N>>n$
2. plug $\tau_i$ into $x(t),y(t)$ to get $\chi_i,\varphi_i$'s, then plot the points $(\chi_i,\varphi_i)$.

eg.

Week 4. June 1

ODEs

eg. (population model) suppose $p'(t)=rp(t),p(t_0)=p_0$ where $r$ is constant, the solution is $p(t)=p_0e^{r(t-t_0)}$.

general first order ODE

defn. if $\frac{dy}{dt}=f(t,y)$ or $y'(t)=f(t,y(t))$ where $t_0\le t\le t_f$, f is a function of two variables t and y, then $f$ is called system dynamics function.

initial value problem (IVP): let $\mathbf{y}=[y_1,...,y_N]$ and $\mathbf{f}(t,y)=[f_1,...,f_N]$ (function of N+1 variables), $t_0\le t\le t_f$ given

$$\left\{\begin{aligned} \frac{d\mathbf{y}}{dt}&=\mathbf{f}(t,\mathbf{y}) \\ \mathbf{y}(t_0)&=\mathbf{a} \end{aligned}\right.$$

what is $\mathbf{y}$?

eg. (Predator-Prey) $x'(t)=x(t)(a-\alpha y(t)),y'(t)=y(t)(-b+\beta x(t))$ where $a,\alpha,b,\beta$ are constants (the populations interact; when there are many prey, predators increase in number ... until they eat too many).

second order ODE

for $t_0\le t\le t_f$, given

$$\left\{\begin{aligned} \frac{d^2y}{dt^2}&=f(t,y,\frac{dy}{dt})\\ y(t_0)&=u_0\\ \frac{dy}{dt}(t_0)&=v_0 \end{aligned}\right.$$

what is $y$? let $y_1=y,y_2=\frac{dy}{dt}$, we can rewrite it to first order ODE (vector form).

eg. pendulum:

$$\left\{\begin{aligned} &\ddot\theta+\frac{c}{m L} \dot\theta+\frac{g}{L} \sin \theta=0 \\ &\theta(0)=\theta_{0} \\ &\dot\theta(0)=v_{0} \end{aligned}\right.$$

let $y_1=\theta,y_2=\dot\theta$, then this system is equivalent to:

$$\left\{\begin{aligned} &\dot y_1=y_2\\ &\dot y_2=-\frac{c}{mL}y_2-\frac{g}{L}\sin(y_1) \\ &y_1(0)=\theta_{0} \\ &y_2(0)=v_{0} \end{aligned}\right.$$

forward euler method

consider first order ODE

$$\left\{\begin{aligned} y'(t)&=f(t,y(t))\\ y(t_0)&=\alpha \end{aligned}\right.$$

suppose at $t=t_{k-1}$, we have $y_{k-1}\approx y(t_{k-1})$. at $t=t_k$, we use $y_{k-1}$ to compute $y_k$.

using geometric

let $h_{k-1}=t_k-t_{k-1}$ be the time step size. by ode, we have $y'(t_{k-1})=f(t_{k-1},y(t_{k-1}))$, we have approximation $y'(t_{k-1})=\text{slope of y at } t_{k-1}\approx\frac{y_k-y_{k-1}}{h_{k-1}}$. we assume our approximation is correct so we replace $y(t_{k-1})$ by $y_{k-1}$ so the approximation is:

$$\begin{aligned} \frac{y_k-y_{k-1}}{h_{k-1}}&=f(t_{k-1},y_{k-1})\\ y_k&=y_{k-1}+h_{k-1}f(t_{k-1},y_{k-1}) \end{aligned}$$

this requires one function evaluation for each time step.

using taylor expansion

local truncation error: assume timestep is constant h and we have a method computing $y_k$. then $y_k-y(t_k)=O(h^{p+1})$. the bigger p the better.

by taylor, we have $y(t_k)=y(t_{k-1}+h_{k-1})=y({t_{k-1}})+h_{k-1}y'(t_{k-1})+\frac{h^2_{k-1}}{2!}y''(\epsilon_k)$ where $t_{k-1}<\epsilon_k<t_k$. substitute the ode we have $y\left(t_{k}\right)=y\left(t_{k-1}\right)+h_{k-1} f\left(t_{k-1}, y\left(t_{k-1}\right)\right)+O(h^{2})$. use the approximation $y(t_{k-1})\rightarrow y_{k-1},y(t_k)\rightarrow y_k$ we get the result. and $y(t_k)-y_k=O(h^2)$.

eg. compute

$$\left\{\begin{aligned} y'(t)&=\sin t-\cos y(t)\\ y(0)&=1 \end{aligned}\right.$$

we use even interval $h=\frac{1}{N}$ (ie $t_0=0,t_1=h,...,t_N=1$).

1. let $y_0=1$
2. for $k=1,2,...,N$, do $y_k=y_{k-1}+h\cdot(\sin t_{k-1}-\cos y_{k-1})$

eg. compute

$$\left\{\begin{aligned} x'(t)&=2x(t)-3y(t)\\ y'(t)&=tx(t)-3y(t)\\ x(0)&=2\\ y(0)&=1 \end{aligned}\right.$$

where $0<t<1$. use interval $h=\frac{1}{N}$ so $t_k=kh$

1. let $[x_0,y_0]=[2,1]$
2. for $k=1,...,N$ do $[x_k,y_k]=[x_{k-1},y_{k-1}]+h\cdot[2x_{k-1}-3y_{k-1},t_{k-1}x_{k-1}-3y_{k-1}]$

Week 5. June 7

trapezoid method

algo.

$$\begin{aligned} y_0&=\alpha\\ y_k&=y_{k-1}+\frac{h_{k-1}}{2}(f(t_{k-1},y_{k-1})+f(t_k,y_k)) \end{aligned}$$

we have to solve for $y_k$ at each step.

using taylor

note by taylor

$$y(t_k)=y(t_{k-1})+hy'(t_{k-1})+\frac{h^2}{2!}y''(t_{k-1})+\frac{h^3}{3!}y'''(\epsilon_k)\\ y'(t_k)=y'(t_{k-1})+hy''(t_{k-1})+\frac{h^2}{2!}y'''(\alpha)$$

isolate $hy''$ we have $hy''(t_{k-1})=y'(t_k)-y''(t_{k-1})-\frac{h^2}{2!}y'''(\alpha)$, so

$$\begin{aligned} y(t_k)&=y(t_{k-1})+hy'(t_{k-1})+\frac{h}{2!}(y'(t_k)-y'(t_{k-1})-\frac{h^2}{2!}y'''(\alpha))+\frac{h^3}{3!}y'''(\epsilon_k)\\ &=y(t_{k-1})+\frac{h}{2!}(y'(t_{k-1})+y'(t_k))+O(h^3)\\ &=y(t_{k-1})+\frac{h}{2!}(f(t_{k-1},y(t_{k-1}))+f(t_k,y(t_k)))+O(h^3) \end{aligned}$$

use the approximation $y(t_{k-1})\rightarrow y_{k-1},y(t_k)\rightarrow y_k$ we get the result.

modified euler method

consider first order ODE

$$\left\{\begin{aligned} y'(t)&=f(t,y(t))\\ y(t_0)&=\alpha \end{aligned}\right.$$

algo. approximate derivative at any point by taking average slope at starting point and approximated endpoint:

$$\begin{aligned} y_0&=\alpha\\ y_k^*&=y_{k-1}+h_{k-1}f(t_{k-1},y_{k-1})\\ y_k&=y_{k-1}+\frac{h_{k-1}}{2}(f(t_{k-1},y_{k-1})+f(t_k,y_k^*)) \end{aligned}$$

using taylor

as before:

$$\begin{aligned} y(t_k)&=y(t_{k-1})+\frac{h}{2!}(f(t_{k-1},y(t_{k-1}))+f(t_k,y(t_k)))+O(h^3)\\ &=y(t_{k-1})+\frac{h}{2}(f(t_{k-1},y(t_{k-1})+f(t_k,y^*_k)))+\frac{h}{2}(-f(t_k,y^*_k)+f(t_k,y(t_k)))+O(h^3) \end{aligned}$$

since from forward euler we have $y(t_k)-y^*_k=O(h^2)$, so from taylor of $f(t_k,y(t_k))$ as function of y: $f(t_k,y(t_k))=f(t_k,y^*_k)+O(y(t_k)-y^*_k)=f(t_k,y^*_k)+O(h^2)$. hence

$$\begin{aligned} y(t_k)&=y(t_{k-1})+\frac{h}{2}(f(t_{k-1},y(t_{k-1}))+f(t_k,y^*_k))+\frac{h}{2}O(h^2)+O(h^3)\\ &=y(t_{k-1})+\frac{h}{2}(f(t_{k-1},y(t_{k-1}))+f(t_k,y^*_k))+O(h^3) \end{aligned}$$

local truncation error

• LTE for forward euler is is $O(h^2)$
• LTE for trapezoid is $O(h^3)$
• LTE for modified euler is $O(h^3)$

in general

given $y_n=$ formula in $y_n,y_{n-1},...$, to determine the LTE

1. replace all $y_k$ by the actual values $y(t_k)$ (assume formula is correct)
2. compare taylor of $y(t_n)$ with taylor at $t_n$ and at $t_{n-1}$
   • remember $f(t_k,y(t_k))=y'(t_k)$ and has a taylor expansion

eg. what is error of $y_k=y_{k-1}+\frac{h}{4}f(t_{k-2},y_{k-2})+\frac{3h}{4}f(t_k,y_k)$.

1. replace:

   $$\begin{aligned} y(t_k)&=y\left(t_{k-1}\right)+\frac{h}{4} f\left(t_{k-2}, y\left(t_{k-2}\right)\right)+\frac{3 h}{4} f\left(t_{k}, y\left(t_{k}\right)\right) \\ &=y\left(t_{k-1}\right)+\frac{h}{4} y^{\prime}\left(t_{k-2}\right)+\frac{3 h}{4} y^{\prime}\left(t_{k}\right) \end{aligned}$$

2. expand centered around $t_{k-1}$

   $$y'(t_k)=y'(t_{k-1})+hy''(t_{k-1})+\frac{h^2}{2!}y'''(t_{k-1})+...\\ y'(t_{k-2})=y'(t_{k-1})-hy''(t_{k-1})+\frac{h^2}{2!}y'''(t_{k-1})+...\\$$

3. plug in and see where there is a difference between our result and expansion of formula

   $$y\left(t_{k}\right)=y\left(t_{k-1}\right)+hy^{\prime}\left(t_{k-1}\right)+\frac{h^{2}}{2}y^{\prime \prime}\left(t_{k-1}\right)+\frac{h^{3}}{2}y^{\prime \prime \prime}\left(t_{k-1}\right)+...\\ y\left(t_{k}\right)=y\left(t_{k-1}\right)+hy^{\prime}\left(t_{k-1}\right)+\frac{h^{2}}{2}y^{\prime \prime}\left(t_{k-1}\right)+\frac{h^{3}}{6}y^{\prime \prime \prime}\left(t_{k-1}\right)+...\text{(actual formula)}$$

it is O(h^3).

timestamp control

how to determine if numerical method likely gives a correct answer?

• obvious: run the method with timestamp h, then run again with timestamp h/2. compare answers.

we can combine forward euler and modified euler:

$$y_k^\mathrm{FE}=y(t_k)+Ah^2+O(h^3)\\ y_k^\mathrm{ME}=y(t_k)+O(h^3)$$

so we have

$$\text{local error}\approx|Ah^2|=|y_k^\mathrm{FE}-y_k^\mathrm{ME}|$$

in general, at every timestamp consider two rules, one for order p, one for order p+1, the local error is the difference between two methods: |local error|=$|Ah^p|=|y_k^{(1)}-y_k^{(2)}|$.

if error is too big then reduce timestep h by 2 and repeat.

eg. a common method: runge kutta methods (RKF45)

$$\begin{aligned} &k_{1}=h f\left(t_{n}, y_{n}\right) \\ &k_{2}=h f\left(t_{n}+\frac{h}{4}, y_{n}+\frac{k_{1}}{4}\right) \\ &k_{3}=h f\left(t_{n}+\frac{3 h}{8}, y_{n}+\frac{3 k_{1}}{32}+\frac{9 k_{2}}{32}\right) \\ &k_{4}=h f\left(t_{n}+\frac{12 h}{13}, y_{n}+\frac{1932 k_{1}}{2197}-\frac{7200 k_{2}}{2197}+\frac{7296 k_{3}}{2197}\right) \\ &k_{5}=h f\left(t_{n}+h, y_{n}+\frac{439 k_{1}}{216}-8 k_{2}+\frac{3680 k_{3}}{513}-\frac{845 k_{4}}{4104}\right) \\ &k_{6}=h f\left(t_{n}+\frac{h}{2}, y_{n}-\frac{8 k_{1}}{27}+2 k_{2}-\frac{3544 k_{3}}{2565}+\frac{1859 k_{4}}{4104}-\frac{11 k_{5}}{40}\right) \\ &y_{n+1}^{*}=y_{n}+\frac{25 k_{1}}{216}+\frac{1408 k_{3}}{2565}+\frac{2197 k_{4}}{4104}-\frac{k_{5}}{5} \text { with error } O\left(h^{4}\right) \\ &y_{n+1}=y_{n}+\frac{16 k_{1}}{135}+\frac{6656 k_{3}}{12825}+\frac{28561 k_{4}}{56430}-\frac{9 k_{5}}{50}+\frac{2 k_{6}}{55} \text { with error } O\left(h^{5}\right) \end{aligned}$$