What's in a title?

What's a differential equation, then?

E.g. relation between unknown function $y$ and its derivative $y'$

What's a differential equation, then? (2)

Relation between unknown function $y$ and its derivatives $y', y'', y''', \dotsc$

Differential equations: Boundary Value Problems

What's in a title? (2)

Solving differential equations

• Specific classes of DEs can be solved analytically

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

  ⇒ $y(t) = \cos(\pi t) + \dfrac{\sinh(t+1)}{\sinh(2)}$

   

   

• Other problems, e.g.

  • Complicated DEs, e.g. $y''=\epsilon \sinh(\epsilon y)$
  • Large DEs, e.g. physics simulations

  have to be solved numerically

  ⇒ $\left\{ -1, -0.8, -0.4, 0.16, 0.7, 1.1,\dotsc, 0 \right\}$

Solving differential equations, numerically

Several families of numerical methods:

• Multistep methods

  A k-step method uses k previous points to proceed to the next

  Example: 2-step method

• Runge-Kutta methods

  Uses s stages to proceed from the current points to the next

  Example: 2-stage method

Comparing methods

The order of a method

What's in a title? (3)

Enter exponential fitting

• Instead of polynomial fitting spaces...

  $$FS=\left\{ 1, t, t^2, \dotsc, t^p \right\}$$

  consider hybrid fitting spaces

  $$FS=\left\{ 1, t, \dotsc, t^K, \class{hilite}{e^{\tfreqa t}, t e^{\tfreqa t}, \dotsc, t^P e^{\tfreqa t}} \right\}$$

  Very accurate results if

  $$y(t) = a_0 + a_1 t + \dots + a_p t^K + \class{hilite}{b_0 e^{\tfreqa t} + b_1 t e^{\tfreqa t} + \dots + b_P t^P e^{\tfreqa t}}$$

• Constructed with the six-step procedure [Ixaru 1997, Ixaru & Vanden Berghe 2004]

  Classical method ⇒ exponentially fitted method

Example: (EF)Euler

• Coefficients are related:

  $$\dfrac{e^{\tfreqah} - 1}{\tfreqah} = 1+\dfrac{1}{2} \tfreqah + \frac{1}{6} (\tfreqah)^2 + \frac{1}{24} (\tfreqah)^{3}+ \frac{1}{120} (\tfreqah)^4 + \order{(\tfreqah)}{5}$$

  $\tfreqa \to 0$ reveals the classical method

Example: (EF)Euler (2)

Application to the liver model

⚫Classical Euler method - ⚫EF Euler method

Symmetric EF aka trigonometric fitting

• Usually, one considers

  \( \def\pma{\class{hilite}{\pm}} \)

  $$FS=\left\{ 1, t, \dotsc, t^K, e^{\pma \tfreqa t}, t e^{\pma \tfreqa t}, \dotsc, t^P e^{\pma \tfreqa t} \right\}$$

  equivalent with ($\tfreqb=i\tfreqa$)

  $$FS=\left\{ 1, t, \dotsc, t^K, \sin(\tfreqb t), \cos(\tfreqb t), \dotsc, t^P \sin(\tfreqb t), t^P \cos( \tfreqb t) \right\}$$

  Very accurate results if

  $$y(t) = a_0 + a_1 t + \dots + a_p t^K + \class{hilite}{b_0 \sin(\tfreqb t) + c_0 \cos(\tfreqb t) + \dots c_P t^P \cos(\tfreqb t)}$$

  i.e. solutions with oscillatory behaviour of frequency $\mu$

  Trigonometrically fitted Euler method

  $$y_{n+1} = \class{hilite}{\cos(\nfreqa)} y_n + h \class{hilite}{\dfrac{\sin(\nfreqa)}{\nfreqa}} f(t_n, y_n)$$

  $FS=\left\{\sin(\tfreqa t), \cos(\tfreqa t)\right\}$

Example: (EF)Euler (3)

Application to a problem with an oscillatory solution

⚫Classical Euler method - ⚫EF Euler method

• Solution: $y(t)=\cos(3t)$

  ⇒ Very accurate results with EFEuler if $\omega=3$

Chapter 2 - 3
Parameter selection for some EF methods

Parameter values

• EF method contains a free parameter, e.g. $\tfreqa$
• An appropriate value can sometimes be derived
  • Pendulum has a known period
  • LRC circuit has a known resonance frequency

  • From the DE:

    $$y''=-\tfreqb^2 y +(\tfreqb^2 - 1) \sin(x)$$

    Frequency $\approx\tfreqb$

• In general?

  $$y'' = \epsilon \sinh(\epsilon y)$$

  ⇒ systematic approach needed

Numerov's method

$$y_{j-1} + a_0 y_j + y_{j+1}= h^2 \left( b_1 f_{j-1} + b_0 f_j + b_1 f_{j+1} \right)$$

Suitable for problems of the form $y''=f(t, y)$

• Parameter $\omega$: annihilate red factor of leading error term

• Derivatives of $y$:

  • Solve with classical method $\Rightarrow \eta^I$
  • Use $\eta^I$ to approximate $\udiff{y}{6}(t_j), \udiff{y}{4}(t_j), \dotsc$

• Solve red factor for $\omega$, # of solutions depends on FS

Numerov's method: determining parameter value(s)

$\left\{1,t,t^2,t^3, e^{\pm\omega t}\right\}$, 1 solution

$\left\{1,t, e^{\pm\omega t}, t e^{\pm\omega t}\right\}$, 2 solutions

$\left\{e^{\pm\omega t}, t e^{\pm\omega t}, t^{2} e^{\pm\omega t}\right\}$, 3 solutions

⇒ Select the root with the smallest modulus

Numerov methods: results

⚫Classical Numerov method

⚫⚫⚫EF Numerov methods

• EF methods yield a smaller error
• EF methods have order of accuracy 6 instead of 4

Usmani's methods

Usmani's methods: results

• Apply the same parameter selection strategy

⚫Classical Usmani method

⚫EF Usmani methods

• EF methods yield smaller errors
• EF methods behave as order $p+2$ instead of $p$

Chapter 4
Deferred correction with EFMIRK methods

Deferred correction?

• Application of a numerical method

  $$\phi(\class{unknown}{\eta})=0 \Rightarrow \eta = \Delta y + \orderh{p}$$

• Numerical method has a residual R

  $$\phi(\Delta y)=R \neq 0$$

• Knowing R would be excellent

  $$\phi(\class{unknown}{\eta})=R \Rightarrow \eta = \Delta y$$

• But an approximation of R can also be useful

  $$\chi \approx R$$

  $$\phi(\class{unknown}{\eta})=\chi \Rightarrow \eta = \Delta y + \orderh{q}$$

Deferred correction + exponential fitting = ?

• Existing DC technique: combination of two Runge-Kutta methods $\phi$ and $\psi$ (2 x )

  $$\phi(\class{unknown}{\eta^I}) = 0$$

  First solution: $\eta^I$ of order p

  $$\phi(\class{unknown}{\eta^{II}}) = -\psi(\eta^I)$$

  Second solution: $\eta^{II}$ of order q

  Right combination of $\phi$ and $\psi \Rightarrow q>p$

• Exponentially fitted variant

  $$ \begin{align*} \phi[\class{hilite}{\tfreqa}](\class{unknown}{\eta^I}) &= 0 \\ \phi[\class{hilite}{\tfreqa}](\class{unknown}{\eta^{II}}) &= -\psi[\class{freq2}{\tfreqb}](\eta^I) \end{align*} $$

  Use free parameters $\class{hilite}{\tfreqa}$ and $\class{freq2}{\tfreqb}$ to annihilate the leading error term

The cure? A forest!

• Runge-Kutta methods can be analysed by means of rooted trees

  ,

  ,

  , ,

  , , , ,

  $\dotsc$

  and B-series

  $$ \begin{align*} B(\ta, y) = a(\emptyset) &+ \, h \ta(\btree{1}{1})F(\btree{1}{1})(y) + \dfrac{h^2}{2} \ta(\btree{2}{1})F(\btree{2}{1})(y) \\ &+ \dfrac{h^3}{3!} \left[ \ta(\btree{3}{1})F(\btree{3}{1}) + \ta(\btree{3}{2})F(\btree{3}{2}{\vphantom{\Huge{3}}}) \right](y) + \dotsc \end{align*} $$

Option 1: a global approach

Example: $\phi$: Trapezoidal rule + $\psi$: Radau I $(s=2)$

$\tau$	$\phi\class{efc}{[\tfreqa]}$			$-\psi\class{efc}{[\tfreqb]}$
	0	0	$\dfrac{1}{12}\tfreqa^2$	0	0	$\dfrac{1}{12}\tfreqa^2$
	0	0	$\dfrac{1}{12}\tfreqa^2$	0	0	$\dfrac{1}{9}\tfreqa^2-\dfrac{1}{36}\tfreqb^2$
,	$-\dfrac{1}{2}$	0		$-\dfrac{1}{2}$	0
,,	$-1, -1$			$-\dfrac{10}{9}, -1$

• Trapezoidal rule: order 2

• Trapezoidal + Radau: order 3

• EF Trapezoidal + EF Radau:

  Can be order 4 with the right $\tfreqa, \tfreqb$

  $$b_4(\phi[\tfreqa] + \psi[\tfreqb]^{\nu[\tfreqa]}) = \dfrac{1}{72} \left[ \left( \tfreqb^2 - \tfreqa^2 \right) F(\btree{2}{1}) + \frac{1}{3} F(\btree{4}{1}) + F(\btree{4}{2}) \vphantom{\Huge{1}} \right]$$

Option 2: a local approach

Example: $\phi$: Trapezoidal rule + $\psi$: Radau I $(s=2)$

$\tau$	$\phi{[\tfreqa]}$			$-\psi{[\tfreqb]}$
	0	0	$\dfrac{1}{12}\tfreqa^2$	0	0	0	0
	0	0		0	0	$\dfrac{1}{36}\tfreqa^2$
,	$-\dfrac{1}{2}$			0	0
,,				$-\dfrac{1}{9}, -\dfrac{1}{3}$

• Annihilate two red diagonals separately

  $$ \begin{align*} b_3(\phi[\tfreqa]^1) &= \frac{1}{12} \left[ \tfreqa^2 F(\btree{1}{1}) - F(\btree{3}{1}) - F(\btree{3}{2}) \vphantom{\Huge{1}} \right]\\ b_4(\phi[\tfreqb]^1) &= \frac{1}{72} \left[ \tfreqb^{2} F(\btree{2}{1}) + \frac{1}{3} F(\btree{4}{1}) + F(\btree{4}{2}) - F(\btree{4}{3}) - F(\btree{4}{4}) \vphantom{\Huge{\frac{1}{1}}} \right] \end{align*} $$

A local approach: results

⚫Deferred correction

⚫Deferred correction + exponential fitting

• The leading error term is not annihilated but minimised: no order gain
• But the errors are smaller

Chapter 5
Two families of EF methods and their stability functions

Stability of a Runge-Kutta method

• Application to $y'=\lambda y$

  Example for $\lambda=-1,\quad y(t)=e^{-t}$

• Stability plots show stable value of $v:=\tfreqTE h$

  

• Order stars reveal the order of the method: order p ⇒ p+1 grey sectors

Stability functions: Runge-Kutta methods

• Stability function of a RK method

  $$R(v)=\dfrac{a_0 + a_1 v + \dotsc + a_m v^m}{1 + b_1 v + \dotsc + b_n v^n}$$

• For a classical method of order p $$R(v) = e^{v} + \order{v}{p+1}$$

  Example: order 4 ⇒ 5 grey sectors

• For an EF method \begin{equation} \left\{e^{\tfreqa t}, t e^{\tfreqa t},\dotsc, t^{p} e^{\tfreqa t}\right\} \subset FS \Rightarrow \left\{ \begin{aligned} R(\nfreqa) &= e^{\nfreqa} \\ R'(\nfreqa) &= e^{\nfreqa} \\ & \vdots \\ R^{(p)}(\nfreqa) &= e^{\nfreqa} \end{aligned} \right. \end{equation}

  Example: $\nfreqa=\pm ( 2 + 2i )$

  Solid lines:$\left|R(v)\right| = \left|e^{v}\right|$

  Dotted lines:$\arg R(v) = \arg e^{v}$

  ⇒ Intersections:$R(v) = e^{v}$

Stability functions: Obreshkoff methods

• Stability function of a symmetric Obreshkoff method

  $$R(\nfreqTE)=\dfrac{a_0 + a_1 \nfreqTE^2 + \dotsc + a_m \nfreqTE^{2m}}{1 + b_1 \nfreqTE^{2} + \dotsc + b_m \nfreqTE^{2m}}$$

• For a classical method of order p $$R(\nfreqTE) = \cosh(\nfreqTE) + \order{\nfreqTE}{p+2}$$

• For an EF method \begin{equation} \left\{e^{\tfreqa t}, t e^{\tfreqa t},\dotsc, t^{p} e^{\tfreqa t}\right\} \subset FS \class{hilite}{\Leftrightarrow} \left\{ \begin{aligned} R(\nfreqa) &= \cosh(\nfreqa) \\ R'(\nfreqa) &= \sinh(\nfreqa) \\ & \vdots \\ R^{(p)}(\nfreqa) &= \udiff{\cosh}{p}(\nfreqa) \end{aligned} \right. \end{equation}

Construction of specific Obreshkoff methods

• Maximal-order Obreshkoff methods

  \begin{equation} R = P^{m,m}_{\left\{ 0^{p_0}, \pm \nfreqa^{p_1},\dotsc,\pm \nfreqb^{p_n} \right\}}[\cosh] \Leftrightarrow FS= \left\{ \begin{aligned} & 1, t, \dotsc, t^{p_0} \\ & e^{\tfreqa t}, t e^{\pm\tfreqa t},\dotsc, t^{p_1} e^{\pm\tfreqa t} \\ & \hphantom{=}\vdots \\ & e^{\tfreqb t}, t e^{\pm\tfreqb t},\dotsc, t^{p_n} e^{\pm\tfreqb t} \end{aligned} \right\} \end{equation}

• P-stable Obreshkoff methods

  \begin{equation} \left\{ \begin{aligned} & R = \dfrac{1}{2} \left[ P(x) + P(-x) \right] \\ & P:=P^{m,m}_{\left\{ 0^{p_0}, \pm \nfreqa^{p_1},\dotsc,\pm \nfreqb^{p_n} \right\}}[\exp] \end{aligned} \right. \Leftrightarrow FS= \left\{ \begin{aligned} & 1, t, \dotsc, t^{p_0} \\ & e^{\tfreqa t}, t e^{\pm\tfreqa t},\dotsc, t^{p_1} e^{\pm\tfreqa t} \\ & \hphantom{=}\vdots \\ & e^{\tfreqb t}, t e^{\pm\tfreqb t},\dotsc, t^{p_n} e^{\pm\tfreqb t} \end{aligned} \right\} \end{equation}

Conclusion