# Mixed precision implicit numerical schemes for systems of ordinary differential equations

Mouhamad Al Sayed Ali, Samuel Bernard, Arsène Marzorati, Jonathan

Rouzaud-Cornabas

► To cite this version:

Mouhamad Al Sayed Ali, Samuel Bernard, Arsène Marzorati, Jonathan Rouzaud-Cornabas. Mixed precision implicit numerical schemes for systems of ordinary differential equations. Numerical Algorithms, 2025, 10.1007/s11075-025-02256-w. hal-05385014

HAL Id: hal-05385014

https://inria.hal.science/hal-05385014v1

Submitted on 22 Jan 2026

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 CC BY 4.0 - Attribution - International License
# Mixed precision implicit numerical schemes for systems of ordinary differential equations

Mouhamad Al Sayed Ali * Samuel Bernard† Arsène Marzorati ‡

Jonathan Rouzaud-Cornabas§

September 1, 2025

## Abstract

Ordinary differential equations (ODEs) are widely used to model complex systems in biology, which result from the interactions of a large number of cells or organisms. This can lead to a substantial system size. These complex interactions can quickly alter their behavior, and some biological systems are stiff when represented as ODEs [61]. Therefore, these stiff ODEs could benefit from implicit numerical schemes. However, each iteration of these schemes involves solving a large nonlinear system, typically using the Newton method. To guarantee the global convergence of this method, we use line search (LS) and trust region (TR) algorithms. In this paper, we introduce a new approach to accelerate the computation of the implicit schemes by using mixed precision arithmetic, combining float and double precision, within the LS and TR algorithms, as well as in the Newton method. This approach aims at balancing the performance of lower precision arithmetic with the accuracy of higher precision arithmetic. We give theoretical results that show the efficiency of our approach. Numerical experiments show that our approach, running in either sequential or in parallel with MPI, is up to twice as fast as the double precision approach with the same level of accuracy. These experiments also show that increasing the size of the ODEs does not impact the quality of our mixed precision solution.

## 1 Introduction

In biology, and more generally, in all fields of scientific computing, including AI, there is a need for more detailed models that are more realistic, while also aiming to simulate larger systems [61]. However, the growth in computing resources and their performance is insufficient to meet the exponentially increasing demands for computational power [55]. On modern architectures, the performance of single precision (usually occupying 32 bits in memory) arithmetic operations is often at least twice as fast as the performance of double precision (64 bits) operations ([1, 42]). Many numerical algorithms are stable but they are costly in execution time when running in double precision. However, performing these algorithms in a lower precision (e.g. single precision arithmetic) can be an option, but it comes at the cost of reduced accuracy. In such case, lowering the arithmetic precision, in some parts of the algorithm, could speed up computations without compromising accuracy. Mixed precision algorithms, which combine low and high arithmetic precisions, could therefore be used to increase performances while maintaining high accuracy. Mixed precision algorithms have become popular in numerical linear algebra [1, 2, 68], iterative schemes [62], machine learning [27, 54], climate and weather model simulation [11, 57, 58, 56, 66, 31] and for numerical integration [22, 38, 39]. A variety of mixed precision methods in numerical linear algebra, both direct and iterative, for problems including matrix factorization, linear systems, have been developed in [42].

<footnote>*Inria Lyon, 69100, Villeurbanne, France. (mouhamad.alsayedali@ul.edu.lb)</footnote>

<footnote>†Université Claude Bernard Lyon 1, CNRS, Ecole Centrale de Lyon, INSA Lyon, Université Jean Monnet, ICJ UMR5208, Inria, Villeurbanne. (samuel.bernard@inria.fr)</footnote>

<footnote>‡Inria Lyon, 69100, Villeurbanne, France. (arsene.marzorati@inria.fr)</footnote>

<footnote>§CITI, INSA Lyon, CNRS, Inria, LIRIS UMR5205, Université Claude Bernard Lyon 1, ECL, Université Lumière Lyon 2, Lyon, France. (jonathan.rouzaud-cornabas@inria.fr)</footnote>
In this paper, we explore the use of mixed precision arithmetic in implicit solvers for scaling ODE
system sizes. Computing these schemes in a double precision require too much time, with respect to
the single precision scheme, but with very good accuracy. In order to accelerate the computation of
the double precision implicit schemes, with little loss of accuracy, we selectively compute some parts of
these schemes in single precision. Although, our approach can be applied for small systems of differential
equations, we are mainly interested in using mixed precision for large systems of ordinary differential
equations. The use of mixed precision for solving large systems of ODEs by explicit numerical schemes
has been studied in [3].

A few studies have investigated the use of mixed precision in the numerical solution of systems of
ordinary differential equations by implicit schemes. These schemes require at each integration step the
computation of the solution of a large nonlinear system. This system is solved by the Newton method. In
[49], the use of reduced arithmetic precision to solve a large nonlinear system by the Newton method is
investigated. Results in [49] show that using single precision in the numerical solution of the linear systems
arising in this method, has very little negative effect on its nonlinear convergence rate. However, this
approach can only be used for solving one nonlinear system and it becomes expensive when solving a great
number of large nonlinear systems. Moreover, the author has not investigated the global convergence of
the Newton method. In [38], a general theoretical additive framework for designing mixed precision for
the implicit Runge-Kutta schemes is provided, and the use of this framework to derive order conditions
for such schemes. The nonlinear systems in these schemes are solved by an implicit solver. Theoretical
results are given that show using mixed precision leverage low precision computation of the implicit solver
while retaining high precision in the overall scheme. However, this approach is only suitable when using
implicit Runge-Kutta schemes for solving small systems of differential equations. Moreover, there is no
discussion on how the nonlinear systems have been solved. In [22], the work presented in [38] is improved
by adding some corrected stages performed at high precision that are cheap explicit computations.
The results presented confirm that these new mixed precision Runge-Kutta schemes provide significant
runtime (and energy) savings per level of error. The nonlinear systems are solved by the Newton method
without investigating the global convergence of this method, and the numerical experiments have been
performed on small systems of ODEs. In [13], mixed precision arithmetic is used in the exponential
Rosenbrock-Euler method for solving ordinary differential equations. This method requires, at each
integration step, the computation of a matrix exponential-vector product. This product is computed
by the Arnoldi's method where the matrix-vector multiplication is performed in single precision. They
numerically showed the efficiency of their approach on a small system of ODE.

No prior research has investigated the use of mixed precision to accelerate the computation of implicit
schemes for solving systems of ordinary differential equations (ODEs). In this paper, we propose an
approach aimed at accelerating the computation of the solution of ODEs, using a great variety of implicit
schemes, by reducing the arithmetic precision in some evaluations of the right-hand side of the ODE. This
reduction in precision enhances speed and reduces memory usage while minimizing the loss of accuracy.

This paper is structured as follows. Section 2 introduces a large class of implicit numerical schemes for solving the ODE, along with the Newton method used to solve the nonlinear systems arising in these schemes. Then, section 2.1 presents our approach that incorporate mixed precision arithmetic in the implicit schemes. In section 2.1.1, we describe a mixed precision line search algorithm and explain its use, alongside mixed precision GMRES, to compute an accurate initial guess for each nonlinear system encountered in the implicit scheme. Section 2.1.2 introduces a trust region algorithm, which also uses mixed precision, to obtain a good initial guess for each nonlinear system in the scheme. Section 2.2 presents two benchmark models, a neural field model to describe the spatiotemporal dynamics of neural interactions in the cortex as a continuous space [12, 67], and a multiscale model for the regulation of the cell cycle by the circadian clock [32]. Section 2.3 presents the implicit schemes that we will use in our numerical experiments. The algorithmic aspect of our mixed implicit schemes is given in section 2.4. The computational environment and the parameters used in our numerical experiments are given in section 3.1. In Section 4.1, we give theoretical results that show the efficiency of our mixed precision approach. Section 4.2 gives numerical results that show that our mixed precision approach applied, on the benchmark models presented in Section 2.2, is up to twice as fast as than double precision approach, while maintaining high accuracy. However, the implicit schemes running in single precision are up to three times faster than those in double precision, but they either fail to provide sufficient accuracy or diverge. Moreover, as the timestep converges to zero, our mixed solution converges to the exact solution while the single solution diverges. We also show that the size of the ODEs does not impact the quality
of our mixed solution.

## 2 Materials and methods

sec:Imp

Consider the system of ODEs

$$ \dot{y}(t) = f(t, y(t)), \quad t_0 \le t \le T, \quad y(t_0) = y^{(0)}, \qquad (1) $$

eqo1

where $y(t) \in \mathbb{R}^n$ and $n$ is an integer.

A class of implicit numerical schemes for solving (1) is given by

$$ G_i(y_i) = y_i - a_i - \beta h f(t_i, y_i) = 0, \quad i = q + 1, \dots, N, \qquad (2) $$

S1

where $y_0 = y^{(0)}$ and $y_1, \dots, y_q$ are assumed to be known, $q \ll N$, $y_i$ is an approximation to $y(t_i)$ with $t_i = t_0 + ih$, $h = \frac{T-t_0}{N}$, $\beta$ is a scalar, and $a_i = F(y_{i-q}, \dots, y_{i-1}, f)$ is a vector that depends on the solutions computed at the previous iterations and on the function $f$.

We assume that the numerical scheme (2) is of order $p$, such that $\max_{0 \le i \le N} \|y_i - y(t_i)\| = \mathcal{O}(h^p)$, where $\|\cdot\|$ denotes the Euclidean norm. We also assume that, for all vectors $z_k, w_k \in \mathbb{R}^n$, $k = 1, \dots, q$, the function $F$ satisfies:

$$ \|F(z_{i-q}, \dots, z_{i-1}, f) - F(w_{i-q}, \dots, w_{i-1}, f)\| \le \sum_{k=1}^{q} |\alpha_k| \|z_{i-k} - w_{i-k}\| + h \sum_{k=1}^{q} |\beta_k| \|f(t_{i-k}, z_{i-k}) - f(t_{i-k}, w_{i-k})\| \qquad (3) $$

F_Prop

for some constants $\alpha_k$ and $\beta_k$, $k = 1, \dots, q$.

Most of implicit numerical schemes satisfy the previous inequality. For example, when $a_i$ is given by the following form:

$$ a_i = \sum_{k=1}^{q} \alpha_k y_{i-k} + h \sum_{k=1}^{q} \beta_k f(t_{i-k}, y_{i-k}), \quad 1 \le k \le q. \qquad (4) $$

vect_ai

Most standard implicit schemes such as implicit Euler, Adams-Moulton and BDF schemes can be written as in (2), and also satisfy (4), as described in [10, 40, 41, 47].

We would like to mention that our approach is also applied to the implicit Runge-Kutta schemes that cannot be written as in (2) (see section 2.3).

Implicit numerical schemes have good stability properties, but require the solution of a nonlinear system (2) at each integration step. The Newton method can be used to solve (2), but to keep the reduce its computational cost, the inexact Newton method is usually preferred [29, 33, 45, 46, 48]. It is given by the following process:

Choose $y_i^{(0)}$ as initial guess for $y_i$

For $k = 0, 1, \dots$, until convergence do:

$$ \text{Solve } G_i'(y_i^{(k)}) s_i^{(k)} = -G_i(y_i^{(k)}) \text{ inexact} \qquad (5) $$

$$ \text{Set } y_i^{(k+1)} = y_i^{(k)} + \tilde{s}_i^{(k)} \qquad (6) $$

where $G_i'(y_i^{(k)})$ is the Jacobian matrix of $G_i$ at $y_i^{(k)}$, and $\tilde{s}_i^{(k)}$ is an approximate solution to (5). However, each iteration of (5)-(6) requires the computation of this Jacobian matrix and also the solution of the linear system (5). In the following section, and in order to reduce the computational cost, we modify the equation (2) by reducing the precision in the evaluation of the function $G_i$.

## 2.1 Mixed precision implicit numerical schemes

sec:Theory

In this section, we show how mixed precision can be used to accelerate the computation of the implicit numerical scheme (2). First, we define the rounding mapping $\text{fl}(y)$ the evaluation of $y$ to a tolerance $\varepsilon$ and is assumed to satisfy the error model $\text{fl}(y) = (I + \delta)y$ where $\|\delta\| \le \varepsilon$, with $\varepsilon$ depends on the precision used ([13, 42]). Our mixed scheme, derived from (2), is then given by

$$ \tilde{G}_i(\tilde{y}_i) = 0, \quad i = q + 1, \dots, N, \qquad (7) $$

S1_L
where $\tilde{G}_i(x) = x - \tilde{a}_i - \beta hf^{(L)}(t_i, x)$, $\tilde{a}_i = F(\tilde{y}_{i-q}, \dots, \tilde{y}_{i-1}, f)$ is performed in a higher precision and obtained by replacing the vectors $y_k$ in $a_i$ by $\tilde{y}_k$. The function $f^{(L)}(t, x) := \text{fl}(f(\text{fl}(t), \text{fl}(x)))$ is defined by performing all the arithmetic operations in the function $f$ in a lower precision. If all operations in (2) are performed in a lower precision, then we obtain the following lower method

$$ \text{fl}(y_i - a_i - \beta hf(t_i, y_i)) = 0, \quad i = q + 1, \dots, N. \qquad (8) $$

The inexact Newton method for (7), is given by

Compute $y_i^{(0)}$ as initial guess for $\tilde{y}_i$

For $k = 0, 1, \dots$, until convergence do:

$$ \text{Solve } \tilde{G}'_i(y_i^{(k)})s_i^{(k)} = -\tilde{G}_i(y_i^{(k)}) \text{ inexact} \qquad (9) $$

$$ \text{Set } y_i^{(k+1)} = y_i^{(k)} + \tilde{s}_i^{(k)}, \qquad (10) $$

where $\tilde{s}_i^{(k)}$ is an approximate solution to (9). Krylov subspace methods [59], for example GMRES [60], can be used to find $\tilde{s}_i^{(k)}$ in the affine subspace $\hat{s}_i^{(k)} + K_i^{(k)}$ where $\hat{s}_i^{(k)}$ is an initial guess for the exact solution and $K_i^{(k)}$ is a Krylov subspace constructed with the matrix $\tilde{G}'_i(y_i^{(k)})$ and the initial residual vector $\tilde{G}'_i(y_i^{(k)})\hat{s}_i^{(k)} + \tilde{G}_i(y_i^{(k)})$. One generally asks that $\tilde{s}_i^{(k)}$ satisfies

$$ \| \tilde{G}'_i(y_i^{(k)}) \tilde{s}_i^{(k)} + \tilde{G}_i(y_i^{(k)}) \| \le \epsilon_i^{(k)} \| \tilde{G}_i(y_i^{(k)}) \|, \qquad (11) $$

with some tolerance threshold $\epsilon_i^{(k)}$. The results in [29] indicate how $\epsilon_i^{(k)}$ must be chosen to ensure local convergence of the inexact Newton method (9)-(10).

Krylov subspace methods have been used, for example, in [29, 16, 19, 20] in the context of Newton's method and in [35, 17, 18, 34, 21, 43, 53] in the context of ODEs. These methods have the advantage of requiring only the multiplication of the Jacobian matrix by a vector, the Jacobian matrix needs not be formed explicitly. A forward finite difference scheme is used to approximate the Jacobian-vector product. Also, good convergence can be expected, provided an efficient preconditioner and/or a good initial guess $\hat{s}_i^{(k)}$ is available ([20, 9]).

If the Newton method converges, it computes an approximation $\tilde{y}_i = y_i^{(k_0)}$ that satisfies

$$ \| \tilde{G}_i(\tilde{y}_i) \| \le \epsilon_{iN}, \qquad (12) $$

where $\epsilon_{iN}$ is a given tolerance threshold.

Theoretical results given by [16, 65] show the local convergence of the inexact Newton method by limiting the residual error in terms of the arithmetic precision used in the Krylov solver.

Like Newton method, the inexact Newton method (9)-(10) converges rapidly provided the initial guess $y_i^{(0)}$ is close enough to $\tilde{y}_i$. Traditionally, $y_i^{(0)}$ is found with the help of a predictor, usually an explicit scheme applied to (1). The drawback is that often the method (9)-(10) does not converge with such a predictor (see [5, 7]). For this reason, the line search backtracking and trust region algorithms can be used to improve global convergence [19, 20, 33] from any initial guess $y_i^{(0)}$. These algorithms generate a sequence of vectors that converge under some conditions to a local minimum or a saddle point of the function $h_i(y) = \frac{1}{2}\|\tilde{G}_i(y)\|^2$. Each vector is obtained from the preceding one by adding a descent direction on which the effectiveness of these algorithms strongly depends.

### 2.1.1 Mixed precision line search algorithm (MP-LS)

Recall that the initial guess $y_i^{(0)}$ has to be chosen to guarantee the convergence of the sequence $(y_i^{(k)})$ to $\tilde{y}_i$. To this end, line search (LS) algorithm can be used. This algorithm computes a sequence of vectors that converges, under some conditions, to a local minimum or a saddle point of the function

$$ h_i(y) = \frac{1}{2} \| \tilde{G}_i(y) \|^2. \qquad (13) $$

We refer to [30] for more details on this algorithm.
The sequence of vectors, $u_i^{(k)}$, computed by the line search algorithm is given by

$$u_i^{(k+1)} = u_i^{(k)} + \lambda_k p_i^{(k)}, \quad k \ge 0, \lambda_k > 0, \qquad (14)$$

where $u_i^{(0)}$ is an initial guess and $p_i^{(k)}$ is a descent direction of $\tilde{G}_i$ at $u_i^{(k)}$ (that is, $\nabla h_i(u_i^{(k)})^T p_i^{(k)} < 0$) and the scalar $\lambda_k$ is chosen so that $u_i^{(k+1)}$ satisfies the Wolf condition

$$h_i(u_i^{(k+1)}) \le h_i(u_i^{(k)}) + \alpha \lambda_k \nabla h_i(u_i^{(k)})^T p_i^{(k)}, \qquad (15)$$

and $\alpha \in (0, 1/2)$ is a parameter typically set to $10^{-4}$.

The computation of $\lambda_k$ uses the backtracking method, which starts with $\lambda_k = 1$ and repeatedly reduces it until an acceptable iterate $u_i^{(k+1)}$ satisfying (15) is found. The descent direction

$$p_i^{(k)} = - \{ \tilde{G}'_i(u_i^{(k)}) \}^{-1} \tilde{G}_i(u_i^{(k)}), \qquad (16)$$

which is referred to as the Newton direction, ensures that (15) is satisfied. However, this direction is expensive to compute. Therefore, using Krylov subspace methods, we look for an approximation $\tilde{p}_i^{(k)}$ to $p_i^{(k)}$ that satisfies

$$\| \tilde{G}'_i(u_i^{(k)}) \tilde{p}_i^{(k)} + \tilde{G}_i(u_i^{(k)}) \| \le \eta_i^{(k)} \| \tilde{G}_i(u_i^{(k)}) \|, \qquad (17)$$

with some tolerance threshold $\eta_i^{(k)} > 0$. Then we have

$$\begin{aligned} \nabla h_i(u_i^{(k)})^T \tilde{p}_i^{(k)} &= -\|\tilde{G}_i(u_i^{(k)})\|^2 + \tilde{G}_i(u_i^{(k)})^T (\tilde{G}'_i(u_i^{(k)}) \tilde{p}_i^{(k)} + \tilde{G}_i(u_i^{(k)})) \\ &\le (\eta_i^{(k)} - 1) \|\tilde{G}_i(u_i^{(k)})\|^2. \end{aligned}$$

In particular, if $\eta_i^{(k)} \ll 1$ then $\tilde{p}_i^{(k)}$ is a good descent direction.

Using this descent direction, by [30, Theorem 6.3.3], we get either $\nabla h_i(u_i^{(k)}) = 0$ for some $k \ge 0$, or $\lim_{k \to +\infty} \|\tilde{G}_i(u_i^{(k)})\| = 0$. If LS algorithm converges, it computes an initial guess $y_i^{(0)} := u_i^{(k_0)}$ that satisfies

$$\| \tilde{G}_i(y_i^{(0)}) \| \le \varepsilon_{LS},$$

where $\varepsilon_{LS}$ is a tolerance threshold.

In conclusion, when $\tilde{p}_i^{(k)}$ satisfies (17) with $0 < \eta_i^{(k)} \ll 1$, then $\tilde{p}_i^{(k)}$ is a descent direction that satisfies, or almost satisfies, (15) with $\lambda_k = 1$. In the latter case, the backtracking method will rapidly find a scalar $\lambda_k$ that satisfies (15) and the sequence $u_i^{(k)}$ converges to $\tilde{y}_i$. To accelerate the computation of $\tilde{p}_i^{(k)}$ we use mixed precision krylov subspace methods [37, 68]. The initial guess used in these methods can be either 0 or the solution of the Petrov-Galerkin projection on the subspace spanned by last r-vectors $f^{(L)}(t_{i-k}, \tilde{y}_{i-k})$, $k = 1, \dots, r$ where $r \ll n$, (see [5] for more details on the use of this projection). The linear systems (9) and (16) are solved using MPGMRES [68]. MPGMRES is a mixed version of the GMRES method that solves a large linear system of the form $Ax = b$, by computing the Hessenberg matrix with single precision. This method is given by the following algorithm:
**Algorithm 1** Mixed Precision Generalized Minimal RESidual (m)-MPGMRES(m) for solving Ax = b

**Input** $x_0$: initial guess, $\epsilon$: convergence criterion, $m$: restart frequency, maximum number of total inner iterations

1: $A^{(L)} = \text{fl}(A)$ \quad $\triangleright$ single precision of A

2: **repeat**

3: $r_0 = b - Ax_0$

4: $\beta_0 = \|r_0\|$

5: **If** $\beta_0/\|b\| \le \epsilon$ **Then return** $x_0$

6: $v_1 = r_0/\beta_0$

7: $\beta_0^{(L)} = \text{fl}(\beta_0)$

8: $v_1^{(L)} = \text{fl}(v_1)$

9: Compute $V_{m+1}^{(L)}$ and $\tilde{H}_m^{(L)}$ (Hessenberg matrix) by the m-step Arnoldi process with $A^{(L)}$ and $v_1^{(L)}$. $\triangleright$ using single precision arithmetic

10: Compute $y_m^{(L)}$ from $\beta_0^{(L)}$ and $\tilde{H}_m^{(L)}$. $\triangleright$ using single precision arithmetic

11: $z_m^{(L)} = V_m^{(L)} y_m^{(L)}$ \quad $\triangleright$ using single precision arithmetic

12: $z_m = \text{Std}(z_m^{(L)})$ \quad $\triangleright$ convert to double precision arithmetic

13: $x_m = x_0 + z_m$

14: $x_0 = x_m$

15: **EndIf**

16: **until** attain the maximum number of total inner iterations

**Output:** $x^{(M)} := x_0$

The algorithm 1 is used for solving the linear systems (5) and (16). If all operations in the algorithm 1 are performed in double precision (resp. single precision) then it will be denoted by DPGMRES(m) (resp. SPGMRES(m)).

Let $\gamma$ be the ratio of the number of repeating of the MPGMRES(m) over that in DPGMRES(m). $\gamma \ge 1$ is generally expected. The authors of [68] have shown that the speedup of MPGMRES(m) over DPGMRES(m) is estimated as

$$ \text{speedup} \approx \frac{2\rho + 6}{\gamma(\rho + 4)} $$

where $\rho = m/(N_{nz}/n)$, with $N_{nz}$ the number of nonzero elements in the matrix A. If $\gamma = 1$ then $1.5 < \frac{2\rho+6}{\rho+4} < 2$ which means no degradation in the convergence behavior, the maximum speedup is up to 2, which is reasonable considering the ratio of bytes of DPGMRES and SPGMRES.

The matrix of the linear systems (5) and (16) is of the form $A = \tilde{G}'_i(x)$. As its computation is too expensive, all matrix-vector multiplications $\tilde{G}'_i(x)v$ will be approximated by

$$ \frac{\tilde{G}_i(x + \sigma v) - \tilde{G}_i(x)}{\sigma} \qquad (18) $$

where $\sigma$ is a small scalar.

**Remark 1**

1. Let $\tilde{p}_{iM}^{(k)}$ be the solution of (17) by MPGMRES(m). Since the MPGMRES is always convergent for some m, there exists an integer $m_0$ such that the solution $\tilde{p}_{iM}^{(k)}$ of MPGMRES($m_0$) satisfies

$$ \|\tilde{G}'_i(u_i^{(k)})\tilde{p}_{iM}^{(k)} + \tilde{G}_i(u_i^{(k)})\| \le \eta_i^{(k)} \|\tilde{G}_i(u_i^{(k)})\|, \quad \eta_i^{(k)} \ll 1, $$

and therefore $\tilde{p}_{iM}^{(k)}$ is a descent direction.

2. The backtracking method, for computing $\lambda_k$ that satisfies (15), starts with $\lambda_k = 1$. If $p_i^{(k)}$ is the Newton direction, then the first iteration of LS algorithm is the Newton iteration. We could say that the goal of LS algorithm is to modify the Newton direction, by reducing its length.

In the next section, we present the trust region (TR) algorithm that also computes a good initial guess for the nonlinear system (7). Then, we show how mixed precision arithmetic can accelerate the computation of this algorithm, and therefore the computation of the initial guess $y_i^{(0)}$.
## 2.1.2 Mixed precision trust region algorithm (MP-TR)

The trust region algorithm generates a sequence of vectors $u_i^{(k)}$ that converges to either a local minimum or a saddle point of the function $h_i$. This sequence is given by $u_i^{(k+1)} = u_i^{(k)} + d_i^{(k)}$ where $d_i^{(k)}$ is the solution of

$$d_i^{(k)} = \underset{\|d\| \le \delta_i^{(k)}}{\operatorname{argmin}} \Psi_i^{(k)}(d), \qquad (19) \qquad \text{model}$$

where $\delta_i^{(k)} > 0$ and $\Psi_i^{(k)}$ is given by

$$\Psi_i^{(k)}(d) = \frac{1}{2} \| \tilde{G}_i(u_i^{(k)}) + \tilde{G}'_i(u_i^{(k)}) d \|^2.$$

The above function can be re-written as

$$\Psi_i^{(k)}(d) = h_i(u_i^{(k)}) + \nabla h_i(u_i^{(k)})^T d + \frac{1}{2} d^T B_i^{(k)} d \qquad (20)$$

where $B_i^{(k)} = \tilde{G}'_i(u_i^{(k)})^T \tilde{G}'_i(u_i^{(k)})$, $\nabla h_i(u_i^{(k)}) = \tilde{G}'_i(u_i^{(k)})^T \tilde{G}_i(u_i^{(k)})$, and $\tilde{G}'_i(u_i^{(k)})^T$ is the transpose of the matrix $\tilde{G}'_i(u_i^{(k)})$.

Since $B_i^{(k)}$ is symmetric and positive definite, the solution of (19) is given by

$$d_i^{(k)} = \begin{cases} -(B_i^{(k)})^{-1} \nabla h_i(u_i^{(k)}) & \text{if } \|(B_i^{(k)})^{-1} \nabla h_i(u_i^{(k)})\| \le \delta_i^{(k)}, \\ -(B_i^{(k)} + \mu_i^{(k)} I_n)^{-1} \nabla h_i(u_i^{(k)}) & \text{if } \|(B_i^{(k)})^{-1} \nabla h_i(u_i^{(k)})\| > \delta_i^{(k)} \end{cases}$$

where $\mu_i^{(k)} \ge 0$ solves

$$\|(B_i^{(k)} + \mu_i^{(k)} I_n)^{-1} \nabla h_i(u_i^{(k)})\| = \delta_i^{(k)}. \qquad (21) \qquad \text{Eq_Mu}$$

Note that $d_i^{(k)}$ is a descent direction of $\tilde{G}_i$ at $u_i^{(k)}$ and coincides with the Newton direction (16) when $\|(B_i^{(k)})^{-1} \nabla h_i(u_i^{(k)})\| \le \delta_i^{(k)}$.

To compute $\mu_i^{(k)}$ that approximately satisfies (21), we will use two methods: the locally constrained optimal hook step method and the dogleg step method (see [30]). The former finds $\mu_i^{(k)}$ such that $\|(B_i^{(k)} + \mu_i^{(k)} I_n)^{-1} \nabla h_i(u_i^{(k)})\| \cong \delta_i^{(k)}$ and takes $u_i^{(k+1)} = u_i^{(k)} - (B_i^{(k)} + \mu_i^{(k)} I_n)^{-1} \nabla h_i(u_i^{(k)})$. The latter makes a piecewise linear approximation to the curve $\mu \to u_i^{(k)} - (B_i^{(k)} + \mu I_n)^{-1} \nabla h_i(u_i^{(k)})$ and takes $u_i^{(k+1)}$ as the point on this approximation such that $\|u_i^{(k+1)} - u_i^{(k)}\| = \delta_i^{(k)}$. Note that both methods necessitate the solution of several linear systems.

In order to reduce the computational cost, we can proceed as in [5] and replace the full quadratic model (20) by the lower-dimensional one

$$\Phi_i^{(k)}(c) = h_i(u_i^{(k)}) + \nabla h_i(u_i^{(k)})^T V_i c + \frac{1}{2} c^T (V_i^T B_i^{(k)} V_i) c, \qquad (22) \qquad \text{tr2}$$

where $V_i$ is an $n \times p$ matrix, with $p \ll n$, whose columns form an orthonormal basis for $V_i$, where $V_i$ is the subspace spanned by last $r$-vectors $f^{(L)}(t_{i-k}, \tilde{y}_{i-k})$, $k = 0, \dots, r-1$ where $r \ll n$. Here, the same subspace $V_i$ is used across all iterations of the trust region algorithm.

However $V_i^T B_i^{(k)} V_i$ is always symmetric positive definite matrix but it is too expensive to compute. To reduce its computational cost, we replace it by $M_i^{(k)} := (L_i^{(k)})^T L_i^{(k)}$, with $L_i^{(k)} \approx \tilde{G}'_i(u_i^{(k)}) V_i$, such that the $l$-th column of $L_i^{(k)}$ is calculated, with single precision, by

$$\frac{\mathrm{fl}(\tilde{G}_i(\mathrm{fl}(u_i^{(k)} + \sigma v_i^{(l)}))) - \mathrm{fl}(\tilde{G}_i(\mathrm{fl}(u_i^{(k)})))}{\sigma}$$

where $v_i^{(l)}$ denotes the $l$-th column of $V_i$, and $\sigma$ is a small scalar. The matrix $L_i^{(k)}$ is then stored in double precision, and the matrix multiplication $(L_i^{(k)})^T L_i^{(k)}$ is also performed in double precision. The next iterate $u_i^{(k+1)}$ is then given by $u_i^{(k)} + V_i c_i^{(k)}$, where $c_i^{(k)}$ is the solution of

$$\min_{\|c\| \le \delta_i^{(k)}} \tilde{\Phi}_i^{(k)}(c), \qquad (23) \qquad \text{TRric}$$
with

$$ \tilde{\Phi}_i^{(k)}(c) = h_i(u_i^{(k)}) + \nabla h_i(u_i^{(k)})^T V_i c + \frac{1}{2} c^T M_i^{(k)} c. \quad (24) $$

phi_tilde

The condition for accepting $u_i^{(k+1)}$ is the one written in (15), namely

$$ h_i(u_i^{(k+1)}) \le h_i(u_i^{(k)}) + \alpha \nabla h_i(u_i^{(k)})^T (u_i^{(k+1)} - u_i^{(k)}), \quad (25) $$

Trust_1

where $\alpha \in (0, \frac{1}{2})$. If $u_i^{(k+1)}$ does not satisfy (25), we reduce $\delta_i^{(k)}$ and return to (23) to compute $c_i^{(k)}$ by the hook step or the dogleg step method. If $u_i^{(k+1)}$ satisfies (25), on must decide whether $\delta_i^{(k+1)}$ should be increased, decreased, or kept the same for the next step $k+1$, see [30]. The above algorithm will be denoted by the method used for solving (23), namely Hookstep or Doglegstep.

The use of mixed precision in the Trust Region (TR) algorithm is not new. In [25], the authors used several precisions (half, single, and double) to solve a minimization problem, of small dimension, using the TR algorithm. They demonstrated that mixed precision can considerably reduce the computational load compared to fixed-precision solvers. They provide practically a criterion to determine whether either we increase or decrease the precision level of the algorithm. However, their approach is limited to solving a single minimization problem of small dimension and is not well-suited for solving a large number of larger minimization problems. Additionally, their approach is not designed for solving ODEs.

Remark_2

**Remark 2** Our approach employs mixed precision at various stages of the implicit schemes. The first stage involves replacing $G_i$ with $\tilde{G}_i$. The second stage applies mixed precision linear solvers to solve the linear systems. Finally, the third stage involves replacing the constrained model, within the trust region algorithm, with other mixed precision model.

## 2.2 Benchmark models

In this section, we present two benchmark models that we will use to show the efficiency of using mixed precision for computing the implicit numerical schemes.

### 2.2.1 Benchmark model 1

For the first benchmark, we consider a mathematical model for the regulation of the cell cycle by the circadian clock [32]. This ODE is of the form (1), defined on $t \in [0, 120]$, where $y = [y^{(1)}, \dots, y^{(d)}]^T$, $y^{(i)} = [y_1^{(i)}, \dots, y_{10}^{(i)}]$, for $i = 1, \dots, d$, and $y^{(i)}$ satisfies the following system

$$ dy_1^{(i)}/dt = \frac{1}{\tau} \left( \nu_{1b}(y_7^{(i)} + \mathbf{\Psi}_1^{(i)})/(k_{1b}(1 + (y_3^{(i)}/k_{1i})^{p_0}) + y_7^{(i)} + \mathbf{\Psi}_1^{(i)}) - k_{1d}y_1^{(i)} \right) $$

$$ dy_2^{(i)}/dt = \frac{1}{\tau} \left( k_{2b}(y_1^{(i)})^q - k_{2d}y_2^{(i)} - k_{2t}y_2^{(i)} + k_{3t}y_3^{(i)} \right) $$

$$ dy_3^{(i)}/dt = \frac{1}{\tau} \left( k_{2t}y_2^{(i)} - k_{3t}y_3^{(i)} - k_{3d}y_3^{(i)} \right) $$

$$ dy_4^{(i)}/dt = \frac{1}{\tau} \left( \nu_{4b}(y_3^{r_0})^{(i)}/(k_{4b}^{r_0} + (y_3^{(i)})^{r_0}) - k_{4d}y_4^{(i)} \right) $$

$$ dy_5^{(i)}/dt = \frac{1}{\tau} \left( k_{5b}y_4^{(i)} - k_{5d}y_5^{(i)} - k_{5t}y_5^{(i)} + k_{6t}y_6^{(i)} \right) $$

$$ dy_6^{(i)}/dt = \frac{1}{\tau} \left( k_{5t}y_5^{(i)} - k_{6t}y_6^{(i)} - k_{6d}y_6^{(i)} + k_{7a}y_7^{(i)} - k_{6a}y_6^{(i)} \right) $$

$$ dy_7^{(i)}/dt = \frac{1}{\tau} \left( k_{6a}y_6^{(i)} - k_{7a}y_7^{(i)} - k_{7d}y_7^{(i)} \right) $$

$$ dy_8^{(i)}/dt = \lambda \left( (k_{impf} + k_{0mpf} \exp(-\eta d))k_{1mpf}^{n_0}/(k_{1mpf}^{n_0} + (y_8^{(i)})^{n_0} + s(y_{10}^{(i)})^{n_0})(1 - y_8^{(i)}) - d_{wee1}y_9^{(i)}y_8^{(i)} \right) $$

$$ dy_9^{(i)}/dt = \lambda \left( \frac{k_{actw}}{k_{actw} + d_{w1}}(c_w + C(y_7 - b_{bmal0}) + b_{bmal0}) + \dots \right. \\ \left. \left( \frac{k_{actw}}{k_{actw} + d_{w1}} - 1 \right)k_{inactw}(y_8^{(i)})^{n_0}y_9^{(i)}/(k_{1wee1}^{n_0} + (y_8^{(i)})^{n_0}) - d_{w2}y_9^{(i)} \right) $$

$$ dy_{10}^{(i)}/dt = \lambda (k_{act}(y_8^{(i)} - y_{10}^{(i)})) $$
These systems are coupled with each other via the variable

$$ \Psi_1^{(i)} = \frac{k_s}{d} \sum_{j=1}^{d} \arctan(y_2^{(j)} - y_2^{(i)}) + \frac{\pi}{2} k_s, \quad i = 1, \dots, d. $$

The constant parameters used in this model and the initial solution $y_0$ are provided in Appendix 1. In the numerical tests, we will consider $d = 10^4$, resulting in a ODE of size $n = 10^5$.

### 2.2.2 Benchmark model 2

For the second benchmark, we consider a neural field model that treats the cortex as a continuous space and describes the spatiotemporal dynamics of the neural activity given long range interactions [12, 67]. The neuronal membrane potential $V(x, t)$ at location $x$ and time $t$ follows the integral-differential equation

$$ \frac{\partial V}{\partial t}(x, t) = I(x, t) - V(x, t) + \int_{-1}^{1} K(|x - z|)S(V(z, t))dz, \quad x \in [-1, 1], t \in [0, 1], \quad (26) $$

where $I(x, t)$ is an external source of stimulus; $S(V)$ is the dependence between the firing rate of the neuron and the membrane potential and $K(|x - z|)$ is the connectivity between neurons at locations $x$ and $z$.

The integral term can be discretized. Let $x_i = -1 + i\Delta x, i = 0, \dots, d$, be a uniform grid of the interval $[-1, 1]$. Then

$$ \int_{-1}^{1} K(|x - z|)S(V(z, t))dz = \sum_{i=0}^{d-1} \int_{x_i}^{x_{i+1}} K(|x - z|)S(V(z, t))dz. $$

In each subinterval $[x_i, x_{i+1}]$, we introduce $k_\ell$ Gaussian nodes: $x_{i,s} = x_i + \frac{\Delta x}{2}(1 + \tau_s), s = 1, \dots, k_\ell$, where $\tau_s$ are the roots, sorted in an increasing order, of the $k_\ell$-th degree Legendre polynomial. Using the composite left rectangle method to evaluate the integral $\int_{x_i}^{x_{i+1}} K(|x - z|)S(V(z, t))dz$ with the points $\{x_{i,s}\}_{1 \le s \le k_\ell}$, and then evaluating (26) at these points, we get a system of ODEs of the form (1), where $n = d \times k_\ell, y(t) = [y_1(t), \dots, y_d(t)]^T$ with

$$ y_i(t) = [y_1^{(i)}(t), \dots, y_{k_\ell}^{(i)}(t)] \approx [V(x_{i,1}, t), \dots, V(x_{i,k_\ell}, t)] \in \mathbb{R}^{k_\ell}. $$

The parameters used in our numerical tests are $K(x) = e^{-x^2}$, $S(x) = I(x, t) = \tanh(x)$, and $k_\ell = 10$. The initial solution $y_0$ is given by $[y_1^{(0)}, \dots, y_n^{(0)}]^T$ where $y_i^{(0)} = \exp(6(i - 0.5n)/n), i = 1, \dots, n$. In the numerical tests, we use $n = 10^6$ ($d = 10^5$ since $k_\ell = 10$).

## 2.3 Implicit numerical schemes

The implicit numerical schemes used for solving the aforementioned benchmarks are the s-steps Backward Differential Formulas (BDFs), the Adams-Moulton of order 2 (AM2) and the Implicit Runge-Kutta (IRK).

The s-step (BDFs) with $s \le 6$ (see [10, 24, 44, 63]) are given by :

$$ (\text{BDF1}) : y_{i+1} - y_i = hf(t_{i+1}, y_{i+1}) \quad (\text{implicit Euler scheme}) $$

$$ (\text{BDF2}) : y_{i+2} - \frac{4}{3}y_{i+1} + \frac{1}{3}y_i = \frac{2}{3}hf(t_{i+2}, y_{i+2}) $$

$$ (\text{BDF3}) : y_{i+3} - \frac{18}{11}y_{i+2} + \frac{9}{11}y_{i+1} - \frac{2}{11}y_i = \frac{6}{11}hf(t_{i+3}, y_{i+3}) $$

$$ (\text{BDF4}) : y_{i+4} - \frac{48}{25}y_{i+3} + \frac{36}{25}y_{i+2} - \frac{16}{25}y_{i+1} + \frac{3}{25}y_i = \frac{12}{25}hf(t_{i+4}, y_{i+4}) $$

$$ (\text{BDF5}) : y_{i+5} - \frac{300}{137}y_{i+4} + \frac{300}{137}y_{i+3} - \frac{200}{137}y_{i+2} + \frac{75}{137}y_{i+1} - \frac{12}{137}y_i = \frac{60}{137}hf(t_{i+5}, y_{i+5}) $$

$$ (\text{BDF6}) : y_{i+6} - \frac{360}{147}y_{i+5} + \frac{450}{147}y_{i+4} - \frac{400}{147}y_{i+3} + \frac{225}{147}y_{i+2} - \frac{72}{147}y_{i+1} + \frac{10}{147}y_i = \frac{60}{147}hf(t_{i+6}, y_{i+6}). $$
The Adams-Moulton of order 2 (AM2) (Crank-Nicolson) scheme (see [40, 52]) is given by

$$y_{i+1} = y_i + \frac{h}{2} (f(t_{i+1}, y_{i+1}) + f(t_i, y_i)).$$

These numerical schemes are almost stable everywhere, and they are especially used for the solution of stiff differential equations. More details about these schemes and their stability can be found in [14, 64].

The s-stage implicit Runge-Kutta (IRK) scheme for solving (1) (see [15, 50, 51]) is given by

$$Y_{ik} - y_i - h \sum_{j=1}^{s} a_{kj} f(t_i + c_j h, Y_{ij}) = 0, \quad \text{for } k = 1, \dots, s, \qquad (27)$$

$$y_{i+1} = y_i + h \sum_{k=1}^{s} d_k f(t_i + c_k h, Y_{ik}), \qquad (28)$$

where $a_{kj}$, $c_j$, $d_k$ are the IRK coefficients, and $Y_{ik}$, $k = 1, \dots, s$, is an approximation of $y(t_i + c_k h)$.

Let $Y_i = [Y_{i1}^T, \dots, Y_{is}^T]^T$. Then, the equations (27) define a nonlinear system of the form

$$G_i(Y_i) = 0, \qquad (29)$$

where $G_i : \mathbb{R}^{n \times s} \to \mathbb{R}^{n \times s}$ is given by

$$G_i(X) = X - \mathbf{1}_s \otimes y_i - h(A_0 \otimes I_n)F_i(X).$$

The symbol $\otimes$ denotes the Kronecker product, $\mathbf{1}_s = (1, \dots, 1)^T \in \mathbb{R}^s$, $A_0 = (a_{ij}) \in \mathbb{R}^{s \times s}$, $F_i(X) = (f(t_i + c_1 h, X_1)^T, \dots, f(t_i + c_s h, X_s)^T)^T$, $I_n$ is the identity matrix of order $n$ and $X = [X_1^T, \dots, X_s^T]^T$.

Here, we consider IRK schemes of orders 4 ($s=2$) and 6 ($s=3$), denoted respectively by IRK4 and IRK6, that are defined with the following parameters (see [23, 24, 41]):

### Implicit Runge-Kutta of order 4 (IRK4)

* $s=2, c_1 = \frac{1}{2} - \frac{\sqrt{3}}{6}, c_2 = \frac{1}{2} + \frac{\sqrt{3}}{6}, d_1 = \frac{1}{2}, d_2 = \frac{1}{2}$.

* The matrix $A_0$ is given by

$$A_0 = \begin{pmatrix} \frac{1}{4} & \frac{1}{4} - \frac{\sqrt{3}}{6} \\ \frac{1}{4} + \frac{\sqrt{3}}{6} & \frac{1}{4} \end{pmatrix}.$$

### Implicit Runge-Kutta of order 6 (IRK6)

* $s=3, c_1 = \frac{1}{2} - \frac{\sqrt{15}}{10}, c_2 = \frac{1}{2}, c_3 = \frac{1}{2} + \frac{\sqrt{15}}{10}, d_1 = \frac{5}{18}, d_2 = \frac{4}{9}, d_3 = \frac{5}{18}$.

* The matrix $A_0$ is given by

$$A_0 = \begin{pmatrix} \frac{5}{36} & \frac{2}{9} - \frac{\sqrt{15}}{15} & \frac{5}{36} - \frac{\sqrt{15}}{30} \\ \frac{5}{36} + \frac{\sqrt{15}}{24} & \frac{2}{9} & \frac{5}{36} - \frac{\sqrt{15}}{24} \\ \frac{5}{36} + \frac{\sqrt{15}}{30} & \frac{2}{9} + \frac{\sqrt{15}}{15} & \frac{5}{36} \end{pmatrix}.$$

## 2.4 Algorithmic aspect of the computation of our mixed scheme

In this section, we write our Mixed Precision algorithm, denoted by MPAINS (Mixed Precision Accelerated Implicit Numerical Schemes), that computes the sequence $(\tilde{y}_i)$ defined by (7). The initial guess $y_i^{(0)}$ for (7) is obtained with MP-LS (part 3.1) or MP-TR (part 3.2). The approximate solution $\tilde{y}_i$ is obtained with the mixed Newton method (part 3.3) starting with $y_i^{(0)}$. The steps 2., (LS2.1), (N3) and 3.5 are optional. They are used only if we compute an initial guess for the MPGMRES method using the Petrov-Galerkin approximation. The parameters $r, \sigma, \epsilon_i^{(k)}, \epsilon_{LS}^{(1)}, \epsilon_{LS}^{(2)}, \epsilon_{IN}$, $kmax_{LS}$ and $kmax_{IN}$ should be fixed by the user.

1. Assume that for $k=1, \dots, q$, $\tilde{y}_k$ is either given or is computed with an $k$-step numerical scheme such that $\max_{1 \le k \le q} \|\tilde{y}_k - y_k\| = O(\epsilon_{IN})$.

Set $i = q + 1$.
2. (Optional) Let $R_i$ be the matrix formed by the last $k_0$ vectors $f^{(L)}(t_{i-k}, \tilde{y}_{i-k})$, $1 \le k \le k_0$, where $k_0 = \min(q+1, r)$. Orthonormalize the columns of $R_i$ in $V_i$.

3. Repeat until $i = N$

### 3.1 Line search for computing $y_i^{(0)}$:

(LS1) Set $k=0$, and compute an initial solution $u_i^{(0)}$ for LS (e.g. by an explicit scheme applied to (1))

(LS2) Repeat until $\|\tilde{G}_i(u_i^{(k)})\| \le \varepsilon_{LS}^{(1)}$ or $\|\tilde{G}_i(u_i^{(k+1)})\| - \|\tilde{G}_i(u_i^{(k)})\| \le \varepsilon_{LS}^{(2)}$ or $k > kmax_{LS}$

(LS2.1) (Optional) Compute an initial guess $\hat{p}_i^{(k)}$ for the linear system (17) by applying the Petrov-Galerkin process on $-u_i^{(k)} + \tilde{a}_i + \text{range}(V_i)$.

(LS2.2) If (LS2.1) is not used set $\hat{p}_i^{(k)} = 0$.

(LS2.3) If $\hat{p}_i^{(k)}$ satisfies (17) then $\tilde{p}_i^{(k)} = \hat{p}_i^{(k)}$. Otherwise, compute an approximation $\tilde{p}_i^{(k)}$ to (17) by MPGMRES($m$), using (18), and starting with $\hat{p}_i^{(k)}$.

(LS2.4) Set $p_i^{(k)} = \tilde{p}_i^{(k)}$ and compute the scalar $\lambda_k$ by the backtracking method such that (15) holds.

(LS2.5) Set $u_i^{(k+1)} = u_i^{(k)} + \lambda_k p_i^{(k)}$, $k = k+1$, $y_i^{(0)} = u_i^{(k)}$.

### 3.2 Trust region: replace the steps (LS2.1) to (LS2.5) above by

(TR1) Compute an approximation $c_i^{(k)}$ to (23) by Hookstep or Doglegstep method.

(TR2) Compute $u_i^{(k+1)} = u_i^{(k)} + V_i c_i^{(k)}$ satisfying (25). Set $k = k+1$, $y_i^{(0)} = u_i^{(k)}$.

### 3.3 Newton method for computing $\tilde{y}_i$:

(N1) Set $k=0$

(N2) Repeat until

$$ \|\tilde{G}_i(y_i^{(k)})\| \le \varepsilon_{IN} \text{ or } k > kmax_{IN} $$

(N3) (Optional) Compute an initial guess $\hat{s}_i^{(k)}$ to the linear systems (9) by applying the Petrov-Galerkin process on $-y_i^{(k)} + \tilde{a}_i + \text{range}(V_i)$.

(N4) If (N3) is not used set $\hat{s}_i^{(k)} = 0$.

(N5) If $\hat{s}_i^{(k)}$ satisfies (11) set $y_i^{(k+1)} = y_i^{(k)} + \hat{s}_i^{(k)}$, $k = k+1$ and go to (N2).

(N6) Else, compute an approximation $\hat{s}_i^{(k)}$ to (9), that satisfies (11), by MPGMRES($m$) using (18) and starting with $\hat{s}_i^{(k)}$.

(N7) Set $y_i^{(k+1)} = y_i^{(k)} + \tilde{s}_i^{(k)}$, $k = k+1$, go to (N2).

### 3.4 Set $\tilde{y}_i = y_i^{(k)}$.

### 3.5 (Optional) Let $k_0$ be the number of columns of $R_i$. If $k_0 < r$ then $R_{i+1} = [R_i, f^{(L)}(t_i, \tilde{y}_i)]$, else $R_{i+1} = [S_i, f^{(L)}(t_i, \tilde{y}_i)]$, where $S_i$ is the matrix formed by the last $r-1$ columns of $R_i$. Orthonormalize the columns of $R_{i+1}$ in $V_{i+1}$ using an updating QR factorization (see [36, p. 594]).

### 3.6 $i = i + 1$.

The final solution $\tilde{y}_N$ and the runtime of this algorithm will be denoted respectively by $y_M$ and $T_M$. If all arithmetic operations in this algorithm are performed in single (or double) precision, the above algorithm, the final solution, and the runtime will be denoted respectively by **SPAINS** "Single Precision Accelerated Implicit Numerical Schemes", $y_S$ and $T_S$ (or **DPAINS** "Double Precision Accelerated Implicit Numerical Schemes", $y_D$ and $T_D$).
# 3 Parameters for the numerical tests

In this section, we present the computational environment and the parameters that we use in our numerical tests.

## 3.1 Computational environment

Sec:Comp

The algorithms **SPAINS**, **MPAINS** and **DPAINS** were implemented in Fortran, in both sequential and parallel, and compiled on GCC 10.4.0 with the optimization option O3. All numerical tests were launched on the computer cluster Gros on Grid5000¹ (Intel Xeon Gold 5220, 2.20GHz). In the parallel mode, each evaluation of the function $f(t, y)$ was computed using parallelization MPI (version 4.1.5), and the numerical tests were launched with 8 hosts (144 processors) on Ethernet (25G) a network. For each numerical test, we ran 4 technical replicates of the simulations. Runtimes (in seconds) were averaged. The relative error with respect to the double solution was computed as the infinite norm of coefficient-wise relative errors $\|\frac{y_D - y_T}{y_D}\|_{\infty}$, where — is the pointwise division of two vectors, with $y_T \in \{y_S, y_M\}$.

## 3.2 Parameters for numerical results with stopping criterion the number of iterations

ed_Iteration

In the numerical tests, in section 4.2.1, and in order to facilitate a better comparison between the algorithms **SPAINS**, **MPAINS** and **DPAINS**, we fix the number of iterations $kmax_{LS}, m$ and $kmax_{IN}$ used in line-search, trust-region, GMRES and in the Newton method. The parameter $\sigma$ in (18) is set to $10^{-4}$. When employing LS algorithm, and since the Newton method is already integrated into the LS algorithm (see remark 1), we omit the steps (3.2) to (3.3) of the algorithm, setting $kmax_{IN} = -1$. Consequently, in this section, the Newton method will not be used with line search, because the initial guess provided by the line search is a sufficiently accurate approximation of the solution to the nonlinear system. For both the **BDFs** and the **AM2** schemes, the step (LS2) will be replaced with "while $k < kmax_{LS}$", where $kmax_{LS} = 5$, and the number of MPGMRES iterations will be set to $m = 10$. However, for the IRK scheme, we will use $kmax_{LS} = 10$ and $m = 5$ since Line Search has not converged with $kmax_{LS} = 5$. The restarted value used in GMRES is set to 1. When using TR algorithm, and for all numerical schemes, we use $kmax_{LS} = kmax_{IN} = m = 5$, and the number of vectors $r$ in the matrix $R_i$ is 2.

## 3.3 Parameters for numerical results with stopping criterion the tolerance

ed_Tolerance

In the previous section, we have fixed the number of iterations used in Line Search, in Trust region and in the Newton method. In section 4.2.2, we only set the tolerance $\epsilon$ to a fixed value. Then, we compute an approximation $\tilde{y}_i$ of each nonlinear system in (7) that satisfies

$$ \|\tilde{G}_i(\tilde{y}_i)\| \le \epsilon. $$

Therefore, the steps (LS2) and (N2) in the three algorithms **DPAINS**, **MPAINS**, and **SPAINS** are respectively replaced by

$$ \text{Repeat until } \|\tilde{G}_i(u_i^{(k)})\| \le \epsilon, $$

$$ \text{Repeat until } \|\tilde{G}_i(y_i^{(k)})\| \le \epsilon. $$

# 4 Results

Sec:Results

This section presents theoretical and numerical results that show the efficiency of our mixed algorithm **MPAINS**.

<footnote n="1">¹https://www.grid5000.fr/</footnote>
## 4.1 Theoretical error analysis

The results of the following theorem, supported by numerical experiments, demonstrate that our algorithm **MPAINS** outperforms **SPAINS** in terms of accuracy. Our algorithm enhances speed without compromising the quality of the numerical scheme (2).

**Theorem 1 (Error analysis of the methods (7) and (8))** Let us assume that the function $f(t, y)$ is Lipschitz with constant $c_0$, and that the time step satisfies $h \le \alpha_0 < \frac{1}{|\beta|c_0(1+\varepsilon)}$, where $\alpha_0$ is a constant independent of $\varepsilon$ (for e. g. $\alpha_0 = \frac{1}{2|\beta|c_0}$ since in fact $\varepsilon < 1$). Then, we have :

1. The error between $y_i$ and the solution $\tilde{y}_i$ of (7) satisfies

$$ \max_{0 \le i \le N} \|y_i - \tilde{y}_i\| = \mathcal{O}(\varepsilon(1 + h^p)). $$

2. The error between $y_i$ and the solution $\tilde{y}_i$ of (8) satisfies

$$ \max_{0 \le i \le N} \|y_i - \tilde{y}_i\| = \mathcal{O}\left(\frac{\varepsilon}{h}(1 + h^p)\right). $$

**Proof 1** 1. We have :

$$ \begin{aligned} \|y_i - \tilde{y}_i\| &= \|(a_i - \tilde{a}_i) + \beta h(f(t_i, y_i) - \text{fl}(f(\text{fl}(t_i), \text{fl}(\tilde{y}_i)))\| \\ &= \|(a_i - \tilde{a}_i) + \beta h(f(t_i, y_i) - f(\text{fl}(t_i), \text{fl}(\tilde{y}_i)))\| + \mathcal{O}(h\varepsilon) \end{aligned} $$

Since $f$ is Lipschitz, we get

$$ \begin{aligned} \|y_i - \tilde{y}_i\| &\le \|a_i - \tilde{a}_i\| + h|\beta|c_0 (|t_i - \text{fl}(t_i)| + \|y_i - \text{fl}(\tilde{y}_i)\|) + \mathcal{O}(h\varepsilon) \\ &\le \|a_i - \tilde{a}_i\| + h|\beta|c_0 (|t_i - \text{fl}(t_i)| + \|(y_i - \tilde{y}_i) + (\tilde{y}_i - \text{fl}(\tilde{y}_i))\|) + \mathcal{O}(h\varepsilon). \end{aligned} $$

Therefore, using the property of the model fl, we get

$$ \begin{aligned} \|y_i - \tilde{y}_i\| &\le \|a_i - \tilde{a}_i\| + h|\beta|c_0 (|t_i - \text{fl}(t_i)| + \|y_i - \tilde{y}_i\| + \varepsilon\|\tilde{y}_i\|) + \mathcal{O}(h\varepsilon) \\ &\le \|a_i - \tilde{a}_i\| + h|\beta|c_0\varepsilon (|t_i| + \|y_i\|) + h|\beta|c_0(1 + \varepsilon)\|y_i - \tilde{y}_i\| + \mathcal{O}(h\varepsilon). \end{aligned} $$

Then

$$ (1 - h|\beta|c_0(1 + \varepsilon)) \|y_i - \tilde{y}_i\| = \|a_i - \tilde{a}_i\| + \mathcal{O}\left(h\varepsilon\left(1 + \max_{0 \le i \le N} \|y_i\|\right)\right). $$

Using (3), and the fact that $\|y_i\| \le \|y_i - y(t_i)\| + \|y(t_i)\| = \mathcal{O}(h^p) + \|y\|_\infty$, we get

$$ (1 - h|\beta|c_0(1 + \varepsilon)) \|y_i - \tilde{y}_i\| = \sum_{k=1}^{q} |\alpha_k| \|y_{i-k} - \tilde{y}_{i-k}\| + h \sum_{k=1}^{q} |\beta_k| \|f(t_{i-k}, y_{i-k}) - f(\text{fl}(t_{i-k}), \text{fl}(\tilde{y}_{i-k}))\| + \mathcal{O}(h\varepsilon(1 + h^p)). $$

By a simple induction, the fact that $h \le \alpha_0 < \frac{1}{|\beta|c_0(1+\varepsilon)}$ and proceeding as before we deduce that

$$ \max_{0 \le i \le N} \|y_i - \tilde{y}_i\| = i \times \mathcal{O}(h\varepsilon(1 + h^p)) = \mathcal{O}(\varepsilon(1 + h^p)). $$

2. Let $\tilde{a}_i = \text{fl}(F(\text{fl}(\tilde{y}_{i-q}), \dots, \text{fl}(\tilde{y}_{i-1}), f))$, then proceeding in the same way of the previous part, we get

$$ (1 - h|\beta|c_0(1 + \varepsilon)) \|y_i - \tilde{y}_i\| = \sum_{k=1}^{q} |\alpha_k| \|y_{i-k} - \text{fl}(\tilde{y}_{i-k})\| + h \sum_{k=1}^{q} |\beta_k| \|f(t_{i-k}, y_{i-k}) - f(\text{fl}(t_{i-k}), \text{fl}(\tilde{y}_{i-k}))\| + \mathcal{O}(h\varepsilon(1 + h^p)). $$

Since $\|y_{i-k} - \text{fl}(\tilde{y}_{i-k})\| \le \|y_{i-k} - \tilde{y}_{i-k}\| + \varepsilon\|\tilde{y}_{i-k}\| = (1 + \varepsilon)\|y_{i-k} - \tilde{y}_{i-k}\| + \mathcal{O}(\varepsilon(1 + h^p))$, we get

$$ \max_{0 \le i \le N} \|y_i - \tilde{y}_i\| = \mathcal{O}\left(\frac{\varepsilon}{h}(1 + h^p)\right). \blacksquare $$

**Remark 3** The previous theorem shows that $y_M$ is a good approximation of $y_D$. It also shows that, when the timestep $h$ decreases to 0, the error between $y_M$ and $y_D$ decreases to $\varepsilon$, while the error between $y_S$ and $y_D$ increases and tends to $+\infty$.
## 4.2 Numerical Results

Sec:results

This section shows the numerical efficiency of our algorithm **MPAINS**. It presents the error and the run-time results when the algorithms **SPAINS**, **MPAINS**, and **DPAINS** are used to solve the benchmark models, given in section 2.2, with the implicit schemes defined in section 2.3.

### 4.2.1 Numerical results with the parameters given in Section 3.2

Sec::Para__1

This section presents the numerical results with the parameters given in section 3.2.

#### Numerical results for benchmark 1

This section presents the numerical results for benchmark 1. Figure 1 displays the outcomes of applying the algorithms **DPAINS**, **MPAINS**, and **SPAINS** to solve benchmark 1 with the implicit schemes BDF1, BDF4 and BDF6, using Line Search, Hookstep and Doglegstep, and with various timesteps. Numerical simulations, presented in Figure 1, show that the **MPAINS** algorithm outperforms **DPAINS**, being 1.4 to 1.6 times faster while it maintains good accuracy. They also show that the **SPAINS** algorithm outperforms **DPAINS**, being 2 to 2.5 times faster but with very poorer accuracy. For some cases, the numerical simulations show a noticeable increase of the error of the single solution. For example, when using BDF1 with any of the algorithms (Line Search, Doglegstep, or Hookstep) and a timestep of 0.012, the error between the single and the double solution is approximately $10^{-1}$, while the error between the mixed and the double solution is on the order of $10^{-12}$. We also deduce, from the regression line, that as the timestep $h$ decreases, the error of the mixed solution decreases while the error of the single solution either stagnates or increases. These trends are consistent with the theoretical results outlined in Theorem 1. The numerical results for the other implicit schemes are shown in Appendix 2.

#### Numerical results for the benchmark 2

Figure 2 presents the numerical results for the benchmark 2. These results show that the **MPAINS** algorithm is as fast as the **SPAINS** algorithm, where both algorithms are up to 1.5 times faster than **DPAINS**. However, the solution of our algorithm **MPAINS** is a good approximation of the double solution, while the solution of **SPAINS** algorithm either provides very poorer accuracy or diverges. The main different results with those of benchmark 1, is that as the timesteps decrease, the **SPAINS** algorithm diverges, while our algorithm **MPAINS** remains convergent. This divergence is influenced by factors such as the number of iterations used in the Newton method, the Line search and the Trust Region algorithm, as well as the number of vectors $r$ used in the subspace $V_i$. Increasing the number of iterations or the number of vectors $r$ in $V_i$ could potentially enable the Newton method to converge. The numerical results for the other implicit schemes are shown in Appendix 2.

### 4.2.2 Numerical results with the parameters given in Section 3.3

Sec::Para__2

This section presents the numerical results with the parameters given in section 3.3.

Figure 3 presents the numerical results when benchmark 1 is solved by BDF1 and BDF2 using the Line Search, with different values of timesteps and different values of $\varepsilon$. For each value of $\varepsilon$, all algorithms **SPAINS**, **MPAINS** and **DPAINS** approximately compute the same solution provided they converge. The numerical results, presented in this Figure, show that as $\varepsilon$ decreases to zero, the algorithm **SPAINS** fails to converge (black rhombuses), while our algorithm **MPAINS** remains convergent and its solution tends to the double solution. The results also show that error between our mixed solution and the double solution is too small with respect to the error between the single and the double solutions. Moreover, when the timestep $h$ is small, the runtime speedup of our mixed algorithm slightly decreases for bigger values of $\varepsilon$, as the nonlinear system $G_i(y)$ converges to a linear system and it becomes easy to solve by the Newton method because $||G_i(y) - (y - a_i)|| = O(h)$. However, as the timestep $h$ increases to 0.8, the runtime speedup of our algorithm **MPAINS** is 1.5 times faster than the algorithm **DPAINS** with very good accuracy. All algorithms diverge when using a tolerance $\varepsilon = 10^{-14}$ with $h = 0.8$, which means the initial guess computed by Line-Search is not good enough for the Newton method. Increasing the number of iterations used in Line-Search could lead to the convergence of the Newton method with this
<figure><img src="image_2.png" /><figcaption>Figure 1: Numerical results for benchmark 1 with Line search(top left), Doglegstep (top right) and Hookstep (bottom left). The x-axis shows the runtime in seconds while the y-axis shows the relative error(log scale) with respect to the double solution. The vertical solid lines indicate the runtime of the double solutions. The dotted lines represent the fitted linear regression curve with respect to the error of both the mixed and the single solutions.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Figure 2: Numerical results for benchmark 2. The filled black rhombuses in the Figure indicate that the Newton method has not converged when solving the nonlinear systems of BDF6.</figcaption></figure>

tolerance. In conclusion, our algorithm **MPAINS** converges whenever the **DPAINS** algorithm does, and with high accuracy, while the **SPAINS** algorithm only converges for large values of $\varepsilon$. Therefore, **MPAINS** provides a reliable and accurate approach for computing solutions to the ODEs.

### 4.2.3 Impact of the size of the ODEs on the numerical results

This section numerically shows how the size of the ODEs impacts the numerical results. Figure 4 presents the numerical results when the algorithms **SPAINS**, **MPAINS**, **DPAINS** are used to solve benchmark 1, with different system sizes. The implicit scheme used is BDF1 with the time step 0.12. The numerical simulations show that the size of the ODE has no effect on the runtime speedup or on the relative error. This indicates that as the size of the ODE increases, the runtime speedup and the relative error remain stable.

Figure 5 shows the numerical results for benchmark 2, with different system sizes, and with different timesteps. Figure 5-right reveals that for all size of ODEs, the algorithm **SPAINS** outperforms **MPAINS** in terms of speed. Figure 5-left shows that the relative error of both single and mixed solutions rises as the ODE size increases, and that **MPAINS** outperforms **SPAINS** in terms of accuracy with the time step 0.01, while the error of the solution of both algorithms is not small enough with the bigger timestep 0.1. When the timestep is reduced to 0.001, the **SPAINS** algorithm fails to converge(black rhombuses) for small sizes of the ODE, while our **MPAINS** algorithm remains convergent with very good speedup and accuracy. As the size increases, the **SPAINS** algorithm converges, but with a very poor level of accuracy.

### 4.2.4 Numerical results of sequential numerical methods

In order to show that using parallelization MPI has no effect on the quality of our mixed precision approach, we have launched numerical tests for benchmark 1 running in sequential on the computer cluster GRID5000 with a timestep of 0.12. Since the code is running sequentially, we have reduced the size of the benchmark 1 to $10^4$, as with a size of $10^5$ the sequential simulation takes too much runtime.
<figure><img src="image_4.png" /><figcaption>Figure 3: Benchmark 1 solved with BDF1 (Left) and BDF2 (right) using Line Search. The x-axis shows the runtime speedup while the y-axis shows the relative error.</figcaption></figure>

BEN1_Fixed

<figure><img src="image_5.png" /><figcaption>Figure 4: Benchmark 1 solved by BDF1. The x-axis shows the size of the benchmark. Relative error (left) and runtime speedup (right).</figcaption></figure>

Ben1_SCALING
<figure><img src="image_6.png" /><figcaption>Figure 5: Benchmark 2 solved by BDF1. Relative error (left) and runtime speedup (right).</figcaption></figure>

The numerical results show that our algorithm **MPAINS** is always up to 1.5 times faster than the **DPAINS** algorithm with very good accuracy, while **SPAINS** algorithm gives very poorer accuracy. We conclude that the parallelization does not affect the effectiveness of our mixed algorithm.

# 5 Conclusion

In this paper, we introduce an algorithm that leverages mixed precision arithmetic to solve (large) systems of ODEs through implicit numerical schemes. Each iteration of these schemes requires the solution of a large nonlinear system. This system is solved by the algorithms Line Search (LS), Trust Region (TR) and Newton. The cost of computation of these algorithms is reduced by performing some Right-Hand Sides (RHS) of the ODE with single precision. Line search and Trust region are employed to calculate a good initial guess for the Newton method, where their efficiency strongly depends on the choice of a descent direction. For LS, where the Newton's direction is commonly used, but requires solving a large linear system, we accelerate the process by employing a mixed precision GMRES method. For the TR algorithm, each iteration involves solving a constrained minimization problem, where we reduce computational costs by projecting the model onto a low-dimensional subspace and performing matrix computations when solving this problem in single precision.

Integrating mixed precision into these algorithms significantly enhance computational efficiency and reduce communication overhead, without sacrificing accuracy. We present theoretical results demonstrating that the error between the double and the mixed solutions is of order $\varepsilon$ where $\varepsilon$ is the machine epsilon of the precision used, while the error between the single and the double solutions is of order $\varepsilon/h$. This means that as $h$ tends to zero the single solution diverges. These theoretical results are validated by numerical experiments on two biological models. These experiments showed that our mixed precision algorithm, running in sequential or in parallel with MPI, achieves a 1.5–2 times speedup compared to its double precision counterpart, while maintaining near equivalent accuracy. The SPAINS algorithm only converges for big values of $\varepsilon$ and timesteps, and its solution (when it converges) is not close enough to the double solution. In this paper, single and double precision were tested but the results are expected to generalize to other precisions such as half and quadruple. The numerical experiments, in both sequential and parallel, show that increasing the size of the ODEs does not impact the quality or the runtime speedup of our mixed solution.

The mixed precision approach given in this paper can also be boosted with other techniques of acceleration such as preconditioning the linear systems, or computing a good initial guess for the linear systems arising in the Line Search algorithm and in the Newton method using the Petrov-Galerkin process (see [4, 5, 6, 7, 8] for more details about these techniques of acceleration).
# 6 Funding

This research received no specific grant from any funding agency in the public, commercial, or not-for-profit sectors.

# References

AA2021

[1] A. Abdelfattah and al., *A survey of numerical linear algebra methods utilizing mixed-precision arithmetic*, International Journal of High Performance Computing Applications, Volume 35, Issue 4, Jul (2021), 344-369.

ALSGK2019

[2] S. Abdulah, H. Ltaief, Y. Sun, M. G. Genton, D. E. Keyes, *Geostatistical modeling and prediction using mixed-precision tile cholesky factorization*, IEEE 26th International Conference on High Performance Computing, Data and Analytics (HiPC), (2019).

ABMB2024

[3] M. Al Sayed Ali, S. Bernard, A. Marzorati, J. Rouzaud-Cornabas, *Mixed precision approaches for explicit numerical methods for ordinary differential equations*, DRAFT, (2025).

AS2008

[4] M. Al Sayed Ali, *Accélération de schémas d'intégration temporelle pour les équations différentielles algébriques linéaires*, Comptes Rendus de l'Académie des Sciences de Paris, Ser. I 346 (2008), 229-232.

AS2009

[5] M. Al Sayed Ali, M. Sadkane, *Acceleration of implicit schemes for large systems of nonlinear ODEs*, Electron. T. Numer. Anal., 35 (2009), 104-117.

AS2012

[6] M. Al Sayed Ali, M. Sadkane, *Improved predictor schemes for large systems of linear ODEs*, Electron. T. Numer. Anal., 39 (2012), 253-270.

AS2019

[7] M. Al Sayed Ali, M. Sadkane, *Acceleration of implicit schemes for large systems of nonlinear differential-algebraic equations*, AIMS Mathematics, 5(1), (2019), 603-618.

AS2021

[8] M. Al Sayed Ali, M. Sadkane, *Acceleration of implicit schemes for large linear systems of differential-algebraic equations*, Journal of Computational and Applied Mathematics, Vol. 389, (2021), 113364.

MAS2025

[9] M. Al Sayed Ali, M. Sadkane, *Acceleration of implicit schemes for large systems of delay differential equations*, accepted for publication to Journal of Computational and Applied Mathematics, (2025).

AP1998

[10] U. M. Ascher, L. R. Petzold, *Computer methods for ordinary differential equations and differential-algebraic equations*, SIAM, Philadelphia, PA, USA, (1998).

ADP2022

[11] J. Ackmann, P. D. Dueben, T. Palmer, P. K. Smolarkiewicz, *Mixed-precision for linear solvers in global geophysical flows*, Journal of Advances in Modeling Earth Systems, Vol. 14, (2022).

A1977

[12] S.L. Amari, *Dynamics of pattern formation in lateral-inhibition type neural fields*, Biol. Cybernet. 27 (2) (1977), 77-87.

BRG2023

[13] C. J. Balos, S. Roberts, D. J. Gardner, *Leveraging Mixed Precision in Exponential Time Integration Methods*, IEEE High Performance Extreme Computing Conference, HPEC 2023, 1-8.

BCTBL2022

[14] L. Beuken, O. Cheffert, A. Tutueva, D. Butusov, V. Legat, *Numerical Stability and Performance of Semi-Explicit and Semi-Implicit Predictor-Corrector Methods*, Mathematics 2022, 10, 2015.

BCP89

[15] K.E. Brenan, S.L. Campbell, L.R. Petzold, *Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations*, SIAM Publications, Philadelphia, PA, 2nd edition, (1996).

B1987

[16] P.N. Brown, *A Local Convergence theory for Combined Inexact-Newton/Finite-Difference Projection Methods*, SIAM J. Numer. Anal., (April 1987), Vol.24, No.2., 407-434.

BH1986

[17] P.N. Brown, A. C. Hindmarsh, *Matrix-free methods for stiff of ODE's*, SIAM J. Numer. Anal., 23 (1986), 610-638.
[18] P.N. Brown, A. C. Hindmarsh, *Reduced storage methods in stiff ODE systems*, J. Appl. Math. Comput., 31 (1989), 40-91.

[19] P. N. Brown and Y. Saad, *Hybrid Krylov methods for nonlinear systems of equations*, SIAM J. Sci. Stat. Comput., 11, (1990), 450-481.

[20] P. N. Brown and Y. Saad, *Convergence theory of nonlinear Newton-Krylov algorithms*, SIAM J. Optimization, Vol. 4, No. 2, (1994), 297-330.

[21] M. Bünttner, B.A. Schmidt and R. Weiner, *W-methods with automatic partitioning by Krylov techniques for large stiff systems*, SIAM J. Numer. Anal., 32, (1995), 260-284.

[22] B. Burnett, S. Gottlieb, Z. J. Grant, A. R. H. Heryudono, *Performance evaluation of Mixed-Precision Runge-Kutta Methods*, IEEE High Performance Extreme Computing Conference (HPEC), (2021), 1-6.

[23] J. C. Butcher, *A history of Runge-Kutta methods*, Applied Numerical Mathematics, vol. 20, number 3, (1996), 247-260.

[24] J. C. Butcher, *Numerical schemes for ordinary differential equations*, John Wiley, (2003).

[25] R. J. Clancy, M. Menickelly, J. Hückelheim, P. Hovland, P. Nalluri, R. Gjini, *TROPHY: Trust Region Optimization Using a Precision Hierarchy*, Computational Science - ICCS 2022. ICCS 2022. Lecture Notes in Computer Science, vol 13350.

[26] M. Croci, G. R. Souza, *Mixed-precision explicit stabilized Runge-Kutta methods for single and multi-scale differential equations*, Journal of Computational Physics 464 (2022), 111139.

[27] D. Das, N. Mellempudi, D. Mudigere, D. Kalamkar, S. Avancha, K. Banerjee, S. Sridharan, K. Vaidyanathan, B. Kaul, E. Georganas, et al., *Mixed precision training of convolutional neural networks using integer operations*, ICLR (Poster), OpenReview.net, (2018).

[28] P.J. Davis, *Interpolation and approximation*, Blaisdell, New York, (1963).

[29] R.S. Dembo, S.C. Eisenstat, and T. Steihaug, *Inexact Newton methods*, SIAM J. Numer. Anal., 19, (1982), 400-408.

[30] J. E. Dennis and R. B. Schnabel, *Numerical schemes for unconstrained optimization and nonlinear equations*, Prentice-Hall, Englewood Cliffs, NJ, (1983).

[31] P.D. Duben, A. Subramanian, A. Dawson, T. Palmer, *A study of reduced numerical precision to make super parameterization more competitive using a hardware emulator in the OpenIFS model*, J. Adv. Model. Earth Syst. 9 (2017), 566-584.

[32] R. El cheikh, S. Bernard, N. El Khatib, *A multiscale modelling approach for the regulation of the cell cycle by the circadian clock*, Journal of Theoretical Biology 426 (2017), 117-125

[33] S.C. Eisenstat, H.F. Walker, *Globally convergent inexact Newton methods*, SIAM J. Optim., 4, (1994), 393-422.

[34] E. Gallopoulos, Y. Saad, *Efficient solution of parabolic equations by Krylov approximation methods*, SIAM J. Sci.Stat. Comput., 13 (1992), 1236-1264.

[35] C.W. Gear, Y.Saad, *Iterative Solution of Linear equations in ODE Codes*, SIAM J. Sci.Stat. Comput., (1983), Vol. 4, No. 4.

[36] G.H. Golub, C.F. Van Loan, *Matrix computations*, 2nd ed., The Johns Hopkins University Press, Baltimore, MD, (1989).

[37] S. Graillat, F. Jézéquel, T. Mary, and R. Molina, *Adaptive Precision Sparse Matrix-Vector Product and Its Application to Krylov Solvers*, SIAM Journal on Scientific Computing, Vol. 46, Iss. 1 (2024).
[38] Z. J. Grant, *Perturbed Runge-Kutta Methods for Mixed Precision Applications*, J. Sci. Comput. 92(1): 6, (2022).

[39] E. Hairer, R.I. McLachlan, A. Razakarivony, *Achieving Brouwer's law with implicit Runge-Kutta methods*, BIT Numer. Math. 48 (2008), 231-243.

[40] E. Hairer, S.P. Norsett, G. Wanner, *Solving ordinary differential equations I: Nonstiff problems*, 2nd ed., Springer-Verlag, Berlin, (1993).

[41] E. Hairer, G. Wanner, *Solving ordinary differential equations II. Stiff and differential-algebraic problems*, 2nd ed., Springer-Verlag, Berlin, (1996).

[42] N. J. Higham, T. Mary, *Mixed precision algorithms in numerical linear algebra*, Acta Numerica, vol. 32, (2022), 1-69.

[43] M. Hochbruck, C. Lubich, H. Selhofer, *Exponential Integrators for Large Systems of Differential Equations*, SIAM J. Sci. Comput. 19, (1998), 1552-1574.

[44] A. Iserles, *A First Course in the Numerical Analysis of Differential Equations*, Cambridge University Press, (1996).

[45] C. T. Kelley, *Iterative Methods for Linear and Nonlinear Equations*, SIAM, (1995).

[46] C. T. Kelley, *Solving Nonlinear Equations with Newton's Method*, Fundam. Algorithms 1, SIAM, Philadelphia, (2003).

[47] C. T. Kelley, *numerical schemes for nonlinear equations*, Acta Numer., 27 (2018), 207-287.

[48] C. T. Kelley, *Notebook for Solving Nonlinear Equations with Iterative Methods: Solvers and Examples in Julia*, IJulia Notebook, (2020).

[49] C. T. Kelley, *Newton's method in mixed precision*, SIAM Rev. 64, (2022), 191-211.

[50] L.O. Jay, *Preconditioning and parallel implementation of implicit Runge-Kutta methods*, Report on Computational Mathematics, No. 150, Department of Mathematics, University of Iowa, (2002).

[51] L.O. Jay, *Inexact Simplified Newton iterations for implicit Runge-Kutta Methods*, SIAM J. Numer. Anal., (2000), vol.38, No.4, 1369-1388.

[52] A. Quarteroni, R. Sacco, F. Saleri, *Matematica Numerica*, Springer Verlag, (2000).

[53] L.O. Jay, *Inexact simplified Newton iterations for implicit Runge-Kutta methods*, SIAM J. Numer. Anal. 38 (2001), 1369-1388.

[54] N. Mellempudi, S. Srinivasan, D. Das, B. Kaul, *Mixed precision training with 8-bit floating point*, arXiv preprint, 12334, (2019).

[55] NIVIDIA: https://www.nvidia.com/fr-fr/about-nvidia/ai-computing/

[56] E.A. Paxton, M. Chantry, M. Klower, L. Saffin, T. Palmer, *Climate modeling in low-precision: effects of both deterministic and stochastic rounding*, J. Climate, 35, (2022), 1215-1229.

[57] M. Klower, P. Duben, T. Palmer, *Number formats, error mitigation, and scope for 16-bit arithmetics in weather and climate modeling analyzed with a shallow water model*, J. Adv. Model. Earth Syst. 12 (2020).

[58] M. Klower, S. Hatfield, M. Croci, P.D. Duben, T.N. Palmer, *Fluid simulations accelerated with 16 bits: approaching 4x speedup on A64FX by squeezing ShallowWaters.jl into Float16*, J. Adv. Model. Earth Syst. 14 (2022).

[59] Y. Saad, *Iterative methods for sparse linear systems*, 2nd ed., SIAM, Philadelphia, PA, (2003).

[60] Y. Saad and M. H. Schultz, *GMRES: A Generalised minimal residual algorithm for solving non-symmetric linear systems*, SIAM J. Sci. Stat. Comp., 7 (1986), 856-869.
[61] P. Städter, Y. Schälte, L. Schmiester, J. Hasenauer, P. L. Stapor, *Benchmarking of numerical integration methods for ODE models of biological systems*, Scientific Reports, Vol. 11, 2696, (2021).

[62] R. Strzodka and D. Göddeke, *Mixed precision methods for convergent iterative schemes*, EDGE, vol. 6, pp. 23-24, 2006.

[63] E. Suli, D. Mayers, *An Introduction to Numerical Analysis*, (2003), Cambridge University Press.

[64] A. Tutueva, T. Karimov, D. Butusov, *Semi-Implicit and semi-Explicit Adams-Bashforth-Moulton Methods*, Mathematics (2020), 8(5), 780.

[65] F. Tisseur, *Newton's method in floating point arithmetic and iterative refinement of generalized eigenvalue problems*, SIAM J. Matrix Anal. Appl. 22, (2001), 1038-1057.

[66] F. Vana, P. Duben, S. Lang, T. Palmer, M. Leutbecher, D. Salmond, G. Carver, *Single precision in weather forecasting models: an evaluation with the IFS*, Mon. Weather Rev. 145 (2017), 495-502.

[67] H.R. Wilson and J.D. Cowan, *Excitatory and inhibitory interactions in localized populations of model neurons*, Biophys. J., 12 (1972), 1-24.

[68] Y. Zhao, T. Fukaya, L. Zhang, T. Iwashita, *Numerical Investigation into the Mixed Precision GMRES(m) Method Using FP64 and FP32*, Journal of Information Processing Vol.30, (2022), 525-537.

# 7 Appendix 1

<table><thead><tr><td>Parameter</td><td>Value</td><td>Unit</td></tr></thead><tbody><tr><td>k<sub>s</sub></td><td>0.1</td><td>unitless</td></tr><tr><td>η</td><td>0.01</td><td>per cell</td></tr><tr><td>k<sub>d</sub></td><td>0.0</td><td>h<sup>-1</sup></td></tr><tr><td>C</td><td>0.4</td><td>nM</td></tr><tr><td>b<sub>bmal0</sub></td><td>0.0</td><td>nM</td></tr><tr><td>p<sub>0</sub></td><td>4</td><td></td></tr><tr><td>ν<sub>1b</sub></td><td>9.0</td><td>nM h<sup>-1</sup></td></tr><tr><td>k<sub>1b</sub></td><td>1.0</td><td>nM</td></tr><tr><td>k<sub>1d</sub></td><td>0.12</td><td>h<sup>-1</sup></td></tr><tr><td>k<sub>1i</sub></td><td>0.56</td><td>nM</td></tr><tr><td>k<sub>2b</sub></td><td>0.3</td><td>nM<sup>-1</sup> h<sup>-1</sup></td></tr><tr><td>k<sub>2d</sub></td><td>0.05</td><td>h<sup>-1</sup></td></tr><tr><td>k<sub>2t</sub></td><td>0.24</td><td>h<sup>-1</sup></td></tr><tr><td>k<sub>3t</sub></td><td>0.02</td><td>h<sup>-1</sup></td></tr><tr><td>q</td><td>2</td><td>unitless</td></tr><tr><td>k<sub>3d</sub></td><td>0.12</td><td>h<sup>-1</sup></td></tr><tr><td>ν<sub>4b</sub></td><td>3.6</td><td>nM<sup>-1</sup> h<sup>-1</sup></td></tr><tr><td>r<sub>0</sub></td><td>3</td><td>unitless</td></tr><tr><td>k<sub>4b</sub></td><td>2.16</td><td>nM<sup>-1</sup> h<sup>-1</sup></td></tr><tr><td>k<sub>4d</sub></td><td>0.75</td><td>h<sup>-1</sup></td></tr><tr><td>k<sub>5b</sub></td><td>0.24</td><td>h<sup>-1</sup></td></tr><tr><td>k<sub>5d</sub></td><td>0.06</td><td>h<sup>-1</sup></td></tr><tr><td>k<sub>5t</sub></td><td>0.45</td><td>h<sup>-1</sup></td></tr><tr><td>k<sub>6t</sub></td><td>0.06</td><td>h<sup>-1</sup></td></tr><tr><td>k<sub>6d</sub></td><td>0.12</td><td>h<sup>-1</sup></td></tr><tr><td>k<sub>6a</sub></td><td>0.09</td><td>h<sup>-1</sup></td></tr><tr><td>k<sub>7a</sub></td><td>0.003</td><td>h<sup>-1</sup></td></tr><tr><td>k<sub>7d</sub></td><td>0.09</td><td>h<sup>-1</sup></td></tr><tr><td>τ<sub>0</sub></td><td>1.0</td><td>unitless</td></tr></tbody></table>
<table><tr><td>std_circadian_clock</td><td>0.05</td><td>unitless</td></tr><tr><td>c</td><td>0.01</td><td>nM</td></tr><tr><td>k<sub>impf</sub></td><td>4.0</td><td>h<sup>-1</sup></td></tr><tr><td>k<sub>0mpf</sub></td><td>6</td><td>h<sup>-1</sup></td></tr><tr><td>k<sub>1mpf</sub></td><td>0.05</td><td>nM</td></tr><tr><td>s</td><td>20.0</td><td>nM</td></tr><tr><td>d<sub>wee1</sub></td><td>5.0</td><td>h<sup>-1</sup></td></tr><tr><td>n<sub>0</sub></td><td>2</td><td>unitless</td></tr><tr><td>k<sub>actw</sub></td><td>1.0</td><td>h<sup>-1</sup></td></tr><tr><td>d<sub>w1</sub></td><td>1.0</td><td>nM</td></tr><tr><td>c<sub>w</sub></td><td>0.5</td><td>nM</td></tr><tr><td>k<sub>inactw</sub></td><td>200.0</td><td>h<sup>-1</sup></td></tr><tr><td>k<sub>1wee1</sub></td><td>0.5</td><td>nM</td></tr><tr><td>d<sub>w2</sub></td><td>1.0</td><td>h<sup>-1</sup></td></tr><tr><td>k<sub>act</sub></td><td>0.01</td><td>h<sup>-1</sup></td></tr><tr><td>target_period</td><td>20.0</td><td>intrinsic period of the cell cycle in hours</td></tr><tr><td>std_cell_cycle</td><td>0.1</td><td>relative variability of the cell cycle periods</td></tr></table>

Para_Cira

<caption>Table 1: The parameters for benchmark 1</caption>

The initial solution $y_0$ is given by :

$$ 
\begin{aligned}
y_1^{(i)}(0) &= 0.1 \times \text{rand}^{(i)} \\
y_2^{(i)}(0) &= 0.2 \times \text{rand}^{(i)} \\
y_3^{(i)}(0) &= 1.8 \times \text{rand}^{(i)} \\
y_4^{(i)}(0) &= 0.4 \times \text{rand}^{(i)} \\
y_5^{(i)}(0) &= 0.5 \times \text{rand}^{(i)} \\
y_6^{(i)}(0) &= 0.6 \times \text{rand}^{(i)} \\
y_7^{(i)}(0) &= 0.1 \times \text{rand}^{(i)} \\
y_8^{(i)}(0) &= 0.1 \times \text{rand}^{(i)} \\
y_9^{(i)}(0) &= 0.1 \times \text{rand}^{(i)} \\
y_{10}^{(i)}(0) &= 0.1 \times \text{rand}^{(i)}
\end{aligned}
$$

for $i = 1, \dots, d$. Parameters $\lambda_0, \lambda$ and $\tau$ are given by

$$ 
\begin{aligned}
\lambda_0 &= 97.4/\text{target\_period} \\
\lambda &= \text{rand}(\lambda_0, \text{std\_cell\_cycle} \times \lambda_0), \\
\tau &= \text{rand}(\tau_0, \text{std\_circadian\_clock} \times \tau_0).
\end{aligned}
$$

The vector *rand* is computed by the subroutine "*call random_number(rand)*."

# 8 Appendix 2

Supp:res
<figure><img src="image_7.png" /><figcaption>Figure 6: Numerical results for benchmark 1.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7: Numerical results for benchmark 2.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 8: Numerical results for benchmark 1. The vertical dashed lines in this Figure indicate that the algorithm **DPAINS** has not converged.</figcaption></figure>

BEN1_IRK46

<figure><img src="image_10.png" /><figcaption>Figure 9: Numerical results for benchmark 2.</figcaption></figure>

BEN2_IRK46