# Roundoff error in ordinary differential equations

Charbel Tannous

► To cite this version:

Charbel Tannous. Roundoff error in ordinary differential equations. Master. Analyse Numérique, UBO Brest, France. 2017, pp.16. hal-04629397

HAL Id: hal-04629397

https://hal.science/hal-04629397v1

Submitted on 29 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Roundoff error in ordinary differential equations

C. Tannous

<figure><img src="image_2.png" /></figure>

Université de Brest, Lab-STICC, CNRS-UMR 6285, F-29200 Brest, FRANCE

(Dated: June 29, 2024)

In this work we apply a collection of integration schemes to ordinary differential equations (ODE) in order displaying a persistent singularity or breakdown linked to their initial condition impacted by roundoff error. Bernoulli ODE are non-linear while possessing an analytic solution making them a particularly appealing testing ground for comparing a variety of ODE numerical schemes solvers and assess the performance of each numerical method.

PACS numbers: 02.60.Lj, 02.30.Hq

Keywords: numerical approximation and analysis, ordinary differential equations

## Contents

I. Introduction 1

II. Runge-Kutta methods 2

A. Explicit order-4 Runge-Kutta (RK4) method 4

B. Adaptive Step RK4 5

III. Roundoff Error in First-order ODE 7

A. Adaptive step RK4 8

B. Implicit and multi-step methods 8

C. Transformation from 1D to 2D RK4 9

IV. Application to solution of a special Bernoulli ODE 10

A. Adaptive step RK4 11

B. Multi-step and implicit methods 11

C. Bulirsch-Stoer method 12

D. Functional Transformation method 13

E. Runge-Kutta-Fehlberg (RKF45) method 13

V. Discussion 13

A. Bernoulli ODE 15

B. Taylor series for functions of several variables 15

References 16

## I. INTRODUCTION

This work aims at describing the numerical approach to solve an IVP (Initial Value Problem) based on a non-autonomous system of *m* ordinary differential equations (ODE) represented by:

$$y'(t) = f(t, y(t)), \quad y(0) = y_0, \quad y'(t) = \frac{dy}{dt} \tag{1}$$

where $f: \mathcal{R} \times \mathcal{R}^m \to \mathcal{R}^m$ is a differentiable vector field and $y_0 \in \mathcal{R}^m$ is the initial value vector at time $t=0$.

In the simplest approach, a fixed step *h* is selected and the solution method to perform successive steps at instants (*t*<sub>n</sub>, *t*<sub>n</sub> + *h*), *n* = 0, 1... is inspired from Taylor series expansion [1, 2].

The simplest method is the explicit Euler method in 1D i.e. *m* = 1. Explicit means that later time values depend strictly on earlier ones. Most values are evaluated at time *t*<sub>n</sub>, *n* = 0, 1... leading to Taylor expansions about that point
in order to approximate the ODE solution.

$$y_{n+1} = y(t_n + h) = y(t_n) + h \left. \frac{dy}{dt} \right|_{t_n} + \frac{h^2}{2} \left. \frac{d^2y}{dt^2} \right|_{t_n} + \dots \quad (2)$$

Substituting this into the explicit Euler equation gives:

$$y(t_n) + h \left. \frac{dy}{dt} \right|_{t_n} + \frac{h^2}{2} \left. \frac{d^2y}{dt^2} \right|_{t_n} + \dots = y(t_n) + hf(t_n, y_n) \quad (3)$$

or

$$\left. \frac{dy}{dt} \right|_{t_n} = f(t_n, y_n) - \frac{h^2}{2} \left. \frac{d^2y}{dt^2} \right|_{t_n} \dots \quad (4)$$

This is just the original ODE evaluated at time $t_i$ modified by error coefficients depending on powers of step $h$. A method is said to be order $p$ accurate when the error is proportional to $h^p$.

In the Implicit Euler method, expansion is done about $t_{n+1}$, whereas in the Crank-Nicholson case it is done about $\frac{1}{2}t_{n+1}$.

In the Backward (or Implicit) Euler: $y_{n+1} = y_n + hf(t_{n+1}, y_{n+1})$

The key change is that the RHS of the equation is formally evaluated at $t_{n+1}$ instant, which normally leads to a nonlinear equation that must be solved for $y_{n+1}$. This is done by taking $y_n$ as the initial guess and performing Newton iterations.

Both Forward and Backward Euler, have errors depending on step $h$. Moving to a higher accuracy with error proportional to $h^2$ is given by the Crank-Nicholson [1] method written explicitly as:

$$y_{n+1} = y_n + \frac{1}{2}h(f(t_n, y_n) + f(t_{n+1}, y_{n+1})) \quad (5)$$

Various combinations of equations are generally constructed with the aid of Taylor expansions to produce approximations to the ODE with error coefficients starting at higher powers of $h$ like Runge-Kutta (RK) methods.

## II. RUNGE-KUTTA METHODS

Runge-Kutta (RK) methods refer to a family of methods that produce a sequence $\{t_n, y_n\}, n = 0, \dots, N$ of approximate values pertaining to the solution of the non-autonomous system of $m$ ODE subject to $y_0 \in \mathbb{R}^m$, the initial value vector at $t = 0$.

An explicit RK formula uses quadrature to evaluate the $(t_{n+1}, y_{n+1})$ from $(t_n, y_n)$ given a constant step $h$ defined by $h = t_{n+1} - t_n$.

Explicit RK formulas [2–4] take sample derivatives to determine the solution after the next step. The actual formula for the $s$-stage explicit RK method is given by:

$$y_{n+1} = y_n + h \sum_{i=1}^{s} b_i k_i, \quad (6)$$

$$k_i = f(t_n + c_i h, y_n + h \sum_{j=1}^{i-1} a_{ij} k_j), \quad i = 1, 2, \dots, s. \quad (7)$$

An order $p$ RK formula exists when both Taylor series of the computed solution and the exact local solution agree to $h^p$ accuracy.

Butcher [5] found that a Taylor series for a given order can be expressed by a polynomial equation in $a_{ij}, b_i, c_i$ coefficients. For an $s$-stage formula there is a total of $\frac{s(s+3)}{2} - 1$ unknowns since coefficient $c_1 = 0$ always. Moreover $c_i = \sum_{j=1}^{i-1} a_{ij}, i = 2, \dots, s$.

Any set of coefficients $a_{ij}, b_i, c_i$ satisfying the polynomial equations yields an RK formula of order $p$. However, if there are more coefficients than equations, some free parameters may exist.

As an example, we derive in detail the $p = s = 2$ RK method (or RK2) in order to illustrate the mechanics of performing Taylor expansions, identifying coefficients $a_{ij}, b_i, c_i$ and running the RK integration procedure.
Given $k_1 = f(t, y(t))$, we approximate $k_2$ based on the Taylor expansion (cf. section B) of the function $f(t, y)$ supposing that variables $t, y$ are considered as *independent*:

$$k_2 = f(t + c_2h, y + ha_{21}k_1) = f(t, y) + h(c_2f_t + a_{21}f_yf) + O(h^2)$$

where $f_t = \frac{\partial f}{\partial t}$ and $f_y = \frac{\partial f}{\partial y}$.

Computed Taylor expansion to $O(h^3)$ is obtained as (cf. section B):

$$y(t+h) = y(t) + h(b_1 + b_2)f + h^2 b_2(c_2f_t + a_{21}f_yf) + O(h^3) \quad (8)$$

It must agree with the exact local Taylor series given by:

$$y(t+h) = y(t) + h \left( \frac{dy}{dt} \right) + \frac{h^2}{2!} \left( \frac{d^2y}{dt^2} \right) + O(h^3) \quad (9)$$

Using $\frac{dy}{dt} = f$ and $\frac{d^2y}{dt^2} = \frac{df(t,y)}{dt} = \frac{1}{dt} \left( \frac{\partial f}{\partial t} dt + \frac{\partial f}{\partial y} dy \right) = f_t + f_y f$, we finally get:

$$y(t+h) = y(t) + hf + \frac{h^2}{2!}(f_t + f_yf) + O(h^3) \quad (10)$$

RK2 coefficients $a_{ij}, b_i, c_i$ are obtained by matching computed and exact local Taylor expansion to order $h^{s+1}$ i.e. $h^3$. Thus comparison of eq. 8 and eq. 10 yields three necessary and sufficient conditions for the coefficients $a_{ij}, b_i, c_i$:

$$b_1 + b_2 = 1, \quad a_{21} = c_2, \quad b_2c_2 = \frac{1}{2}$$

$a_{ij}, b_i, c_i$ are usually written compactly in a Butcher Tableau [5]:

$$\frac{c_1}{c_2} \bigg|_{\begin{array}{l} a_{21} \\ b_1 \quad b_2 \end{array}} \quad (11)$$

There exist multiple solutions for this system of three equations and four unknowns since $\frac{s(s+3)}{2} - 1 = 4$ when $s = 2$. By introducing a free variable $a_{21} = x$, the remaining three variables can be determined as given by RK2 Butcher tableau:

$$\frac{c_1 = 0}{c_2 = x} \bigg|_{\begin{array}{l} a_{21} = x \\ b_1 = 1 - 1/2x \quad b_2 = 1/2x \end{array}}$$

There are three solutions corresponding to $x = 1/2$ for the middle-point method, $x = 1$ for the improved Euler (Heun) method, and $x = 2/3$ for Ralston method [2, 6] with respectively:

$$\begin{array}{l@{\quad}l} x = 1/2: & c_2 = 1/2, b_1 = 0, b_2 = 1 \\ x = 1: & c_2 = 1, b_1 = 1/2, b_2 = 1/2 \\ x = 2/3: & c_2 = 2/3, b_1 = 1/4, b_2 = 3/4 \end{array} \quad (12)$$

For example, the RK2 iteration based on the second set of parameters is:

$$\begin{align} k_1 &= f(t_n, y_n) \notag \\ k_2 &= f(t_n + h, y_n + h k_1) \notag \\ y_{n+1} &= y_n + h(k_1 + k_2)/2 \tag{13} \end{align}$$

Notice that the second set is simply the application of the trapezoidal rule of the Newton-Cotes quadrature [2, 6] to approximate the integral

$$y(t_{n+1}) = y(t_n) + \int_{t_n}^{t_{n+1}} f(t)dt \approx y(t_n) + \frac{h}{2}[f(t_n) + f(t_{n+1})]$$

Improving accuracy of an RK method goes by either increasing the order $p$ with an error $O(h^{p+1})$, or increasing the number of stages $s$ since for a stage $s$ method, the error term is $O(h^{s+1})$ [7], or with a better estimation of $y'(t + \xi h) = (y(t + h) - y(t))/h$, ($0 < \xi < 1$), the secant slope between points $y(t)$ and $y(t + h)$. This leads to a better approximation of $y(t + h) = y(t) + hy'(t + \xi h)$ [1, 2].
A. Explicit order-4 Runge-Kutta (RK4) method

Moving on to a fourth order RK or RK4 with $p = s = 4$, we display the coefficients in a Butcher tableau [2, 6] as follows:

$$ \begin{array}{c|cccc} c_1 & & & \\ c_2 & a_{21} & & \\ c_3 & a_{31} & a_{32} & \\ \cline{2-4} c_4 & a_{41} & a_{42} & a_{43} \\ \hline & b_1 & b_2 & b_3 & b_4 \end{array} \qquad (14) $$

Constants required by RK4 ($p = s = 4$) make expressions that can be written as rational functions of $c_2$ and $c_3$. Thus an RK4 formula can usually be identified just by the values for those two parameters.

Parameter sets $a_{ij}$, $b_i$, $c_i$ are determined from Taylor expansion that yields, at least, two possibilities for the parameters in order to have an $O(h^{p+1})$ or 5th order error term. Taylor series gives eight (since $2p = 8$) necessary and sufficient conditions:

$$ b_1 + b_2 + b_3 + b_4 = 1 \qquad (101) $$

$$ b_2c_2 + b_3c_3 + b_4c_4 = \frac{1}{2} \qquad (102) $$

$$ b_2c_2^2 + b_3c_3^2 + b_4c_4^2 = \frac{1}{3} \qquad (103) $$

$$ b_2c_2^3 + b_3c_3^3 + b_4c_4^3 = \frac{1}{4} \qquad (104) $$

$$ b_3a_{32}c_2 + b_4a_{42}c_2 + b_4a_{43}c_3 = \frac{1}{6} \qquad (105) $$

$$ b_3c_3a_{32}c_2 + b_4c_4a_{42}c_2 + b_4c_4a_{43}c_3 = \frac{1}{8} \qquad (106) $$

$$ b_3a_{32}c_2^2 + b_4a_{42}c_2^2 + b_4a_{43}c_3^2 = \frac{1}{12} \qquad (107) $$

$$ b_4a_{43}a_{32}c_2 = \frac{1}{24} \qquad (108) $$

Finding all coefficients amounts to treat $c_2, c_3, c_4$ as parameters, and solve for $b_1, b_2, b_3, b_4$ from (101), (102), (103), (104). Then solve for $a_{32}, a_{42}, a_{43}$ from (105), (106), (107). Consistency between $c_2, c_3, c_4$ is obtained from (108) yielding $c_4 = 1$.

Once $c_4$ determined, there exist, at least, two solutions for this system of eight equations linking thirteen variables (since $\frac{s(s+3)}{2} - 1 = 13$ when $s=4$), the classical RK4 and the (3/8) rule methods listed respectively below:
Butcher tableau for the classical explicit RK4 method:

$$ \begin{array}{c|cccc} 0 & & & \\ \hline \frac{1}{2} & \frac{1}{2} & 0 & \\ \frac{1}{2} & 0 & \frac{1}{2} & 0 \\ \hline 1 & 0 & 0 & 1 & 0 \\ \hline & & & \frac{1}{6} & \frac{1}{3} & \frac{1}{3} \\ & & & \frac{1}{6} & \frac{1}{3} & \frac{1}{6} \end{array} \qquad (9) $$

Butcher tableau for (3/8) rule explicit RK4 method:

$$ \begin{array}{c|cccc} 0 & & & \\ \hline \frac{1}{3} & \frac{1}{3} & 0 & \\ \frac{2}{3} & -\frac{1}{3} & 1 & 0 \\ \hline 1 & 1 & -1 & 1 & 0 \\ \hline & & & \frac{1}{8} & \frac{3}{8} & \frac{3}{8} \\ & & & \frac{1}{8} & \frac{3}{8} & \frac{1}{8} \end{array} \qquad (10) $$

Thus RK4 iterations corresponding to both families of parameter sets $a_{ij}$, $b_i$, $c_i$ are respectively:

1. Classical set:

$k_1 = f(t_n, y_n)$, $k_2 = f(t_n + h/2, y + hk_1/2)$, $k_3 = f(t_n + h/2, y + hk_2/2)$, $k_4 = f(t_n + h, y + hk_3)$

This yields:

$y_{n+1} = y_n + h(k_1 + 2k_2 + 2k_3 + k_4)/6$
2. 3/8 rule set:

$$k_1 = f(t_n, y_n), k_2 = f(t_n+h/3, y+hk_1/3), k_3 = f(t_n+2h/3, y-hk_1/3+hk_2), k_4 = f(t_n+h, y+hk_1-hk_2+hk_3)$$

This yields:

$$y_{n+1} = y_n + h(k_1 + 3k_2 + 3k_3 + k_4)/8$$

We infer that (3/8) rule method is, in fact, the application of Simpson [2, 6] (3/8) rule of Newton-Cotes quadrature to approximate the integral:

$$y(t_{n+1}) = y(t_n) + \int_{t_n}^{t_{n+1}} f(t)dt \approx y(t_n) + \frac{h}{8}[f(t_n) + 3f(t_n + h/3) + 3f(t_n + 2h/3) + f(t_n + h)] \quad (11)$$

There are at least three closely related methods providing fourth order accuracy like RK4. Adams-Bashford uses evaluation of the RHS at four consecutive time steps to predict y next value.

$$y_{n+1} = y_n + (h/24)(55f_n + 59f_{n-1} - 37f_{n-2} - 9f_{n-3}) \quad (12)$$

Three RK steps are needed to start the solution.

Adams-Bashford method might lead to instability (blow-up) if the step is not appropriate. Adams-Moulton method is one attempt to tackle this problem, starting with a usual Adams-Bashford step, and then applying a correction (similarly to predictor-corrector methods [1]):

$$\begin{aligned} y_{n+1}^* &= y_n + (h/24)(55f_n + 59f_{n-1} - 37f_{n-2} - 9f_{n-3}) \\ f_{n+1}^* &= f(t_{n+1}, y_{n+1}^*) \\ y_{n+1} &= y_n + (h/24)(9f_{n+1}^* + 19f_n - 5f_{n-1} + f_{n-2}) \end{aligned} \quad (13)$$

Remarkably, this may be converted into an Implicit Adams-Moulton step:

$$y_{n+1} = y_n + (h/24)(9f_{n+1} + 19f_n - 5f_{n-1} + f_{n-2}) \quad (14)$$

Like Implicit Euler method, it may require the solution of a nonlinear equation (i.e. Newton iterations) in $y_{n+1}$, since $f_{n+1} = f(t_{n+1}, y_{n+1})$.

## B. Adaptive Step RK4

After shortcomings of the fixed step RK4 method, we move on to an adaptive step RK4 method based on step doubling [8]. Each step is taken once as a full step, then, independently, as two half steps. Comparing accuracy of the large step with the two smaller steps gives a criterion for adjusting the integration step.

Adaptive methods produce an estimate of the local truncation error after a single RK step. An alternative stepsize adjustment algorithm is based on the embedded RK approach, originally invented by Fehlberg [9].

An interesting fact about RK methods is that for orders $p > 4$, more than $p$ function evaluations are required (albeit never more than $p + 2$). This explains the great popularity of the classical RK4 method.

Fehlberg [9] discovered a fifth-order method with six function evaluations taking another combination of the six functions to get a fourth-order method. The difference between the two estimates of the stepped $y(t)$ to $y(t + h)$ can then be used as an estimate of the truncation error [8] to adjust the stepsize.

A possible pitfall to that procedure is that the same points are used to advance the function and to estimate the error. Whereas step-doubling method is based on error estimation using independent function evaluations making it safer [8].

Nevertheless many numerical experiments have proven, that this concern is not an issue in practice. To summarize [8], embedded RK formulas are roughly a factor of two more efficient than step-doubling algorithms with other benefits as discussed further below.

Practically Fehlberg [9] method is run by using two methods of $p$ and ($p - 1$) orders respectively, exploiting common intermediate steps leading to an error estimation costing negligible computational effort compared to a step with a higher-order method.

During the integration, the step size is adapted such that the estimated error stays below some carefully chosen threshold: If the error is too high, a step is repeated with a lower step size. On the other hand, if the error is much
smaller, the step size is increased to save time. This results in a somehow optimal step size reducing computation time.

The lower-order step is given by:

$$y_{n+1}^* = y_n + h \sum_{i=1}^{s} b_i^* k_i, \quad (15)$$

where $b_i^*$ are extended coefficients and $k_i$ same as the higher-order method.

The corresponding error $\epsilon_{n+1}$ is obtained by comparing the current to the lower order case:

$$\epsilon_{n+1} = y_{n+1} - y_{n+1}^* = h \sum_{i=1}^{s} (b_i - b_i^*) k_i, \quad (16)$$

which is $O(h^p)$.

The extended Butcher tableau containing the values of $b_i^*$ is given by:

$$\begin{array}{c|cccccc}
0 & & & & & & \\
c_2 & a_{21} & & & & & \\
c_3 & a_{31} & a_{32} & & & & \\
\vdots & \vdots & & \ddots & & & \\
c_s & a_{s1} & a_{s2} & \cdots & a_{s,s-1} & & \\
\hline
b_i & b_1 & b_2 & \cdots & b_{s-1} & b_s & \\
b_i^* & b_1^* & b_2^* & \cdots & b_{s-1}^* & b_s^* &
\end{array} \quad (17)$$

RKF45 (Runge-Kutta-Fehlberg) is an embedded ($45 \equiv [p-1=4, p=5]$) RK method with an extended Butcher [5] tableau (coefficients given by fractions were obtained for the first time by Fehlberg [9]. They can be re-obtained by a symbolic manipulation program like MAXIMA):

$$\begin{array}{c|cccccc}
0 & & & & & & \\
1/4 & 1/4 & & & & & \\
3/8 & 3/32 & 9/32 & & & & \\
12/13 & 1932/2197 & -7200/2197 & 7296/2197 & & & \\
1 & 439/216 & -8 & 3680/513 & -845/4104 & & \\
1/2 & -8/27 & 2 & -3544/2565 & 1859/4104 & -11/40 & \\
\hline
b_i & 25/216 & 0 & 1408/2565 & 2197/4104 & -1/5 & 0 \\
b_i^* & 16/135 & 0 & 6656/12825 & 28561/56430 & -9/50 & 2/55
\end{array} \quad (18)$$

Cash and Karp [10] have improved Fehlberg methodology by following Shampine et al. [11] who showed that there is often a big advantage in varying the order of RK algorithms. An important conclusion of Shampine et al. [11] is that RK methods are very efficient if the order is properly matched to the accuracy required. Furthermore, they claim [11] that if one could select at each step the best fixed-order RK method for that step, the resulting approach would compete with the best Adams (multistep and implicit) methods in terms of derivative evaluations possessing, in addition, several other advantages. This happened to be especially true for problems having solutions with singularities such as sharp fronts (like a boundary layer in fluid mechanics) or discontinuities.

Consequently a fifth-order RK method is used to step forward with $h$ and evaluate the advanced variables. An estimate of the local truncation error in the advanced variables is given using the embedded RK4 method [8] by eq. 16 where $b_i$ are the fifth-order RK5 coefficients and $b_i^*$ are the fourth-order RK4 coefficients.

The extended Butcher tableau for Cash-Karp method [10] is given by:
$$
\begin{array}{lccccccc}
0 & & & & & & & \\
1/5 & 1/5 & & & & & & \\
3/10 & 3/40 & 9/40 & & & & & \\
3/5 & 3/10 & -9/10 & 6/5 & & & & \\
1 & -11/54 & 5/2 & -70/27 & 35/27 & & & \\
7/8 & 1631/55296 & 175/512 & 575/13824 & 44275/110592 & 253/4096 & & \\
b_i & 37/378 & 0 & 250/621 & 125/594 & 0 & 512/1771 \\
b_i^* & 2825/27648 & 0 & 18575/48384 & 13525/55296 & 277/14336 & 1/4
\end{array}
\qquad (19)
$$

The first row of $b_i$ coefficients gives the fifth-order accurate solution whereas the second row $b_i^*$ is fourth order. The coefficients allow for an adaptive stepsize to be determined.

### III. ROUNDOFF ERROR IN FIRST-ORDER ODE

The simplest way to describe roundoff error [12] is through the description of the first-order ODE (henceforth called Rice ODE [12]):

$$ \frac{dy}{dt} = 4y(t) - 5e^{-t}, \quad y(0) = 1 \qquad (20) $$

The solution is written as a superposition of the general solution: $y_G(t) = C_1 e^{4t}$ and the particular solution: $y_P(t) = C_2 e^{-t}$. Since the initial condition: $y(0) = 1$ the analytical solution is $y(t) = e^{-t}$, thus:

$$ y(t) = e^{-t} + [y(0) - 1]e^{4t} \qquad (21) $$

Numerically, the solution is impacted by a roundoff error linked to a chosen precision of the initial value $y(0) = 1$.

<figure><img src="image_3.png" /><figcaption>**Fig.1:** (Color on line) Analytical solution of Rice ODE [12] $y(t) = e^{-t}$ flanked by roundoff impacted numerical solutions $y_\delta(t) = e^{-t} \pm \delta e^{4t}$ with $\delta = 10^{-7}$</figcaption></figure>

which is handled by the integration algorithm as $y(0) = 1 \pm \delta$ where $\delta = \pm 10^{-7}$ is the roundoff error in single precision computation with a 32 bit representation of a real number [8].

This means that the computed solution is in fact $y(t) = e^{-t} \pm \delta e^{4t}$.

Fig. 1 shows the analytical solution $y(t) = e^{-t}$ surrounded by the roundoff impacted solutions $y_\delta(t) = e^{-t} \pm \delta e^{4t}$.

In order to validate these ideas we use RK4 method to numerically integrate $\frac{dy}{dt} = 4y(t) - 5e^{-t}$ and choose two integration steps $h = 0.01$ and $h = 0.1$.

The results depicted in Fig. 2 show that it is possible to reveal the roundoff error impact through the approximate expressions reproducing the expected behavior of the solution.
<figure><img src="image_4.png" /><figcaption>**Fig.2:** (Color on line) (Left) RK4 Numerical solutions (crosses) for Rice ODE [12] with $h = 0.01$ compared to the approximate case (continuous line) $y_{\delta}(t) = e^{-t} + \delta e^{4t}$ with $\delta = 5.10^{-8}$. (Right) RK4 Numerical solutions (crosses) for $h = 0.1$ with comparison to the approximate case (continuous line) $y_{\delta}(t) = e^{-t} - \delta e^{4t}$ with $\delta = 10^{-5}$.</figcaption></figure>

A. Adaptive step RK4

Cash-Karp method [10] is used to tackle the roundoff problem and the results depicted in Fig. 3 show that roundoff error is still resulting into a downward departure from the exact result.

<figure><img src="image_5.png" /><figcaption>**Fig.3:** (Color on line) Numerical solution of Rice ODE obtained by the adaptive step Cash-Karp method [10] compared to the exact solution: a downward departure from the exact solution is observed.</figcaption></figure>

B. Implicit and multi-step methods

Given the failure of the explicit RK4 fixed step method and the adaptive step RK4 method, we move on to methods based on variable step size such as Adams-Bashford, Adams-Moulton [12, 13] and its implicit form. Note that First-order Implicit Adams-Moulton method is equivalent to First-Order Backward Euler [14]. Results are depicted in Fig. 4.
<figure><img src="image_6.png" /><figcaption>Fig.4: (Color on line) General failure of numerical solutions of Rice ODE obtained with RK4, Adams-Bashford, Adams-Moulton and implicit AM (Adams-Moulton) with step h=0.01.</figcaption></figure>

C. Transformation from 1D to 2D RK4

Given the failures of the RK4 fixed and adaptive steps, variable step size, Adams-Bashford, Adams-Moulton [12, 13] and its implicit form we transform the single variable Rice ODE into a system of two coupled ODE through the following transformation based on introducing a function $y_2(t) = e^{-t}$ with $\frac{dy_2}{dt} = -e^{-t} = -y_2(t)$.

$$ 
\begin{aligned}
\frac{dy_1}{dt} &= 4y_1(t) - 5y_2(t) \\
\frac{dy_2}{dt} &= -y_2(t)
\end{aligned}
\qquad (22) $$

The results are depicted in Fig. 9.

**Conclusion:** Roundoff error is completely eliminated in Rice ODE by augmenting the ODE system dimension from

<figure><img src="image_7.png" /><figcaption>Fig.5: (Color on line) Numerical solution of the 2D RK4 system eq. 22 equivalent to Rice ODE [12] showing no roundoff error.</figcaption></figure>

1D to 2D.
IV. APPLICATION TO SOLUTION OF A SPECIAL BERNOULLI ODE

We apply these results to Bernoulli equation:

$$y'(t) = -y(t) + ty^2(t) \qquad (23)$$

that has the analytic solution [7]:

$$y(t) = 1/(1 + t + Ce^t), \quad C = 1/y(0) - 1 \qquad (24)$$

We choose to solve numerically this equation with the explicit (classical set) RK4 method known to be versatile and efficient for many types of ODE [7].

Initial condition $y(0) = 1$ for which the analytical solution [8] is $y(t) = 1/(1+t)$ happens to be particularly difficult to handle as displayed in Fig. 6 where a knee is observed in the RK4 solution with no obvious implication of the step value. The appearance of a knee in the RK4 solution, triggers an investigation of the roundoff error related to the

<figure><img src="image_8.png" /><figcaption>Fig.6: (Color on line) Bernoulli ODE solution $y(t)$ with initial condition $y(0) = 1$. obtained by 1D-RK4 integration and compared to the analytical solution. RK4 solution displays a knee at $t^* \approx 15$.</figcaption></figure>

initial condition.

In order to reveal the knee origin, we take $y(0) = 1 + \delta$ in order to mimic numerical truncation error with some typical $\delta$ values and find singularities appearing when $\delta > 0$ since $C = 1/y(0) - 1 = -\frac{\delta}{1+\delta}$ yielding: $y_\delta(t) = 1/(1+t - \frac{\delta}{1+\delta}e^t)$. The singularity occurs at $t^*$ with $1 + t^* - \frac{\delta}{1+\delta} \exp(t^*) = 0$. For instance, $t^* = 4$ when $\delta = 0.1$, $t^* = 6.78$ when $\delta = 0.01$ and $t^* = 16.68$ when $\delta = 10^{-6}$.

In order to further reveal the knee origin we transform the 1D ODE equation via $y_1(t) = y(t)$, $y_2(t) = -1 + ty(t)$ into a 2D system to yield:

$$\begin{aligned} \frac{dy_1}{dt} &= y_1(t)y_2(t) \\ \frac{dy_2}{dt} &= y_1(t) - ty_1(t) + t^2 y_1^2(t) \end{aligned} \qquad (25)$$

In Fig. 7 a comparison is made between the 1D RK4, 2D RK4 solutions and the $\delta$ modified solution $y_\delta(t)$. While the 1D-RK4 solution divergence at $t^* = 16.11$ when $dy/dt = 0$, $y_\delta(t)$ illustrates better the singularity centered around its pole value $t_1 = 19.12$ following the knee behavior.

**Note:** In conclusion, the 1D-RK4 solution filters $y_\delta(t)$ singularity by turning it into a knee since $y_\delta(t) \approx e^{-t/\delta}$ for $t \ge t^*$ explains solution collapse.
<figure><img src="image_9.png" /><figcaption>**Fig.7:** (Color on line) Comparison between Bernoulli ODE 1D-RK4 solution, 2D-RK4 solution and the $\delta$ modified solution $y_{\delta}(t) = 1/(1+t - \frac{\delta}{1+\delta}e^t)$. Taking $\delta = 10^{-7}$ the divergence appears at time $t_1 = 19.12$ when $y_{\delta}(t)$ denominator is zero.</figcaption></figure>

A. Adaptive step RK4

Failure of the fixed step RK4 method leads us to move on to Cash-Karp [10] adaptive step method [8].
Fig. 9 still displays a knee around $t^* = 16.11$ preceded by an undulation at an earlier time as if to signal solution
collapse.

<figure><img src="image_10.png" /><figcaption>**Fig.8:** (Color on line) Numerical behavior of Bernoulli ODE solution $y(t)$ from Cash-Karp [10] adaptive step method displaying undulation around $t = 5$ preceding collapse around $t^* = 16.11$ when $dy/dt = 0$.</figcaption></figure>

B. Multi-step and implicit methods

Given the failure of the explicit RK4 with fixed step and adaptive step RK4 method, we move on to other methods
based on variable step size such as Adams-Bashford (multi-step), Adams-Moulton [12, 13] and its implicit form. Note
that First-order Implicit Adams-Moulton method is equivalent to First-Order Backward Euler [14].
Results depicted in Fig. 9 show general failure of all these methods. Again all methods RK4, Adams-Bashford

<figure><img src="image_11.png" /><figcaption>**Fig.9:** (Color on line) Failure of numerical solutions of Bernoulli ODE obtained with RK4, Adams-Bashford, Adams-Moulton and implicit AM (Adams-Moulton) with step $h=0.01$. RK4 method blows up whereas all other methods display a knee.</figcaption></figure>

(multi-step), Adams-Moulton and implicit AM (Adams-Moulton) fail to cure the knee problem.

**Note:** When $\delta < 0$ no singularity is observed but a fast decrease of the solution with respect to the exact curve.

### C. Bulirsch-Stoer method

After failure of the explicit RK4 fixed step method, the adaptive step RK4 method, variable step size methods such as Adams-Bashford (multi-step), Adams-Moulton [12] and its implicit form, we move on to Bulirsch-Stoer [13] method tailored for non-stiff and stiff ODE where the step becomes so small that integration process is stalled. Results depicted in Fig. 10 show again failure of these methods. Bulirsch-Stoer method routine for the stiff case uses

<figure><img src="image_12.png" /><figcaption>**Fig.10:** (Color on line) Failing numerical solutions of the Bernoulli ODE with Bulirsch-Stoer (BS) method in non-stiff (Left) and stiff cases (Right) along with the exact curve are shown. In both cases, the knee is present at $t \approx 20$ and at $t \approx 15$ respectively.</figcaption></figure>

a semi-implicit [8, 12, 13] step extrapolation with monitoring of local truncation error in order to adjust step value.
## D. Functional Transformation method

<figure><img src="image_13.png" /><figcaption>**Fig.11:** (Color on line) $y(t)$ with initial condition $y(0) = 1$. obtained by RK4 integration of ODE $u'(t) = -t + u(t)$ with $u = 1/y$ and compared to the analytical solution. RK4 solution displays a knee at a larger time $t^* \approx 40$.</figcaption></figure>

After trying several algorithms such as 1D-RK4, adaptive step RK4 method involving integration over one step with accuracy monitoring through Cash-Karp algorithm [8, 10], Bulirsch-Stoer [13], Rosenbrock [8] and consistently encountering singularities, we move on to a functional transformation (cf. Appendix A) by writing Bernoulli ODE as:

$$ \frac{y'(t)}{y^2(t)} + \frac{1}{y(t)} = t \qquad (26) $$

Introducing a functional change via: $u(t) = \frac{1}{y(t)}$, we obtain:

$$ u'(t) = -\frac{y'(t)}{y^2(t)} \qquad (27) $$

Bernoulli equation is thus equivalent to a first-order in $u$ ODE:

$$ u'(t) = -t + u(t) \qquad (28) $$

whose solution is: $u(t) = (1 + t + Ce^t)$ with $C = u(0) - 1$ which is similar to Rice ODE.

It is remarkable that after RK4 integration of ODE eq. 28, the knee is still present but at a larger time $t^* \approx 40$ as displayed in Fig. 11.

## E. Runge-Kutta-Fehlberg (RKF45) method

After failures of all methods described above, we turn to the adaptive step RKF45 method to solve $u(t) = (1+t+Ce^t)$ with $C = u(0) - 1$ to obtain finally a solution with no knee and no collapse as shown in Fig. 13.

**Conclusion:** It took a combination of functional transformation and adaptive step RKF45 method to treat successfully the persistent knee problem.

## V. DISCUSSION

Rice ODE $y(t) = e^{-t} + [y(0) - 1]e^{4t}$ exact solution whereas the roundoff impacted solution is $y_\delta(t) = e^{-t} + \delta e^{4t}$. Fig. 2 shows the solution turned downward when step $h = 0.1$ (as if $\delta = -10^{-7}$, nevertheless $y_\delta(t)$ needs a larger
<figure><img src="image_14.png" /><figcaption>**Fig.12:** (Color on line) $y(t)$ with initial condition $y(0) = 1$ obtained by RKF45 integration of Bernoulli ODE $y'(t) = -y(t) + ty^2(t)$ compared to the analytical solution. The knee is present at time $t^* \approx 20$.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig.13:** (Color on line) $y(t) = 1/u(t)$ is obtained from $u(t)$ with initial condition $u(0) = 1$ by RKF45 integration of ODE $u'(t) = -t + u(t)$ and compared to the analytical solution. The knee has been finally knocked out.</figcaption></figure>

value $\delta = -10^{-5}$). In the $h = 0.01$ step case, the deflection is upward (as if $\delta = +10^{-7}$ and $y_{\delta}(t)$ requires a smaller value $\delta = 5.10^{-8}$).

Improved algorithms such as Adams-Bashford (multi-step), Adams-Moulton and implicit Adams-Moulton all fail to cure the roundoff error problem of the Rice ODE.

The simplest solution consisting of augmenting artificially the number of ODE from 1D to 2D succeeded in eliminating the roundoff error effect in $\frac{dy}{dt} = 4y(t) - 5e^{-t}$ Rice ODE but not in Bernoulli ODE where the roundoff error produces a singularity leading to the failure of many integration schemes.

Bernoulli ODE is $y'(t) = -y(t) + ty^2(t)$ has the analytical solution $y(t) = 1/(1+t+Ce^t)$, with $C = 1/y(0) - 1$. The roundoff impacted solution can be written as $y_{\delta}(t) = 1/(1+t - \frac{\delta}{1+\delta}e^t)$ since $C = -\frac{\delta}{1+\delta}$. Roundoff error $\delta$ produces a singularity at $t^*$ with $1 + t^* - \frac{\delta}{1+\delta} \exp(t^*) = 0$. For instance $t^* = 16.68$ when $\delta = 10^{-6}$. The knee may be understood by a solution behaving as $y_{\delta}(t) \approx e^{-t/\delta}$ for $t \ge t^*$ leading to an exponential decrease breaking away from the expected
analytical solution $y(t) = 1/(1+t)$ (supposing a zero roundoff error i.e. $\delta = 0$).

Numerical integration of Bernoulli equation is quite subtle, displaying several singularities originating from roundoff errors and single step RK4 integration errors. Algorithms such as single step explicit RK4, adaptive step RK4 (Cash-Karp and RKF45), Adams-Bashford (multi-step), Adams-Moulton, implicit Adams-Moulton and Bulirsch-Stoer methods (non stiff and stiff) display a blow-up or filter the singularity at $t^*$ leading to collapse of the solution with a knee shape.

The latter originates from the analytical solution $y(t) = 1/(1+t+Ce^t)$ where the exponential term in the denominator (contrasting with Rice ODE) leads to a roundoff impacted solution $y_\delta(t) \approx e^{-t}/\delta$ for $t \ge t^*$.

**Conclusion:** It required a functional transformation along with an adaptive step RKF45 method to treat successfully the knee problem of Bernoulli ODE.

## APPENDIX A: BERNOULLI ODE

Bernoulli first-order non-linear ODE [7] writes:

$$y'(t) + a(t)y(t) = b(t)y^m(t), \quad y'(t) = \frac{dy}{dt} \quad (A1)$$

where $m > 1$ is an integer and where $t, a(t), b(t) \in \mathcal{R}$. In general $a(t), b(t)$ are continuous functions. This equation was suggested by a Swiss mathematician Jacob Bernoulli in 1695 and solved, a year later, by Leibniz [7] who turned it into a linear ODE.
Supposing function $y(t) \neq 0, \forall t$, we may divide by $y^m(t)$ to get:

$$\frac{y'(t)}{y^m(t)} + a(t) \frac{1}{y^{m-1}(t)} = b(t) \quad (A2)$$

Introducing a functional change via: $u(t) = \frac{1}{y^{m-1}(t)}$, we obtain:

$$u'(t) = \frac{(1-m)y'(t)}{y^m(t)} \quad (A3)$$

Bernoulli equation is thus equivalent to a first-order in $u$ ODE: $\frac{1}{1-m}u'(t) + a(t)u(t) = b(t)$ with solution $y = u^{1/(1-m)}$.

## APPENDIX B: TAYLOR SERIES FOR FUNCTIONS OF SEVERAL VARIABLES

Given a function $f(\mathbf{r})$ with $\mathbf{r} = (x_1, x_2...x_n) \in \mathcal{R}^n$, a small displacement $\epsilon \in \mathcal{R}^n$ yields [15]:

$$f(\mathbf{r} + \boldsymbol{\epsilon}) = e^{(\boldsymbol{\epsilon} \cdot \nabla)} f(\mathbf{r}) \quad (B1)$$

where $\nabla_i = \frac{\partial}{\partial x_i}, i = 1...n$ that is:

$$f(\mathbf{r} + \boldsymbol{\epsilon}) = \sum_{i=0}^{\infty} \frac{1}{i!} (\boldsymbol{\epsilon} \cdot \nabla)^i f(\mathbf{r}) \quad (B2)$$

Applying this to $f(\mathbf{r}) = f(t, y)$ with $n = 2$ and $x_1 = t, x_2 = y$, we get $\boldsymbol{\epsilon} \cdot \nabla = \epsilon_1 \frac{\partial}{\partial t} + \epsilon_2 \frac{\partial}{\partial y}$ and consequently:

$$f(\mathbf{r} + \boldsymbol{\epsilon}) = \sum_{i=0}^{\infty} \frac{1}{i!} \left( \epsilon_1 \frac{\partial}{\partial t} + \epsilon_2 \frac{\partial}{\partial y} \right)^i f(\mathbf{r}) \quad (B3)$$

This results into the following expressions:

$$\begin{aligned} k_1 &= f(t, y) \\ k_2 &= f(t + c_2 h, y + h a_{21} k_1) \end{aligned} \quad (B4)$$
Thus $\epsilon_1 = c_2h$ and $\epsilon_2 = ha_{21}k_1$ and the operator $(\epsilon_1 \frac{\partial}{\partial t} + \epsilon_2 \frac{\partial}{\partial y})$ becomes: $h(c_2 \frac{\partial}{\partial t} + a_{21}k_1 \frac{\partial}{\partial y})$.

Since the operator is $\propto h$ it suffices to Taylor expand $k_2$ to $O(h^2)$, in order to obtain $y(t+h) = y(t)+h(b_1k_1+b_2k_2)$ to $O(h^3)$:

$$k_2 = f(t, y) + h(c_2 f_t + a_{21} f_y) + O(h^2) \quad (B5)$$

Substituting $k_2$ produces successively:

$$\begin{aligned}
y(t+h) &= y(t) + h(b_1k_1 + b_2k_2) \\
&= y(t) + hb_1f + hb_2[f + hc_2f_t + ha_{21}f_yf + O(h^2)] \\
&= y(t) + h(b_1 + b_2)f + h^2b_2(c_2f_t + a_{21}f_yf) + O(h^3)
\end{aligned} \quad (B6)$$

recovering the result announced in eq. 8.

[1] D.E. Potter, *Computational physics*, John Wiley & Sons, London, New York [1973].

[2] J.D. Lambert, *Numerical Methods for Ordinary Differential Systems*, John Wiley & Sons, New York, (1991).

[3] S. D. Conte, C. de Boor, *Elementary Numerical Analysis*, McGraw-Hill, New York, (1980).

[4] E. Hairer, G. Wanner, *Solving Ordinary Differential Equations II*, Springer-Verlag, New York, (1991).

[5] J.C. Butcher, *The Numerical Analysis of Ordinary Differential Equations: Runge-Kutta and General Linear Methods*, John Wiley & Sons, Chichester, (1987).

[6] T.E. Hull and R.L. Johnston, *Optimum Runge-Kutta methods*, Math. Comp., **18**, 306 (1964).

[7] L. Collatz, *The numerical treatment of differential equations*, third edition translated from a supplemented version of the second German edition by P. G. Williams, New York, Springer-Verlag, (1966).

[8] W. H. Press, W. T. Vetterling, S. A. Teukolsky and B. P. Flannery, *Numerical Recipes in C: The Art of Scientific Computing* Third Edition, Cambridge University Press, New-York (2007).

[9] E. Fehlberg, *Classical fifth-, sixth-, seventh-, and eighth-order Runge-Kutta formulas with stepsize control*, NASA Technical Report R-287 (1968).

[10] J.R. Cash and A.H. Karp, *ACM Transactions on Mathematical Software*, **16**, 201 (1990).

[11] L. F. Shampine, M. K. Gordon and J. Wisniewski, *A Variable order Runge-Kutta codes* In Computational Techniques for Ordinary Differential Equations Conference (Univ. of Manchester, 1978), I Gladwell and D. K. Sayers, Eds. Academic Press, London, p. 83-101 (1980).

[12] J. R. Rice, *Numerical Methods, Software and Analysis*, Second Edition, Academic Press, New-York (1993).

[13] J. Stoer and R. Bulirsch, *Introduction to Numerical Analysis*, Second Edition, Springer-Verlag, (New York, 1992).

[14] J. H. Mahaffy, Pennsylvania State University lecture notes (2018).

[15] L. D. Landau and E. M. Lifshitz, *Quantum Mechanics*, Vol. 3, Third Edition, Pergamon, Oxford (1977). See also G. F. Arfken and H. J. Weber, *Mathematical methods for physicists*, Sixth Edition, Elsevier Academic Press (London, 2005).