$\def\udiff#1#2{{#1}^{\left(#2\right)}} \def\stabs#1#2#3{{#1}^{#2,#2}\left(#3\right)} \def\etaf#1#2{\eta_{#1}(#2)} \def\tfreqa{\class{param}{\omega}} \def\tfreqah{\tfreqa h} \def\nfreqa{\tfreqa_h} \def\nfreqsqa{\Omega_h} \def\tfreqb{\class{param}{\mu}} \def\tfreqbh{\class{param}{\tfreqa h}} \def\nfreqb{\class{param}{\tfreqa_h}} \def\nfreqsqb{\Mu_h} \def\tfreqTE{\class{param}{\lambda}} \def\nfreqTE{\class{param}{\tfreqTE_h}} \def\mcO{\mathcal{O}} \def\order#1#2{\mcO\left(#1^{#2}\right)} \def\orderh#1{\order{h}{#1}} \def\Bbseries#1#2{B^\star\left( \boldsymbol{#1},#2 \right)} \def\evalg#1#2#3#4{#3\mathopen{}#1#4#2\mathclose{}} \def\eval#1#2{\evalg{\left(}{\right)}{#1}{#2}}$

# Application of exponential fitting techniques to numerical methods for solving differential equations

## Davy Hollevoet

Promotor: prof. dr. Marnix Van Daele

Vakgroep Toegepaste Wiskunde, Informatica en Statistiek Application of

techniques to

for solving

## What's a differential equation, then?

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

Breakdown of alcohol by liver

• High concentration of alcohol ⇒ liver very active ⇒ fast breakdown
• Low concentration of alcohol ⇒ liver at ease ⇒ slow breakdown
• breakdown rate ∝ concentration of alcohol

$$\cssId{de1}{y'} \cssId{de2}{= \alpha y}$$

change in alcohol concentration
alcohol concentration

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

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

$$F(t, y, y', y'', \dotsc, \udiff{y}{q})=0$$

A q-th order differential equation

Motion of a pendulum

$$y'' + 2 \zeta \omega_0 y' + \omega_0^2 y = 0$$

$\zeta$: damping

$\omega_0$: undamped frequency

⇒ Initial Value Problem

## Differential equations: Boundary Value Problems

Free hanging rope

$$(y'')^2 - \mu^2 (y')^2 = \mu^2$$

Subject to boundary conditions

$$y(-1)=a, \quad y(1)=b$$

Application of

techniques to

for solving

## 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

• Euler's method
$E=\orderh{1}$ 1st order method
• AB2 $E=\orderh{2}$ 2nd order method
• AB3 $E=\orderh{3}$ 3rd order method
How do the methods perform when $h \to 0$?

## The order of a method

#### Global error

• $E=\orderh{p} \Rightarrow$ method of order p

#### Local error

• Error at each step
• $e_k=\orderh{p+1}$
• Accumulation of local errors $e_1, e_2, e_3, \dotsc \Rightarrow E=\orderh{p}$
• Depends on the problem:

$$e_k = c h^{p+1} \class{highlite}{D^{p+1} y(\zeta_k)}$$

• If $$y(t) = a_0 + a_1 t + a_2 t^2 + \dots + a_p t^p,$$ then $e_k=E=0$

• It is said that $$FS=\left\{ 1, t, t^2, \dotsc, t^p \right\}$$

Application of

techniques to

for solving

## 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

#### Classical Euler method

$$y_{n+1} = y_n + h f(t_n, y_n)$$

$FS=\left\{1, t\right\}$

#### Exponentially fitted Euler method

$$y_{n+1} = y_n + h \class{hilite}{\dfrac{e^{\tfreqah} - 1}{\tfreqah}} f(t_n, y_n)$$

$FS=\left\{1, e^{\tfreqa t}\right\}$

• 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

• Solution: $y(t)=e^{-2 t}$

⇒ Very accurate results with EFEuler if $\omega=-2$

## 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 - 3Parameter 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)$

$\left\{1,t,t^2,t^3,t^4,t^5\right\}$

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

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

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

$e_j=h^6 c_0 \udiff{y}{6}(t_j) + \orderh{8}$

$e_j=h^6 c_1 \class{erroperator}{\left[ \udiff{y}{6} - \tfreqa^2 \udiff{y}{4} \right](t_j)} + \orderh{8}$

$e_j=h^6 c_2 \class{erroperator}{\left[ \udiff{y}{6} - 2 \tfreqa^2 \udiff{y}{4} + \tfreqa^4 \udiff{y}{2} \right](t_j)} + \orderh{8}$

$e_j=h^6 c_3 \class{erroperator}{\left[ \udiff{y}{6} - 3 \tfreqa^2 \udiff{y}{4} + 3 \tfreqa^4 \udiff{y}{2} - \tfreqa^6 y \right](t_j)} + \orderh{8}$

• 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 $$y_{j-2} + a_1 y_{j-1} + a_0 y_j + a_1 y_{j+1} + y_{j+2} = h^4 \left[ \class{b2}{b_2} f_{j-2} + \class{b1}{b_1} f_{j-1} + b_0 f_j + \class{b1}{b_1} f_{j+1} + \class{b2}{b_2} f_{j+2} \right]$$

Suitable for problems of the form $\udiff{y}{4} + f(t) y = g(t)$

$M=10$ $M=8$ $M=6$
$\class{b2}{b_2}=0$ $\class{b2}{b_2}=\class{b1}{b_1}=0$
$\orderh{6}$ $\orderh{4}$ $\orderh{2}$

Exponentially fitted variants:

$$\{0^8, \pm\tfreqa^1\}$$ $$\{0^6, \pm\tfreqa^2\}$$ $$\{0^4, \pm\tfreqa^3\}$$ $$\{0^2, \pm\tfreqa^4\}$$ $$\{\pm\tfreqa^5\}$$
$$\{0^6, \pm\tfreqa^1\}$$ $$\{0^4, \pm\tfreqa^2\}$$ $$\{0^2, \pm\tfreqa^3\}$$ $$\{\pm\tfreqa^4\}$$
$$\{0^4, \pm\tfreqa^1\}$$ $$\{0^2, \pm\tfreqa^2\}$$ $$\{\pm\tfreqa^3\}$$

## 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 4Deferred 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!

$\def\btree#1#2{\class{btree}{\class{placeholder}{\cssId{#1-#2}{.}}}} \def\treef#1{{\bf #1}} \def\ta{\treef{a}}$
• 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]$$

## A global approach: results

Deferred correction

Deferred correction + exponential fitting  • We can increase the order of accuracy
• But the errors are not always smaller

## 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 5Two 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

Application of

techniques to

for solving

### differential equations

• Parameter selection strategy for a family of EF Numerov methods
• Parameter selection strategy for a family of EF Usmani methods
• Combination of EF and DC and some approaches to find parameter values
• A few interesting properties of stability function of RK and Obreshkoff methods

/

#