MATH 417 Lecture 16

« previous | Thursday, March 20, 2014 | next »

Differential Equations

Initial value problem ODE:

${\begin{cases}{\frac {\mathrm {d} y}{\mathrm {d} t}}=f(t,y(t))\\y(a)=\alpha \end{cases}}$

for $t\in \left[a,b\right]$.

Discretize the interval into $n$ points, and let $h={\frac {b-a}{n}}$

$t_{i}=a+h\,i$ for $i\in 0,1,\ldots ,n$

We are given a method

${\begin{cases}w_{0}=\alpha \\w_{i+1}=w_{i}+h\,\phi (t_{i},w_{i})\end{cases}}$

Last Time: Euler

Let $\phi (t,y(t))=f(t,y(t))$

1. Given a current approximation $w_{i}$, we want ${\begin{cases}y'(t)=f(t,y(t))\\y(t_{i})=w_{i}\end{cases}}$
2. Approximate $f(t,y(t))\approx f(t_{i},w_{i})$
3. $w'(t)=f(t_{i},w_{i})$ by ODE definition
4. $w(t_{i}+1)=w_{i}+h\,f(t_{i},w_{i})$

Error

$\left|y_{i}-w_{i}\right|\leq {\frac {h\,M}{2L}}\,\left(\mathrm {e} ^{L\,\left(t_{i}-a\right)}-1\right)$

In order for $\left|y(b)-w_{n}\right|\in O(h)$, we must have $\left|y''(t)\right|\leq M$ and the IVP must be well-posed:

1. $f(t,y)$ is continuous in $D=\left\{(t,y)~\mid ~a\leq t\leq b,y\right\}$
2. $f$ must be Lipschitz continuous: $\left|f(t,x)-f(t,y)\right|\leq L\,\left|x-y\right|$ for all $(t,x),(t,y)\in D$
3. $D$ must be a convex set: for every two points in the set, the line between them must lie entirely in the set.

Lipschitz Continuity

There is an easy way to verify Lipschitz continuity using partial derivatives:

Rewrite definition as

$\left|{\frac {f(t,x)-f(t,y)}{x-y}}\right|\leq L$

Then $\left|{\frac {\partial f}{\partial y}}(t,\xi )\right|\leq L$

Euler Method Roundoff Error

Each step has a small $\delta _{i}$, so the accumulated error is $\delta$

On most machines, $\delta \approx 10^{-16}$

Then the error equation becomes $\left|y(b)-w_{n}\right|\in O\left(h+{\frac {\delta }{h}}\right)$

Don't make $h$ too small, otherwise accumulated error $\delta$ will be more than the error of $h$

Cost of evaluating is number of steps × local cost

Let error term be bounded by $\epsilon$.

For euler, the cost is ${\frac {b-a}{\epsilon }}\cdot 1$

Euler is an order-1 method

Taylor Series Based Methods

${\begin{aligned}y(t_{i}+h)=y(t_{i})+h\,y'(t_{i})+{\frac {h^{2}}{2}}\,y''(t_{i})+{\frac {h^{3}}{3!}}\,y'''(t_{i})+{\mbox{error}}\\h\cdot T^{(n)}(t,y)=h\,\left(y'(t_{i})+{\frac {h}{2}}\,y''(t_{i})+{\frac {h^{2}}{6}}y'''(t_{i})+\dots \right)\end{aligned}}$

Let's calculate for multiple values of $n$:

n = 1

${\begin{aligned}h\,T^{(1)}(t,y)&=h\,y'(t_{i})\\&=h\,f(t,y)\end{aligned}}$

(this is Euler's method)

Our numerical method is then $w_{i+1}=w_{i}+h\,f(t_{i},w_{i})$

n = 2

${\begin{aligned}h\,T^{(2)}(t,y)&=h\,y'(t_{i})+{\frac {h^{2}}{2}}\,y''(t_{i})\\&=h\,f(t,y)+{\frac {h^{2}}{2}}\,{\frac {\mathrm {d} }{\mathrm {d} t}}\left(f(t,y(t))\right)\\&=h\,f+{\frac {h^{2}}{2}}\left({\frac {\partial f}{\partial t}}+{\frac {\partial f}{\partial y}}\cdot y'(t)\right)\\&=h\,f(t,y)+{\frac {h^{2}}{2}}\left(f_{t}(t,y)+f_{y}(t,y)\cdot f(t,y)\right)\end{aligned}}$

Our numerical method is then $w_{i+1}=w_{i}+h\,f(t_{i},w_{i})+{\frac {h^{2}}{2}}\left(f_{t}(t_{i},w_{i})+f_{y}(t_{i},w_{i})\cdot f(t_{i},w_{i})\right)$

Local Truncation Error

At each step, how far off is our approximation?

$\tau _{i+1}(h)={\frac {y_{i+1}-\left(y_{i}+h\,\phi (t_{i},y_{i})\right)}{h}}$

For Taylor-series-based methods, $\phi =T^{(n)}(t_{i},y_{i})$, and

$\tau _{i+1}(h)={\frac {h^{n+1}}{(n+1)!}}\,y^{(n+1)}(\xi )\cdot {\frac {1}{h}}=C\,h^{n}$

---

n = 3

$h\,T^{(3)}(t,y)=h\,y'(t_{i})+{\frac {h^{2}}{2}}\,y''(t_{i})+{\frac {h^{3}}{6}}\,y'''(t_{i})$

Then

${\begin{aligned}w_{i+1}&=w_{i}+h\,T^{(3)}(t_{i},w_{i})\\&=\ldots \\&=w_{i}+h\,\left(h\,f+{\frac {h^{2}}{2}}\,\left(f_{t}+f_{y}\,f\right)+{\frac {h^{3}}{6}}\,\left(f_{tt}+f_{ty}\,f+f\,\left(f_{yt}+f_{yy}\,f\right)+f_{t}\left(f_{t}+f_{y}\,f\right)\right)\right)\end{aligned}}$

Numerical Integration Methods

Given $w_{i}$, we want to construct $w_{i+1}$

So far, we have seen:

• Euler: $w_{i+1}=w_{i}+h\,f(t_{i},w_{i})$
• Taylor: $w_{i+1}=w_{i}+h\,T^{(n)}(t_{i},w_{i})$

When we evaluate the exact solution to an ODE, we solve

$y_{i+1}=y(t_{i}+h)=w_{i}+\int _{t_{i}}^{t_{i}+h}\,f(t,y(t))\,\mathrm {d} t$

Let the integrand term be represented by $g(t)$ (assume these values are approximated with accuracy $O(h^{2})$):

$y_{i+1}=y(t_{i}+h)=w_{i}+\int _{t_{i}}^{t_{i}+h}\,g(t)\,\mathrm {d} t$

Modified Euler Method

Euler's method is analogous to evaluating the integral by using the left-riemann summ (a horrible approximation)

Let's use the trapezoidal rule instead:

$w_{i+1}=w_{i}+{\frac {h}{2}}\left(g(t_{i})+g(t_{i+1})\right)=w_{i}+{\frac {h}{2}}\left(f(t_{i},w_{i})+f(t_{i+1},w_{i}+h\,f(t_{i},w_{i}))\right)$

If we let $z_{i+1}=w_{i}+h\,f(t_{i},w_{i})$ be the standard Euler method, we can simplify the expression to:

$w_{i+1}=w_{i}+{\frac {h}{2}}\,\left(f(t_{i},w_{i})+f(t_{i+1},z_{i+1})\right)$

Midpoint Method

Let's try using the midpoint method; this requires evaluation at half-time-steps

$w_{i+1}=w_{i}+h\,f\left(t_{i}+{\frac {h}{2}},w_{i}+{\frac {h}{2}}\,f(t_{i},w_{i})\right)$

Let $z_{i+{\frac {1}{2}}}=w_{i}+{\frac {h}{2}}\,f(t_{i},w_{i})$

$w_{i+1}=w_{i}+h\,f\left(t_{i}+{\frac {h}{2}},z_{i+{\frac {1}{2}}}\right)$

Heun's Method

Let's use Gaussian Quadrature with $\int _{0}^{1}f(x)\,\mathrm {d} x\approx {\frac {1}{4}}\,f(0)+{\frac {3}{4}}\,f\left({\frac {2}{3}}\right)$.

Then $w_{i+1}=w_{i}+{\frac {h}{4}}\,\left(f(t_{i},w_{i})+3\,f\left({\frac {t_{i}+{\frac {2h}{3}}}{t_{i}+{\frac {2}{3}}}},{\frac {w_{i}+{\frac {2h}{3}}\,f(t_{i},w_{i})}{w_{i}+{\frac {2}{3}}}}\right)\right)$