# Acceleration of convergence for numerical sequences

Clément Vincent, Mohamed Khelifi

► To cite this version:

Clément Vincent, Mohamed Khelifi. Acceleration of convergence for numerical sequences. CEA - Commissariat à l'énergie atomique et aux énergies alternatives. 2023. hal-04207550

HAL Id: hal-04207550

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

Submitted on 14 Sep 2023

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 CC0 1.0 - Universal - International License
# Acceleration of convergence for numerical sequences

Clément Vincent<fnref n="1" />, Mohamed Khelifi<fnref n="2" />

## Contents

<table><tr><td>Introduction</td><td>1</td></tr><tr><td>1 Acceleration of convergence</td><td>2</td></tr><tr><td>1.1 Order of convergence</td><td>2</td></tr><tr><td>1.2 Comparison between two sequences</td><td>4</td></tr><tr><td>2 The ε-algorithm</td><td>4</td></tr><tr><td>2.1 The scalar ε-algorithm</td><td>4</td></tr><tr><td>2.1.1 Presentation</td><td>4</td></tr><tr><td>2.1.2 Theory</td><td>9</td></tr><tr><td>2.1.3 Implementation of the algorithm</td><td>17</td></tr><tr><td>2.1.4 Application to iterative methods</td><td>21</td></tr><tr><td>2.2 The vectorial ε-algorithm</td><td>26</td></tr><tr><td>2.2.1 Presentation</td><td>26</td></tr><tr><td>2.2.2 Theory</td><td>27</td></tr><tr><td>2.2.3 Implementation of the algorithm</td><td>28</td></tr><tr><td>2.2.4 Application to the fixed point method</td><td>28</td></tr><tr><td>3 Application to concrete problems</td><td>30</td></tr><tr><td>3.1 Stationary diffusion problem</td><td>31</td></tr><tr><td>3.1.1 Problem 1 : κ constant</td><td>31</td></tr><tr><td>3.1.2 Problem 2 : κ discontinuous</td><td>39</td></tr><tr><td>3.2 Non-stationary diffusion problem</td><td>43</td></tr><tr><td>3.2.1 Problem 3 : κ constant</td><td>44</td></tr><tr><td>3.2.2 Problem 4 : κ discontinuous</td><td>45</td></tr><tr><td>4 Conclusion</td><td>48</td></tr><tr><td>Appendix A Laptop information</td><td>49</td></tr><tr><td>Appendix B Particular rule</td><td>50</td></tr><tr><td>Appendix C Graphs for the problem 1 on the mesh 512 × 512</td><td>52</td></tr><tr><td>Appendix D Graphs for the problem 2 on the mesh 512 × 512</td><td>53</td></tr><tr><td>References</td><td>54</td></tr></table>

<footnote n="1">¹Université de Lille, Villeneuve d'Ascq, France, clement.vincent2@cea.fr</footnote>

<footnote n="2">²CEA, DAM, DIF, F-91297 Arpajon, France, mohamed.khelifi@cea.fr</footnote>
# Introduction

Nowadays, numerical simulation has an important role in scientific research. Indeed, the development of computers has led to the resolution of a growing number of increasingly complex and time-consuming problems. To reduce these costs, High-Performance Computing (HPC) was born. Its aim is to perform complex computations and massive data processing at very high speeds, by combining the power of several hundred processors. In some cases, HPC can be combined with other tools to accelerate problem solving.

Indeed, some problems require the resolution of non-linear systems. Very often, the methods used for these resolutions are so-called iterative methods (Newton's method, fixed point method, etc...). When these methods converge, it sometimes happens that the convergence is slow. Convergence acceleration methods are used to improve the speed of the convergence.

The acceleration of convergence is a field in its own right in numerical analysis. Examples include Aitken's $\Delta^2$ method and Richardson's extrapolation process, respectively popularized by mathematicians Alexander Aitken and Lewis Fry Richardson at the beginning of the XXth century. In the second half of the XXth century, other methods came into being, such as the $\rho$-algorithm and the $\varepsilon$-algorithm, due to the mathematician Peter Wynn, or the $\theta$-algorithm, due to the mathematician Claude Brezinski. Each of these methods are presented in [Bre77] and [Bre78]. We will study and focus on just one algorithm: the $\varepsilon$-algorithm.

We'll give some definitions to introduce the acceleration of convergence for sequences. Then, we will present the $\varepsilon$-algorithm. Firstly, we'll study the $\varepsilon$-algorithm for accelerating scalar sequences. Secondly, we'll look at the $\varepsilon$-algorithm for accelerating vector sequences. Finally, we'll apply this algorithm to concrete problems. All programs were run on the same laptop (Appendix [A]). The computations were not performed on a supercomputer. But, when it was possible, the resources of the laptop were used to the maximum, by parallelizing the codes over 8 cores.

When we solve a numerical problem by using an iterative method, we want to generate a sequence $(S_n)_n$ which converges to the solution of the problem when $n$ tends to infinity.

The principle of a convergence acceleration method is to transform a sequence $(S_n)_n$ into a sequence $(V_n)_n$ of the same nature, but converging faster to the limit when it's possible. The interest of convergence acceleration methods will depend on the speed with which the initial sequence converges. Indeed, if $(S_n)_n$ converges slowly, the gain in terms of computation and therefore computation time can be appreciable.

The aim of this study is to apply a convergence acceleration method to a sequence generated by an iterative method, in order to speed up the convergence of the method to the solution of the problem. To do so, we will use the $\varepsilon$-algorithm. The choice of this method holds in the simplicity of its formulation and implementation.

To achieve this objective, we will first give definitions to introduce the acceleration of the convergence of sequences. Secondly, we will present the scalar $\varepsilon$-algorithm. This algorithm accelerates the convergence of real numbers sequences. We will apply this algorithm to various iterative methods and analyse the results obtained. Then, we will apply the $\varepsilon$-algorithm to solve non-linear equations on the form $F(\mathbf{X}) = 0$, where we have $F: \mathbb{R}^n \rightarrow \mathbb{R}^n$, with $n \in \mathbb{N}^*$. We will therefore need to accelerate vector sequences. We'll look at the vectorial $\varepsilon$-algorithm. Finally, we will apply the $\varepsilon$-algorithm to stationary and non-stationary diffusion equations.
# 1 Acceleration of convergence

In this part, we'll give some definitions which characterise the order of convergence of a sequence and which allow us to compare the speed of convergence between two sequences [Bre77].

In the two following sections, we will only consider sequences of positive or zero real numbers which tend to zero. In our context, it isn't restrictive. In fact, we want to accelerate the convergence of sequences. So, let $(E, d)$ be a metric space and $(S_n)_n$ a sequence of elements of $E$ which converges to $S$. Quantities $d(S_n, S)$ are positive or zero real numbers and the sequence $(d(S_n, S))_n$ converges to zero.

## 1.1 Order of convergence

To talk about the acceleration of convergence of a sequence, we need to be able to compare the speed of convergence of two sequences. So let us begin by recalling the comparison relationships between two sequences.

Let $(u_n)_n$ and $(v_n)_n$ two sequences of real numbers which tends to zero when $n$ tends to infinity.

We denote,

$$v_n = O(u_n)$$

if $\exists N \in \mathbb{N}$ and $C > 0$, $\forall n > N$, $|v_n| \le C|u_n|$.

And we denote,

$$v_n = o(u_n)$$

if $\forall \varepsilon > 0, \exists N \in \mathbb{N}, \forall n > N, |v_n| < \varepsilon|u_n|$. In other words, $\lim_{n \to \infty} \frac{v_n}{u_n} = 0$.

The order of convergence of a sequence is defined by :

**Definition 1.** The sequence $(u_n)_n$ is of order $r \in \mathbb{R}_+^*$ if :

$$u_{n+1} = O(u_n^r) \text{ and } u_n^r = O(u_{n+1}).$$

**Remark 1.** It means there exists $A$ and $B$, with $0 < A \le B < +\infty$, such that, $\forall n > N$,

$$A \le \frac{u_{n+1}}{u_n^r} \le B.$$

In addition, we can prove that :

**Theorem 1.** If it exists, $r$ is unique.

*Proof.* Let suppose there exists two strictly positive real numbers $p$ and $r$ with $p \neq r$ such that the sequence is of order $p$ and $r$.

We have,

$$u_{n+1} = O(u_n^r) \text{ and } u_n^r = O(u_{n+1}),$$

and,

$$u_{n+1} = O(u_n^p) \text{ and } u_n^p = O(u_{n+1}).$$
So, there exists $N_{r_1}, N_{r_2}, N_{p_1}, N_{p_2} \in \mathbb{N}$ and $C_1, C_2, C_3, C_4 > 0$ such that,

$$ \forall n > N_{r_1}, u_{n+1} < C_1 u_n^r \quad \text{and} \quad \forall n > N_{r_2}, u_n^r < C_2 u_{n+1}, $$

$$ \forall n > N_{p_1}, u_{n+1} < C_4 u_n^p \quad \text{and} \quad \forall n > N_{p_2}, u_n^p < C_3 u_{n+1}. $$

We can write that, for all $n > \max(N_{r_1}, N_{r_2}, N_{p_1}, N_{p_2})$,

$$ u_{n+1} \le C_1 u_n^r = C_1 u_n^{r-p} u_n^p \le C_1 C_3 u_n^{r-p} u_{n+1}. $$

So, we have,

$$ 1 \le C_1 C_3 u_n^{r-p}. $$

Let suppose $r > p$.

As $u_n$ is a sequence which tends to zero and $r > p$, we have that $u_n^{r-p}$ tends to zero when $n$ tends to infinity. It's absurd because of the result above.

Now, let suppose $r < p$. We can write that, for all $n > \max(N_{r_1}, N_{r_2}, N_{p_1}, N_{p_2})$,

$$ u_{n+1} \le C_4 u_n^p = C_4 u_n^{p-r} u_n^r \le C_4 C_2 u_n^{p-r} u_{n+1}. $$

So, we have,

$$ 1 \le C_4 C_2 u_n^{p-r}. $$

In using the same argument than the previous one but with $r < p$, we can say that $u_n^{p-r}$ tends to zero when $n$ tends to infinity. It's also absurd.

Finally, $r$ and $p$ must be equal.

□

There exists another definition to characterise the order of a sequence which is :

**Definition 2.** If we have $u_{n+1} = O(u_n^r)$, we say that the sequence $(u_n)_n$ is at least of order $r$.

**Remark 2.** The previous definition doesn't guarantee the uniqueness of $r$.

We define another term which can help us to know how converges a sequence.

**Definition 3.** The asymptotic error coefficient is the number

$$ C = \limsup_{n \to \infty} \frac{u_{n+1}}{u_n^r}. $$

The notions of order and asymptotic error coefficient are not purely theoretical. They are closely related to the number of exact digits obtained. Since the sequence $(u_n)_n$ is positive and converges to zero, $u_n$ represents the absolute error.

Let $R = -\log_{10} C$.

It shows in [Bre77] (pages 4-5) that if $r = 1$, we add about $R$ exact significant digits by going from $u_n$ to $u_{n+1}$. However, if $r > 1$, the number of exact significant digits is approximately multiplied by $r$ by going from $u_n$ to $u_{n+1}$. It shows it's interesting to have an order of convergence bigger than 1.
## 1.2 Comparison between two sequences

Let $(u_n)_n$ and $(v_n)_n$ two sequences of positive real numbers which converge to zero. We will give two definitions which allow us to compare their convergence speed.

**Definition 4.** $(u_n)_n$ converges like $(v_n)_n$ if:

$$u_n = O(v_n) \quad \text{and} \quad v_n = O(u_n).$$

Moreover, we can say that $(v_n)_n$ converges better than $(u_n)_n$ if the real number $C$ defined by

$$C = \limsup_{n \to \infty} \frac{v_n}{u_n}$$

is strictly inferior to one.

**Definition 5.** $(v_n)_n$ converges faster than $(u_n)_n$ if $v_n = o(u_n)$.

A method which accelerates the convergence of a sequence is defined by:

**Definition 6.** Let $T$ a method which transforms the sequence $(u_n)_n$ to another sequence $(v_n)_n$. If $v_n = o(u_n)$, we say that we accelerate the convergence and the method $T$ is a convergence acceleration method.

# 2 The $\varepsilon$-algorithm

The $\varepsilon$-algorithm is due to the English mathematician Peter Wynn (1931-2017). His research was mainly focused on the theory of approximation. In particular, he studied Padé's approximants and their applications to numerical methods for improving the order of convergence of sequences of real numbers.

## 2.1 The scalar $\varepsilon$-algorithm

In this section, we'll restrict ourselves to the study of real-values sequences. These sequences will be denoted by $(S_n)_n$, where $n \in \mathbb{N}$. If $(S_n)_n$ converges, we'll denote by $S$ its finite limit.

Now, let us take another sequence $(V_n)_n$ which also converges to $S$. We say that $(V_n)_n$ converges faster than $(S_n)_n$ if

$$\lim_{n \to \infty} \frac{|V_n - S|}{|S_n - S|} = 0.$$

By analogy with the previous section, in denoting by $(v_n)_n$ the sequence defined by $v_n = |V_n - S|$, $\forall n \in \mathbb{N}$, and $(u_n)_n$ the sequence defined by $u_n = |S_n - S|$, $\forall n \in \mathbb{N}$, we have $v_n = o(u_n)$.

### 2.1.1 Presentation

The $\varepsilon$-algorithm looks like this. We have to compute quantities with two indices $\varepsilon_k^{(n)}$ in using the recurrence:

$$\begin{align}
\varepsilon_{-1}^{(n)} &= 0 & \varepsilon_{0}^{(n)} &= S_n & n &= 0, 1, \dots \nonumber \\
\varepsilon_{k+1}^{(n)} &= \varepsilon_{k-1}^{(n+1)} + \frac{1}{\varepsilon_k^{(n+1)} - \varepsilon_k^{(n)}} & n, k &= 0, 1, \dots \tag{1}
\end{align}$$
These values can be put in a double-entry table called the $\epsilon$ table. The lower index remains constant in a column, while the upper index remains constant on a descending diagonal :

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

Elements of the sequence $(S_n)_n$ are located in the second column.

In theory, the $\epsilon$ table is infinite. However, in practice, we only know or calculate a limited number of elements in the sequence $(S_n)_n$. For example, if we only know seven elements of the sequence $(S_n)_n$, then, the $\epsilon$ table will be the same as the previous one without the dotted lines.

Looking at the previous $\epsilon$ table, we note that the expression (1) connects the quantities located at the four vertices of a rhombus as follows :

<figure><img src="image_3.png" /></figure>
The previous diagram tells us that if we know $\varepsilon_{k-1}^{(n+1)}$, $\varepsilon_k^{(n)}$ and $\varepsilon_k^{(n+1)}$, the recurrence (1) allows us to determine $\varepsilon_{k+1}^{(n)}$. In this sense, we say that the $\varepsilon$-algorithm is a rhombus algorithm.

The recurrence (1) allows us to progress from left to right and from top to bottom in the $\varepsilon$ table from the initial conditions $\varepsilon_{-1}^{(n)} = 0$ and $\varepsilon_0^{(n)} = S_n$ for $n = 0, 1, \dots$. Indeed, if we know $S_0$ and $S_1$, we can compute $\varepsilon_1^{(0)}$. Then, if we know the next element in the initial sequence $S_2$, we can compute $\varepsilon_1^{(1)}$ and $\varepsilon_2^{(0)}$:

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

**Remark 3.** For each new element of the initial sequence $(S_n)_n$, we can construct an ascending diagonal on the $\varepsilon$ table.

**Remark 4.** To compute $\varepsilon_k^{(n)}$, we need to know $S_n, S_{n+1}, \dots, S_{n+k}$, so $k+1$ elements of the initial sequence.

**Remark 5.** From the mathematical relation (1), we have:

$$ \varepsilon_2^{(n)} = \frac{S_n S_{n+2} - S_{n+1}^2}{S_{n+2} - 2S_{n+1} + S_n} \quad n = 0, 1, \dots $$

This is the Aitken's $\Delta^2$ process, which is a convergence acceleration method.

The $\varepsilon$-algorithm from P. Wynn is a generalisation of the Aitken's $\Delta^2$ process [Bre77] [Bre78].

Indeed, Aitken's $\Delta^2$ process corresponds to the specific column constituted by the $\varepsilon_2^{(n)}$ elements from the $\varepsilon$ table. The sequence $(\varepsilon_2^{(n)})_n$ is the accelerated sequence of $(S_n)_n$ obtained by the $\varepsilon$-algorithm.

The theory of the $\varepsilon$-algorithm shows that quantities with an odd lower index (the odd columns) are only intermediate calculations ([Bre78], page 11). Only quantities with an even lower index (the even columns) can be integrated in the new sequence created by the $\varepsilon$-algorithm. With these quantities, we compare the convergence of the new sequence and the initial one. By eliminating the odd columns, a table that contains only quantities with an even lower index is obtained:
$\epsilon_0^{(0)}$

$\epsilon_{-2}^{(2)} = \infty$

$\epsilon_0^{(1)}$

$\epsilon_2^{(0)}$

$\epsilon_{-2}^{(3)} = \infty$

$\epsilon_0^{(2)}$

$\epsilon_2^{(1)}$

$\epsilon_4^{(0)}$

$\epsilon_{-2}^{(4)} = \infty$

$\epsilon_0^{(3)}$

$\epsilon_2^{(2)}$

$\epsilon_{-2}^{(5)} = \infty$

$\epsilon_0^{(4)}$

The initial conditions $\epsilon_{-1}^{(n)} = 0$ in (1) are replaced by the initial conditions $\epsilon_{-2}^{(n)} = \infty$ for all $n$. The even lower index quantities in the previous table can be calculated without using the odd lower index quantities.

Indeed, the following recurrence, known as the rule of the cross, is used:

$$ \epsilon_{-2}^{(n)} = \infty \qquad \epsilon_{0}^{(n)} = S_n \qquad n = 0, 1, \dots $$

(2)

$$ \frac{1}{\epsilon_{2k+2}^{(n-1)} - \epsilon_{2k}^{(n)}} + \frac{1}{\epsilon_{2k-2}^{(n+1)} - \epsilon_{2k}^{(n)}} = \frac{1}{\epsilon_{2k}^{(n+1)} - \epsilon_{2k}^{(n)}} + \frac{1}{\epsilon_{2k}^{(n-1)} - \epsilon_{2k}^{(n)}} \qquad n, k = 0, 1, \dots $$

This expression can be found in using (1).

In writing the recurrence (1) for $n$ and $n-1$,

$$ (a) \quad : \quad \epsilon_{k+1}^{(n)} = \epsilon_{k-1}^{(n+1)} + \frac{1}{\epsilon_k^{(n+1)} - \epsilon_k^{(n)}} $$

$$ (b) \quad : \quad \epsilon_{k+1}^{(n-1)} = \epsilon_{k-1}^{(n)} + \frac{1}{\epsilon_k^{(n)} - \epsilon_k^{(n-1)}} $$

By subtracting (b) from (a),

$$ (c) \quad : \quad \epsilon_{k+1}^{(n)} - \epsilon_{k+1}^{(n-1)} = \epsilon_{k-1}^{(n+1)} - \epsilon_{k-1}^{(n)} + \frac{1}{\epsilon_k^{(n+1)} - \epsilon_k^{(n)}} - \frac{1}{\epsilon_k^{(n)} - \epsilon_k^{(n-1)}} $$

But,

$$ \epsilon_k^{(n)} = \epsilon_{k-2}^{(n+1)} + \frac{1}{\epsilon_{k-1}^{(n+1)} - \epsilon_{k-1}^{(n)}} $$

It implies that,

$$ \epsilon_{k-1}^{(n+1)} - \epsilon_{k-1}^{(n)} = \frac{1}{\epsilon_k^{(n)} - \epsilon_{k-2}^{(n+1)}} $$
In the same way, by applying the recurrence process of the $\epsilon$-algorithm to $\epsilon_{k+2}^{(n-1)}$,

$$ \epsilon_{k+2}^{(n-1)} = \epsilon_k^{(n)} + \frac{1}{\epsilon_{k+1}^{(n)} - \epsilon_{k+1}^{(n-1)}} $$

It implies that,

$$ \epsilon_{k+1}^{(n)} - \epsilon_{k+1}^{(n-1)} = \frac{1}{\epsilon_{k+2}^{(n-1)} - \epsilon_k^{(n)}} $$

By replacing the previous values in (c), we have

$$ \frac{1}{\epsilon_{k+2}^{(n-1)} - \epsilon_k^{(n)}} = \frac{1}{\epsilon_k^{(n)} - \epsilon_{k-2}^{(n+1)}} + \frac{1}{\epsilon_k^{(n+1)} - \epsilon_k^{(n)}} - \frac{1}{\epsilon_k^{(n)} - \epsilon_k^{(n-1)}} $$

which is equivalent to the rule of the cross,

$$ \frac{1}{\epsilon_{k+2}^{(n-1)} - \epsilon_k^{(n)}} + \frac{1}{\epsilon_{k-2}^{(n+1)} - \epsilon_k^{(n)}} = \frac{1}{\epsilon_k^{(n+1)} - \epsilon_k^{(n)}} + \frac{1}{\epsilon_k^{(n-1)} - \epsilon_k^{(n)}} $$

The name "the rule of the cross" comes from the fact that this formula links quantities located at the four vertices and the center of a cross :

$$ \epsilon_{2k}^{(n-1)} $$

$$ \epsilon_{2k-2}^{(n+1)} $$

$$ \epsilon_{2k}^{(n)} $$

$$ \epsilon_{2k+2}^{(n-1)} $$

$$ \epsilon_{2k}^{(n+1)} $$

The computation of $\epsilon_{2k+2}^{(n-1)}$ is done using the recurrence (2), so from $\epsilon_{2k-2}^{(n+1)}$, $\epsilon_{2k}^{(n-1)}$, $\epsilon_{2k}^{(n)}$ and $\epsilon_{2k}^{(n+1)}$. We progress through the $\epsilon$ table with only the even-numbered columns, from the left to the right and from the top to the bottom, starting from the initial conditions $\epsilon_{-2}^{(n)} = \infty$ and $\epsilon_0^{(n)} = S_n$ for $n = 0,1, \dots$. It is possible to symbolise the previous diagram with the first letters of the cardinal points to simplify the notations :

<figure><img src="image_5.png" /></figure>
and the rule of the cross is rewriting,

$$ \frac{1}{N-C} + \frac{1}{S-C} = \frac{1}{W-C} + \frac{1}{E-C}. \qquad (3) $$

**Remark 6.** The computation of $\varepsilon_{2k}^{(n)}$, for $n$ and $k$ fixed, requires the knowledge of $S_n$, $S_{n+1}$, ..., $S_{n+2k}$.

Now, let us determine the number of arithmetic operations of the $\varepsilon$-algorithm. To do this, let us take the expression (1). As it has already been said, we are just interested in quantities in the $\varepsilon$ table with an even inferior index to see the convergence of the new sequence generated by the $\varepsilon$-algorithm. Moreover, to accelerate the convergence of a sequence, we'll see that the most interesting new sequence to take may be the quantities located on the first downward diagonal $(\varepsilon_n^{(0)})_n$.

Let us check the number of arithmetic operations used to compute an element $\varepsilon_{2N}^{(0)}$, with $N \in \mathbb{N}$. As seen in the remark 4, we need to know $2N + 1$ elements of the initial sequence to compute $\varepsilon_{2N}^{(0)}$. So, in using the first $\varepsilon$ table, the number of $\varepsilon_k^{(n)}$ that must be computed is equal to

$$ \sum_{i=1}^{2N} i = \frac{2N(2N+1)}{2} = N(2N+1). $$

Computing each of these terms requires 2 additions/subtractions and 1 multiplication/division. This gives a total of arithmetic operations equals to

$2N(2N+1)$ additions,

and,

$N(2N+1)$ divisions.

It means that the total number of arithmetic operations will be proportional to $N^2$. In other words, the complexity of the algorithm is quadratic : $O(N^2)$.

### 2.1.2 Theory

The $\varepsilon$-algorithm is a generalisation of Aitken's $\Delta^2$ process [Bre77] [Bre78]. This generalisation was made possible by the Shanks transformation, discovered by the mathematician Daniel Shanks and described below. The Shanks transformation is a non-linear method for accelerating the convergence of numerical sequences. It requires the values of some determinants, so it is difficult to implement. Indeed, computing a determinant on a computer requires a large number of arithmetic operations and a good precision is rarely achieved.

The $\varepsilon$-algorithm from P. Wynn is a recursive procedure that avoids computing the determinants which appear in the Shanks transformation.

Recall we denote by $(S_n)_n$ the initial sequence which converges to the finite limit $S$.
## The Shanks transformation

The Shanks transformation consists in transforming the sequence $(S_n)_n$ into a sequence denoted by $(e_k(S_n))_n$. This transformation is defined by the ratio of two determinants

$$e_k(S_n) = \frac{\begin{vmatrix} S_n & \cdots & S_{n+k} \\ S_{n+1} & \cdots & S_{n+k+1} \\ \vdots & \ddots & \vdots \\ S_{n+k} & \cdots & S_{n+2k} \end{vmatrix}}{\begin{vmatrix} \Delta^2 S_n & \cdots & \Delta^2 S_{n+k-1} \\ \vdots & \ddots & \vdots \\ \vdots & \ddots & \vdots \\ \Delta^2 S_{n+k-1} & \cdots & \Delta^2 S_{n+2k-2} \end{vmatrix}} = \frac{\begin{vmatrix} S_n & \cdots & S_{n+k} \\ \Delta S_n & \cdots & \Delta S_{n+k} \\ \vdots & \ddots & \vdots \\ \Delta S_{n+k-1} & \cdots & \Delta S_{n+2k-1} \end{vmatrix}}{\begin{vmatrix} 1 & \cdots & 1 \\ \Delta S_n & \cdots & \Delta S_{n+k} \\ \vdots & \ddots & \vdots \\ \Delta S_{n+k-1} & \cdots & \Delta S_{n+2k-1} \end{vmatrix}}$$

with

$$\Delta^0 S_n = S_n \text{ for all } n \text{ and } \Delta^{k+1} S_n = \Delta^k S_{n+1} - \Delta^k S_n \text{ for } n, k = 0, 1, \dots$$

The goal of D. Shanks was to construct a transformation that would give the value of $S$ when the sequence $(S_n - S)_n$ satisfies a linear, homogeneous difference equation with constant coefficients and of order $k$ [Bel]. In other words, Shanks wanted to construct a transformation that gives $S$ for any sequences $(S_n)_n$ satisfying

$$\sum_{i=0}^{k} a_i (S_{n+i} - S) = 0 \quad \forall n \in \mathbb{N},$$

with $k \in \mathbb{N}$ and $\forall i \in [0, k-1]$, $a_i \in \mathbb{R}$ and $a_k \in \mathbb{R}^*$.

The condition to obtain the limit $S$ with the Shanks transformation is given by the following theorem:

**Theorem 2.** A necessary and sufficient condition to have $e_k(S_n) = S$ for all $n \ge N$ is that the sequence $(S_n)_n$ satisfies:

$$\sum_{i=0}^{k} a_i (S_{n+i} - S) = 0 \quad \forall n \ge N$$

with $\sum_{i=0}^{k} a_i \neq 0$.

*Proof.* Available in [Bre77], pages 40-41. □

## Properties of the $\varepsilon$-algorithm

The fundamental property of the $\varepsilon$-algorithm is the one which directly links it to the Shanks transformation. This fundamental property is due to P. Wynn and is the following one:
**Theorem 3.**

$$ \varepsilon_{2k}^{(n)} = e_k(S_n) \quad \text{and} \quad \varepsilon_{2k+1}^{(n)} = \frac{1}{e_k(\Delta S_n)} \qquad n, k = 0, 1, \dots $$

*Proof.* Available in [Bre77], pages 44-46. □

Thus, we see that the even columns of the $\varepsilon$-algorithm are identical to the Shanks transformation. We can also see that the odd lower index quantities are the inverses of the Shanks transformation applying on the sequence $(\Delta S_n = S_{n+1} - S_n)_n$. From a practical point of view (when we are trying to speed up the convergence of a sequence), the odd columns of the $\varepsilon$-algorithm are only intermediate computations with no particular meaning.

The $\varepsilon$-algorithm is a non-linear method from sequence to sequence. Let us consider two convergent sequences $(S_{1,n})_n$ and $(S_{2,n})_n$ which respectively converge to $S_1$ and $S_2$. We denote by $(e_k(S_{1,n}))_n^k$ and $(e_k(S_{2,n}))_n^k$ their accelerated sequence obtained by the Shanks transformation. We have, for $k$ and $n$ fixed,

$$ e_k(S_{1,n} + S_{2,n}) \neq e_k(S_{1,n}) + e_k(S_{2,n}). $$

So, if we want to compute the terms of the accelerated sequence of $(S_{1,n} + S_{2,n})_n$, we can't add the terms of the accelerated sequence of $(S_{1,n})_n$ to the terms of the accelerated sequence of $(S_{2,n})_n$.

**Example 1.** Let $S_{1,n} = \frac{n}{(n+1)^2} + 5$ which converges to 5 and $S_{2,n} = \frac{1}{n+1}$ which converges to 0. The sequence $(S_{1,n} + S_{2,n})_n$ also converges to 5. The $\varepsilon$-algorithm gives :

<table><thead><tr><th colspan="5">Example of non-linearity for the ε-algorithm</th></tr><tr><th>ε<sub>2k</sub><sup>(0)</sup></th><th>e<sub>k</sub>(S<sub>1,0</sub>)</th><th>e<sub>k</sub>(S<sub>2,0</sub>)</th><th>e<sub>k</sub>(S<sub>1,0</sub>) + e<sub>k</sub>(S<sub>2,0</sub>)</th><th>e<sub>k</sub>(S<sub>1,0</sub>+S<sub>2,0</sub>)</th></tr></thead><tbody><tr><td>ε<sub>2</sub><sup>(0)</sup></td><td>5.250000e+00</td><td>2.499999e-01</td><td>5.5e+00</td><td>5.431818e+00</td></tr><tr><td>ε<sub>4</sub><sup>(0)</sup></td><td>5.163867e+00</td><td>1.111111e-01</td><td>5.274978e+00</td><td>5.207078e+00</td></tr><tr><td>ε<sub>6</sub><sup>(0)</sup></td><td>5.063802e+00</td><td>6.250000e-02</td><td>5.126302e+00</td><td>5.112758e+00</td></tr><tr><td>ε<sub>8</sub><sup>(0)</sup></td><td>5.033262e+00</td><td>4.000000e-02</td><td>5.073262e+00</td><td>5.065332e+00</td></tr></tbody></table>

However, there is linearity when we go to the limit. Indeed, the following general result states :

**Theorem 4.** Let $(S_{1,n})_n$, $(S_{2,n})_n$ and $(S_{1,n}+S_{2,n})_n$ three convergent sequences which respectively converge to $S_1$, $S_2$ and $S_3$. Let $T$ a transformation from the space of convergent sequences to the same space, and such that their accelerated convergent sequences $(T(S_{1,n}))_n$, $(T(S_{2,n}))_n$ and $(T(S_{1,n}+S_{2,n}))_n$ converge respectively to the limits $S_1$, $S_2$ and $S_3$. Then,

$$ \lim_{n \to \infty} T(S_{1,n} + S_{2,n}) = \lim_{n \to \infty} T(S_{1,n}) + \lim_{n \to \infty} T(S_{2,n}). $$

*Proof.* We have,

$$ \lim_{n \to \infty} (S_{1,n} + S_{2,n}) = \lim_{n \to \infty} S_{1,n} + \lim_{n \to \infty} S_{2,n} = \lim_{n \to \infty} T(S_{1,n}) + \lim_{n \to \infty} T(S_{2,n}). $$
But,

$$ \lim_{n \to \infty} (S_{1,n} + S_{2,n}) = \lim_{n \to \infty} (T(S_{1,n} + S_{2,n})). $$

Then,

$$ \lim_{n \to \infty} T(S_{1,n} + S_{2,n}) = \lim_{n \to \infty} T(S_{1,n}) + \lim_{n \to \infty} T(S_{2,n}). $$

□

There is also linearity for multiplication by a constant and for the addition of a constant to all the terms of a sequence. Expressed in Shanks notation, this property can be written as:

**Theorem 5.** Let $a \neq 0$ and $b$ two arbitrary real constants, $(S_n)_n$ a convergent sequence and $(V_n)_n$ a sequence defined by $V_n = aS_n + b$. Then:

$$ e_k(V_n) = ae_k(S_n) + b $$

$$ \frac{1}{e_k(\Delta V_n)} = \frac{1}{ae_k(\Delta S_n)}. $$

*Proof.* To show the first equality, let use the expression of the Shanks transformation with determinants:

$$ e_k(V_n) = \frac{\begin{vmatrix} V_n & \cdots & V_{n+k} \\ \Delta V_n & \cdots & \Delta V_{n+k} \\ \vdots & \cdots & \vdots \\ \Delta V_{n+k-1} & \cdots & \Delta V_{n+2k-1} \end{vmatrix}}{\begin{vmatrix} 1 & \cdots & 1 \\ \Delta V_n & \cdots & \Delta V_{n+k} \\ \vdots & \cdots & \vdots \\ \Delta V_{n+k-1} & \cdots & \Delta V_{n+2k-1} \end{vmatrix}} = \frac{\begin{vmatrix} aS_n + b & \cdots & aS_{n+k} + b \\ a\Delta S_n & \cdots & a\Delta S_{n+k} \\ \vdots & \cdots & \vdots \\ a\Delta S_{n+k-1} & \cdots & a\Delta S_{n+2k-1} \end{vmatrix}}{\begin{vmatrix} 1 & \cdots & 1 \\ a\Delta S_n & \cdots & a\Delta S_{n+k} \\ \vdots & \cdots & \vdots \\ a\Delta S_{n+k-1} & \cdots & a\Delta S_{n+2k-1} \end{vmatrix}} $$

⇔

$$ e_k(V_n) = \frac{a^k \begin{vmatrix} S_n & \cdots & S_{n+k} \\ \Delta S_n & \cdots & \Delta S_{n+k} \\ \vdots & \cdots & \vdots \\ \Delta S_{n+k-1} & \cdots & \Delta S_{n+2k-1} \end{vmatrix}}{a^{k-1} \begin{vmatrix} 1 & \cdots & 1 \\ \Delta S_n & \cdots & \Delta S_{n+k} \\ \vdots & \cdots & \vdots \\ \Delta S_{n+k-1} & \cdots & \Delta S_{n+2k-1} \end{vmatrix}} + \frac{b \begin{vmatrix} b & \cdots & b \\ \Delta S_n & \cdots & \Delta S_{n+k} \\ \vdots & \cdots & \vdots \\ \Delta S_{n+k-1} & \cdots & \Delta S_{n+2k-1} \end{vmatrix}}{b \begin{vmatrix} 1 & \cdots & 1 \\ \Delta S_n & \cdots & \Delta S_{n+k} \\ \vdots & \cdots & \vdots \\ \Delta S_{n+k-1} & \cdots & \Delta S_{n+2k-1} \end{vmatrix}} = ae_k(S_n) + b. $$

We find the second equality by applying the definition of Shanks transformation to $\Delta V_n$.

□
**Remark 7.** In denoting $\varepsilon_{2k}^{(n)}$ the $\varepsilon$-algorithm terms for the sequence $(S_n)_n$ and $\tilde{\varepsilon}_{2k}^{(n)}$ those for the sequence $(V_n)_n$, the above equations in theorem 5 can be rewritten as $\tilde{\varepsilon}_{2k}^{(n)} = a\varepsilon_{2k}^{(n)} + b$ and $\tilde{\varepsilon}_{2k+1}^{(n)} = \varepsilon_{2k+1}^{(n)}/a$.

Convergence theorems

The following results comes from [Bre77] and [Bre78], where we can find convergence theorems on the $\varepsilon$-algorithm concerning :

- convergence theorems for sequences of well-defined form, for example : $S_n = S + \sum_{i=1}^{\infty} \alpha_i \lambda_i^n$,

- convergence theorems for classes of sequences with particular properties,

- convergence theorems for sequences whose formation law is known, like $S_{n+1} = f(S_n, S_{n-1}, \dots, S_{n-k})$,

- convergence theorems for the $2k$ column of the $\varepsilon$ table when some properties of the previous columns are known.

In this section, we will restrict ourselves to convergence results for certain classes of sequences : totally monotonic sequences and totally oscillating sequences.

Let us start by defining what a totally monotonic sequence is :

**Definition 7.** The sequence $(S_n)_n$ is said to be totally monotonic if :

$$(-1)^k \Delta^k S_n \geq 0 \quad n, k = 0, 1, \dots$$

In other words, $(S_n)_n$ is totally monotonic if :

$$ \begin{aligned} S_0 &\geq S_1 \geq S_2 \geq \dots \geq 0 \\ \Delta S_0 &\leq \Delta S_1 \leq \Delta S_2 \leq \dots \leq 0 \\ \Delta^2 S_0 &\geq \Delta^2 S_1 \geq \dots \geq 0 \end{aligned} $$

and so on.

**Remark 8.** In theoretical physics, totally monotonic sequences often appear as coefficients of series

$$f(x) = \sum_{i=0}^{\infty} (-1)^i c_i x^i$$

where $(c_n)_n$ is a totally monotone sequence. That's why totally monotonic sequences can be an important class of sequences.

For totally monotonic sequences, the following convergence result can be demonstrated :
**Theorem 6.** If we apply the ε-algorithm to a sequence $(S_n)_n$ which converges to S and if there exists two constants $a \neq 0$ and $b$ such that the sequence $(aS_n + b)_n$ is completely monotone, then

$$ \lim_{n \to \infty} \varepsilon_{2k}^{(n)} = S \quad k = 0, 1, \dots, \text{fixed} $$

$$ \lim_{k \to \infty} \varepsilon_{2k}^{(n)} = S \quad n = 0, 1, \dots, \text{fixed} $$

*Proof.* Available in [Bre77] page 72 (Theorem 49) and page 75 (Theorem 50). □

**Remark 9.** In the particular case where $a = 1$ and $b = 0$, it means that the sequence $(S_n)_n$ is completely monotonic and the ε-algorithm applied to this sequence converges in $n$ and $k$ to the limit $S$. It converges in $n$ for even columns. It converges in $k$ along the descending diagonals, in keeping only the even lower indices.

We can also establish the following theorem related to the ε-algorithm :

**Theorem 7.** If we apply the ε-algorithm to a totally monotonic sequence $(S_n)_n$, then, for any $n, k \in \mathbb{N}$ fixed:

$$ 0 \le \varepsilon_{2k+2}^{(n)} \le \varepsilon_{2k}^{(n)} $$

$$ 0 \le \varepsilon_{2k}^{(n+1)} \le \varepsilon_{2k}^{(n)} $$

$$ 0 \le \varepsilon_{2k+2}^{(n)} \le \varepsilon_{2k}^{(n+1)} $$

$$ 0 \le \varepsilon_{2k+2}^{(n)} \le \varepsilon_{2k}^{(n+2)} $$

*Proof.* Available in [Bre77] page 72 (Lemma 7), page 74 (Lemma 9) and page 75 (Lemma 9 bis). □

The theorem can be interpreted by the following diagram :

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

where the arrows go from the the furthest term to the limit $S$ to the nearest term.

If only a finite number of terms of the totally monotone sequence $(S_n)_n$ are known, then the situation is as follows :
<figure><img src="image_7.png" /></figure>

It means that if we want to get closer to the limit S for a totally monotonic sequence, we have to progress from top to bottom and from left to right in the ε table. The elements on the first descending diagonal appear to be the closest to the limit S.

Based on this observation, we can quote the following theorem ([Bre78], page 83) :

**Theorem 8.** Let $(S_n)_n$ be a totally monotone sequence which converges to S. Of all the elements in the ε table that can be calculated from $S_0, S_1, \dots, S_p$; the nearest to S is :

$$ \varepsilon_p^{(0)} \quad \text{if p is an even number} $$

$$ \varepsilon_{p-1}^{(1)} \quad \text{if p is an odd number} $$

Let us now study the convergence of the ε-algorithm for totally oscillating sequences.
Let us give a definition of a totally oscillating sequence :

**Definition 8.** The sequence $(S_n)_n$ is said to be totally oscillating if the sequence $((-1)^n S_n)_n$ is totally monotonic.

For these types of sequences, we have the following convergence result :

**Theorem 9.** If we apply the ε-algorithm to a sequence $(S_n)_n$ which converges to S and if there exists two constants $a \neq 0$ and $b$ such that the sequence $(aS_n + b)_n$ is totally oscillating, then :

$$ \lim_{n \to \infty} \varepsilon_{2k}^{(n)} = S \quad k = 0, 1, \dots \text{fixed} $$

$$ \lim_{k \to \infty} \varepsilon_{2k}^{(n)} = S \quad n = 0, 1, \dots \text{fixed.} $$

*Proof.* Available in [Bre77] page 77 (Theorem 51) and page 80 (Theorem 51 bis). □

**Remark 10.** In the particular case where $a = 1$ and $b = 0$, it means that the sequence $(S_n)_n$ is completely oscillating and that the ε-algorithm applied to this sequence converges in n and k to the limit S. It converges in n for even columns. It converges in k along the descending diagonals, in keeping only the even lower indices.

To know how we must progress in the ε table to accelerate the convergence of totally oscillating sequence, the following theorem states :
**Theorem 10.** If we apply the ε-algorithm to a totally oscillating sequence $(S_n)_n$, then, for $n, k = 0, 1, \dots$, we have:

$$ 
\begin{aligned}
0 &\le (-1)^n \varepsilon_{2k+2}^{(n)} \le (-1)^n \varepsilon_{2k}^{(n)} \\
(-1)^n \Delta \varepsilon_{2k}^{(n)} &\le (-1)^n \Delta \varepsilon_{2k+2}^{(n)} \le 0 \\
0 &\le (-1)^n \varepsilon_{2k+2}^{(n)} \le (-1)^{n+2} \varepsilon_{2k}^{(n+2)}
\end{aligned}
$$

where $\Delta \varepsilon_{2p}^{(q)} = \varepsilon_{2p}^{(q+1)} - \varepsilon_{2p}^{(q)}$.

*Proof.* Available in [Bre77] page 78 (Lemma 17) and page 80 (Lemma 17 bis). □

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

The situation can be interpreted with the following diagram where the arrows go from the furthest term of the limit to the closest term of the limit $S$:

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

There is no connection between quantities with an even upper index and those with an odd upper index.

As for the totally monotone sequences, if we want to get closer to the limit $S$ for a totally monotonic sequence, we have to progress from top to bottom and from left to right in the $\varepsilon$ table.

According to [Bre78] (page 85), we can show that if $(S_n)_n$ is totally monotonic or totally oscillating, then

$$ \lim_{n \to \infty} \frac{S_{n+1} - S}{S_n - S} = a \qquad (4) $$

with $0 < a \le 1$ for totally monotonic sequences and $-1 \le a < 1$ for totally oscillating sequences.

Theorems 6 and 9 are convergence theorems, not convergence acceleration theorems. However, we find in [Bre77] (pages 80-81):

**Theorem 11.** If we apply the ε-algorithm on a totally oscillating or monotonic non-logarithmic sequence ((4) with $a < 1$), then:

$$ \lim_{n \to \infty} \frac{\varepsilon_{2k}^{(n)} - S}{S_{n+2k} - S} = 0 \quad \text{for } k = 1, 2, \dots \text{ fixed} $$

and,

$$ \lim_{k \to \infty} \frac{\varepsilon_{2k}^{(n)} - S}{S_{n+2k} - S} = 0 \quad \text{for } n = 1, 2, \dots \text{ fixed.} $$
It means that we have acceleration of the convergence of monotonic or oscillating sequences when we apply the ε-algorithm on them. Moreover, as seen in theorems 7 and 10, the further down on the right we go, the closer we are to the limit of the sequence $(S_n)_n$.

Based on this observation, the accelerated sequence generated by our algorithm will be the quantities located on the first downward diagonal with an even lower index $(\epsilon_{2k}^{(0)})_k$.

### 2.1.3 Implementation of the algorithm

#### Presentation

We want to construct the sequence $(\epsilon_{2k}^{(0)})_k$ for $k = 0, 1, ...$

The algorithm implemented is based on the expression (1) of the ε-algorithm, which we recall here:

$$ 
\begin{array}{l@{\qquad}l@{\qquad}l}
\epsilon_{-1}^{(n)} = 0 & \epsilon_{0}^{(n)} = S_n & n = 0, 1, \dots \\
\epsilon_{k+1}^{(n)} = \epsilon_{k-1}^{(n+1)} + \frac{1}{\epsilon_{k}^{(n+1)} - \epsilon_{k}^{(n)}} & n, k = 0, 1, \dots &
\end{array}
$$

As seen in the remark 3, when a new element $S_k$ ($k \in \mathbb{N}$) of the initial sequence is known, an ascending diagonal of the ε table can be computed until the element $\epsilon_k^{(0)}$ is reached. If $k$ is an even number, $\epsilon_k^{(0)}$ is an element of the accelerated sequence.

To compute this ascending diagonal in using the above recurrence, the previous ascending diagonal must be kept in memory. Let us take an example with the ε table. If we want to compute the new element represented by a red point, we have to know elements represented by green ones:

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

When we reach an element $\epsilon_k^{(0)}$ with $k$ an even number, we check if this element of the accelerated sequence is enough close to the limit or not. For that, we need to define a stopping criterion and a tolerance associated (denoted by tol). If the limit $S$ is known, the stopping criterion will be

$$ |\epsilon_k^{(0)} - S| < tol. $$
If the limit $S$ isn't known, the stopping criterion will be the absolute error between two consecutive terms of the accelerated sequence

$$|\varepsilon_k^{(0)} - \varepsilon_{k-2}^{(0)}| < tol$$

or the relative error

$$\frac{|\varepsilon_k^{(0)} - \varepsilon_{k-2}^{(0)}|}{|\varepsilon_{k-2}^{(0)}|} < tol.$$

If the stopping criterion isn't satisfied, we compute another term of the initial sequence $(S_n)_n$ and construct an ascending diagonal until the stopping criterion is verified for the accelerated sequence or the initial sequence.

The main structure of the code is known. However, some cancellation issues can appear with the formulation (1) of the $\varepsilon$-algorithm.

## Cancellation issues

The system IEEE754 is a standard for the representation of floating-point numbers in binary. This representation only uses a limited number of bits. For example, the representation of double precision numbers is made with 64 bits. There are 1 bit for the sign, 11 bits for the exponent, and 52 bits for the mantissa:

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

Because of the limited number of bits, doing operations with floating-point numbers can reveal truncation errors. So, a series of computation can lead to important computation errors, due to the accumulation of rounding errors.

One of the most well-known source of errors is the cancellation phenomenon. This phenomenon appears when we subtract two terms that are very close to each other. Indeed, their bits representation will be similarly the same. So, when we'll subtract these two terms, there will be hardly any significant bits left. The result will be small compared to the two operands. Let us see an example to highlight the cancellation phenomenon.

**Example 2.** Let us take a variable $\delta$. We want to compute the quantity $\frac{(1+\delta)-1}{\delta}$.

Analytically, for any $\delta \neq 0$, the previous expression is equal to 1. Now, let us compute on the laptop this quantity in double precision by diminishing $\delta$:

<table><thead><tr><th>Iteration</th><th>δ</th><th>(1+δ)-1/δ</th><th>Exact result</th></tr></thead><tbody><tr><td>1</td><td>1.23456789101112130e+00</td><td>9.99999999999999778e-01</td><td>1.0e+00</td></tr><tr><td>2</td><td>1.23456789101112130e-01</td><td>9.99999999999999112e-01</td><td>1.0e+00</td></tr><tr><td>3</td><td>1.23456789101112144e-02</td><td>1.00000000000000622e+00</td><td>1.0e+00</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>14</td><td>1.23456789101112205e-13</td><td>1.00000009102136289e+00</td><td>1.0e+00</td></tr><tr><td>15</td><td>1.23456789101112208e-14</td><td>1.00719433628050936e+00</td><td>1.0e+00</td></tr><tr><td>16</td><td>1.23456789101112212e-15</td><td>1.07913678887197428e+00</td><td>1.0e+00</td></tr><tr><td>17</td><td>1.23456789101112227e-16</td><td>1.79856131478662351e+00</td><td>1.0e+00</td></tr><tr><td>18</td><td>1.23456789101112227e-17</td><td>0.00000000000000000e+00</td><td>1.0e+00</td></tr></tbody></table>
Smaller the $\delta$ is, closer quantities $1 + \delta$ and $1$ are. Then, truncation errors appear with the cancellation phenomenon when the $\delta$ diminishes. For a small $\delta$, the results obtained are far from the analytic result. Let us remark that the quantity computed is equal to zero when the $\delta$ is smaller than the unit round-off in double precision ($\approx 1.11 \times 10^{-16}$). Indeed, in this case, $\delta$ is rounded to zero.

In the expression (1), the cancellation phenomenon can appear in the subtraction $\varepsilon_k^{(n+1)} - \varepsilon_k^{(n)}$, for $n$ and $k$ fixed. However, we compute the quotient of this difference. So, if the terms $\varepsilon_k^{(n+1)}$ and $\varepsilon_k^{(n)}$ are very closed, the previous subtraction will give a small number, and the quotient of the difference will give a big number. The cancellation phenomenon and rounding errors will be propagated at the level of these large values. Then, we will lose a lot of precision in the computation of $\varepsilon_{k-1}^{(n+1)}$ and the following terms in the $\varepsilon$-algorithm will be very imprecise.

To minimise the rounding errors, the computation of the $\varepsilon$-algorithm terms will be done in quadruple precision [Quad]. Here is the comparison between the floating-point numbers in single, double and quadruple precision :

<table><thead><tr><th>Type</th><th>Size</th><th>Significand</th><th>Exponent</th><th>Unit Round-off</th><th>Range</th></tr></thead><tbody><tr><td>Single</td><td>32 bits</td><td>23+1 bits</td><td>8 bits</td><td>2<sup>-24</sup> &approx; 5.96 &times; 10<sup>-8</sup></td><td>10<sup>&pm;38</sup></td></tr><tr><td>Double</td><td>64 bits</td><td>52+1 bits</td><td>11 bits</td><td>2<sup>-53</sup> &approx; 1.11 &times; 10<sup>-16</sup></td><td>10<sup>&pm;308</sup></td></tr><tr><td>Quadruple</td><td>128 bits</td><td>112+1 bits</td><td>15 bits</td><td>2<sup>-113</sup> &approx; 9.63 &times; 10<sup>-35</sup></td><td>10<sup>&pm;4932</sup></td></tr></tbody></table>

The advantage of working in quadruple precision is to keep enough information even if the two terms $\varepsilon_k^{(n+1)}$ and $\varepsilon_k^{(n)}$ are very closed. Indeed, as we see in the above array, zero machine is reached around $1.0 \times 10^{-34}$. The math library `libquadmath` provides this quadruple precision for coding in C and C++.

Let us go back to the previous example and look at the impact of quadruple precision on truncation errors.

**Example 3.** Let us compute the quantity $\frac{(1+\delta)-1}{\delta}$ in quadruple precision. When $\delta$ decreases, we obtain:

<table><thead><tr><th>Iteration</th><th>δ</th><th>(1+δ)-1/δ</th><th>Exact result</th></tr></thead><tbody><tr><td>1</td><td>1.23456789101112130e+00</td><td>1.00000000000000000000e+00</td><td>1.0e+00</td></tr><tr><td>2</td><td>1.23456789101112130e-01</td><td>1.00000000000000000000e+00</td><td>1.0e+00</td></tr><tr><td>3</td><td>1.23456789101112144e-02</td><td>1.00000000000000000000e+00</td><td>1.0e+00</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>20</td><td>1.23456789101112249e-19</td><td>1.00000000000000000000e+00</td><td>1.0e+00</td></tr><tr><td>21</td><td>1.23456789101112264e-20</td><td>9.999999999999996891e-01</td><td>1.0e+00</td></tr><tr><td>22</td><td>1.23456789101112283e-21</td><td>9.999999999999996891e-01</td><td>1.0e+00</td></tr><tr><td>23</td><td>1.23456789101112287e-22</td><td>1.0000000000000077671e+00</td><td>1.0e+00</td></tr><tr><td>24</td><td>1.23456789101112293e-23</td><td>9.999999999992976729e-01</td><td>1.0e+00</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>32</td><td>1.23456789101112335e-31</td><td>9.99962094705972859e-01</td><td>1.0e+00</td></tr><tr><td>33</td><td>1.23456789101112357e-32</td><td>9.98402091438100769e-01</td><td>1.0e+00</td></tr><tr><td>34</td><td>1.23456789101112354e-33</td><td>9.36001960723219395e-01</td><td>1.0e+00</td></tr><tr><td>35</td><td>1.23456789101112362e-34</td><td>1.56000326787203236e+00</td><td>1.0e+00</td></tr><tr><td>36</td><td>1.23456789101112357e-35</td><td>0.00000000000000000000e+00</td><td>1.0e+00</td></tr></tbody></table>
In this example, the truncation errors seems to appear when we reach a $\delta$ closed to $10^{-20}$. When we worked in double precision, the quantity computed was equal to zero at this stage. As with the double precision, when the $\delta$ decreases, the truncation errors increase. But, the quantity computed is equal to the analytic one until the variable $\delta$ is near to $10^{-20}$. Let us remark that the quantity computed is equal to zero when $\delta$ is smaller than the unit round-off in quadruple precision ($\approx 9.63 \times 10^{-35}$).

We have just described the method we use to avoid the cancellation phenomenon as much as possible. However, P. Wynn has established a particular rule to avoid the cancellation phenomenon when $\varepsilon_k^{(n+1)}$ and $\varepsilon_k^{(n)}$ are very closed or equal ([Bre78], pages 316-317). To explain the particular rule, let us place ourselves in the following context.

Let us suppose $\varepsilon_{k-2}^{(n+1)}$ and $\varepsilon_{k-2}^{(n+2)}$ both become equal to the same quantity $b$. Then, $\varepsilon_{k-1}^{(n+1)}$ becomes infinite (because $\varepsilon_{k-1}^{(n+1)} = \varepsilon_{k-3}^{(n+2)} + \frac{1}{\varepsilon_{k-2}^{(n+2)} - \varepsilon_{k-2}^{(n+1)}}$, $\varepsilon_k^{(n)}$ and $\varepsilon_k^{(n+1)}$ become equal to the quantity $b$, and $\varepsilon_{k+1}^{(n)}$ is indeterminate :

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

To overcome these errors of cancellation and division by zero, Wynn established the following particular rule when the previous situation happens. We have to compute the relation

$$ \varepsilon_{k+1}^{(n)} = \frac{a}{\left(1 + \frac{a}{\varepsilon_{k-1}^{(n+1)}}\right)} $$

where $a$ is defined by
$$a = \frac{\varepsilon_{k-1}^{(n+2)}}{\left(1 - \frac{\varepsilon_{k-1}^{(n+2)}}{\varepsilon_{k-1}^{(n+1)}}\right)} + \frac{\varepsilon_{k-1}^{(n)}}{\left(1 - \frac{\varepsilon_{k-1}^{(n)}}{\varepsilon_{k-1}^{(n+1)}}\right)} - \frac{\varepsilon_{k-3}^{(n+2)}}{\left(1 - \frac{\varepsilon_{k-3}^{(n+2)}}{\varepsilon_{k-3}^{(n+1)}}\right)}$$

This formulation allows us to ignore the unfortunate terms $\varepsilon_{k-2}^{(n+1)}$ and $\varepsilon_{k-2}^{(n+2)}$. But, we have to keep in memory terms $\varepsilon_{k-1}^{(n+1)}$, $\varepsilon_{k-1}^{(n+2)}$, $\varepsilon_{k-1}^{(n)}$ and $\varepsilon_{k-3}^{(n+2)}$ to compute $\varepsilon_{k+1}^{(n)}$. It is possible to show that this particular rule can be reduced to the rule of the cross given by relations (2) or (3) (Appendix [B]).

**Remark 11.** We decide to implement the algorithm with the quadruple precision, and without the particular rule.

**Remark 12.** When we check if the stopping criterion is correct for the accelerated sequence, we cast the term $\varepsilon_k^{(n)}$ in double precision. So, the stopping criterion for the initial and accelerated sequence is computed in double precision.

### 2.1.4 Application to iterative methods

In this section, we will investigate the performance of the scalar $\varepsilon$-algorithm by applying it to the following iterative methods: the fixed point method and Newton's method.

#### Fixed point method

The goal of fixed point problems is to find a real number $x \in I$ ($I \subset \mathbb{R}$) such that, for $f: I \to \mathbb{R}$,

$$f(x) = x.$$

A method for solving a fixed point problem is to construct a sequence generated by

$$S_0 \text{ given,} \\ S_{n+1} = f(S_n), \quad \forall n \in \mathbb{N},$$

until the sequence $(S_n)_n$ converges. The stopping criterion is $|f(S_n) - S_n| \le \delta$ with $\delta$ a small positive real number.

The construction of fixed-point iterations has the advantage of being simple and easy to implement. However, the convergence of the sequence generated to the solution isn't guaranteed (it depends on some properties on the function $f$) and is often linear with a slow convergence speed.

That's why we want to see if we can accelerate the convergence of such sequence with the $\varepsilon$-algorithm.

**First example:** Let us take the function $f(x) = \frac{\sin(\pi x)}{2} e^{\frac{x}{2}}$. There are different solutions to the fixed point problem which are $x_1 = -0.38611395805...$, $x_2 = 0$ and $x_3 = 0.62920719609...$.
For the initial point, let us put $S_0 = 0.3$. We expect that the fixed point method converges to $S = x_3$.

The following results are obtained :

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

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

The convergence of the initial and accelerated sequences are also available in the convergence table below :

<table><caption>Convergence table for the initial sequence and the accelerated sequence</caption><thead><tr><th>n</th><th>S<sub>n</sub></th><th>Error : |S<sub>n</sub> - f(S<sub>n</sub>)|</th><th>ε-algorithm : ε<sub>n</sub><sup>(0)</sup></th><th>Error : |ε<sub>n</sub><sup>(0)</sup> - f(ε<sub>n</sub><sup>(0)</sup>)|</th></tr></thead><tbody><tr><td>0</td><td>3.000000000000000e-01</td><td>1.699718235069e-01</td><td>3.000000000000000e-01</td><td>1.699718235069e-01</td></tr><tr><td>2</td><td>6.296333878767816e-01</td><td>6.548130983241e-04</td><td>3.102104234431190e+00</td><td>3.845643016906e+00</td></tr><tr><td>4</td><td>6.293293183707927e-01</td><td>1.875039580524e-04</td><td>6.293264217362983e-01</td><td>1.830553348863e-04</td></tr><tr><td>6</td><td>6.292421577536804e-01</td><td>5.366887496338e-05</td><td>6.292071781446587e-01</td><td>2.755194783432e-08</td></tr><tr><td>8</td><td>6.292172024493184e-01</td><td>1.535967387544e-05</td><td>6.292071960975309e-01</td><td>4.944489262470e-12</td></tr><tr><td>10</td><td>6.292100597927642e-01</td><td>4.395683337321e-06</td><td>6.292071960943095e-01</td><td>2.220446049250e-16</td></tr></tbody></table>

The accelerated sequence generated with the ε-algorithm converges to $S = x_3$ and faster than the initial one. Indeed, 47 iterations are needed to reach the tolerance $(5 \times 10^{-16})$ with the sequence of the fixed point method. But, it only needs 10 terms of the previous sequence to reach the tolerance with the ε-algorithm. So, it divides the number of terms calculated with the fixed point method by almost 5.

Let us take another example which requires more iterations to converge.

Second example : We want to resolve the fixed point problem with $f(x) = -4.2x^3 + 5.8x^2 - 1.8x + 0.2$.

There is only one solution which is $S \approx 8.5697085 \times 10^{-02}$. We take for the initial point $S_0 = 0.8$.

<table><caption>Convergence table for the initial sequence and the accelerated sequence</caption><thead><tr><th>n</th><th>S<sub>n</sub></th><th>Error : |S<sub>n</sub> - f(S<sub>n</sub>)|</th><th>ε<sub>n</sub><sup>(0)</sup></th><th>Error : |ε<sub>n</sub><sup>(0)</sup> - f(ε<sub>n</sub><sup>(0)</sup>)|</th></tr></thead><tbody><tr><td>0</td><td>8.000000000000000e-01</td><td>4.784000000000e-01</td><td>8.000000000000000e-01</td><td>4.784000000000e-01</td></tr><tr><td>2</td><td>8.129372487679925e-02</td><td>8.451432447825e-03</td><td>-1.612456612135621e-01</td><td>8.198969350766e-01</td></tr><tr><td>4</td><td>8.213717340793169e-02</td><td>6.818316226940e-03</td><td>8.281633340325150e-02</td><td>5.508229790718e-03</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>20</td><td>8.505112400731149e-02</td><td>1.228294336187e-03</td><td>8.569708510823789e-02</td><td>2.289834988289e-15</td></tr><tr><td>22</td><td>8.517548877698161e-02</td><td>9.915083507983e-04</td><td>8.569708510823796e-02</td><td>2.400857290751e-15</td></tr><tr><td>24</td><td>8.527593542072517e-02</td><td>8.003683452830e-04</td><td>8.569708510823660e-02</td><td>1.804112415015e-16</td></tr></tbody></table>
<figure><img src="image_15.png" /></figure>

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

We can have the same conclusion than the previous example. However, in this case, the sequence generated by the fixed point method reached the tolerance after 285 iterations. We divide by almost 12 the number of terms computed with the ε-algorithm.

In the end, our two examples show that the ε-algorithm can speed up the convergence of sequences generated by fixed point method. The functions that have been taken allow us to have a fixed point method that converges linearly. However, what happens when we work with methods with high-order used for solving non-linear problems? Does the ε-algorithm speed up the sequences obtained by these methods?

### Newton's method

We use Newton's method when we want to find a real number $x \in I$ ($I \subset \mathbb{R}$) such that, for $f: I \to \mathbb{R}$,

$$f(x) = 0.$$

Newton's method is an iterative method given by

$$S_0 \text{ given,} \\ S_{n+1} = S_n - \frac{f(S_n)}{f'(S_n)}, \quad \forall n \in \mathbb{N},$$

until the sequence $(S_n)_n$ converges. The stopping criterion is $|f(S_n)| \le \delta$ with $\delta$ a small positive real number.

Newton's method has quadratic convergence (order 2) under certain assumptions, and this is its main advantage. However, it requires the derivative of the function $f$. It implies that the derivative exists.

Let us see if we can improve the order of convergence of a sequence generated by the Newton's method with the ε-algorithm; even if the convergence is already quadratic.

**First example:** The following example comes from [Vig], pages 38-39.

Let us take the function $f(x) = -4.1x^3 + 5.6x^2 - 1.7x + 0.2$ and the initial point $S_0 = 0$. The unique solution to $f(x) = 0$ is $x = 1$. We have:
<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /><figcaption>Convergence table for the initial sequence and the accelerated sequence</figcaption></figure>

<table><thead><tr><th>n</th><th>S<sub>n</sub></th><th>Error : |f(S<sub>n</sub>)|</th><th>ε<sub>n</sub><sup>(0)</sup></th><th>Error : |f(ε<sub>n</sub><sup>(0)</sup>)|</th></tr></thead><tbody><tr><td>0</td><td>0.000000000000000e+00</td><td>2.0000000000000e-01</td><td>0.000000000000000e+00</td><td>2.0000000000000e-01</td></tr><tr><td>2</td><td>2.458285756381450e-01</td><td>5.959995025545e-02</td><td>-1.313862607980254e+00</td><td>2.139942795158e+01</td></tr><tr><td>4</td><td>1.631009797981288e-01</td><td>5.391005695061e-02</td><td>1.408421474279540e-01</td><td>6.019815951406e-02</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>26</td><td>1.053601918301910e+00</td><td>1.699670105510e-01</td><td>1.375480190250052e-01</td><td>6.144776085674e-02</td></tr><tr><td>27</td><td>1.005772466291446e+00</td><td>1.638694739604e-02</td><td>//</td><td>//</td></tr><tr><td>28</td><td>1.000078126852174e+00</td><td>2.187960835388e-04</td><td>8.744677820497137e-02</td><td>9.142166481555e-02</td></tr><tr><td>29</td><td>1.000000014601469e+00</td><td>4.088411559299e-08</td><td>//</td><td>//</td></tr><tr><td>30</td><td>1.000000000000001e+00</td><td>1.498801083243e-15</td><td>9.845798021575279e-02</td><td>8.299444603691e-02</td></tr></tbody></table>

We can see that the original sequence generated by the Newton's method has a big search phase and that the convergence phase only begins at the 26th iteration. By definition, the original sequence fails to converge during the search phase. Therefore, accelerated sequence fails to do so. Indeed, the $\varepsilon$-algorithm works like a predictive algorithm. So, if the initial sequence doesn't converge, the $\varepsilon$-algorithm won't be able (except in few cases [Bre78], pages 48-49) to converge and accelerate it. Furthermore, the convergence phase seems to be too fast. The $\varepsilon$-algorithm fails to accelerate convergence. After 30 iterations, the accelerated sequence doesn't even show signs of convergence.

Let us take another example to see if the problem comes from the search phase.

Second example : Let us take the function $f(x) = \frac{\sin(\pi x)}{2} e^{\frac{x}{2}} - x$.

The initial point $S_0 = 1$ is taken and the solution we are supposed to converge to is $S = 0.62920719609...$.

We have,

<table><caption>Convergence table for the initial sequence and the accelerated sequence</caption><thead><tr><th>n</th><th>S<sub>n</sub></th><th>Error : |f(S<sub>n</sub>)|</th><th>ε<sub>n</sub><sup>(0)</sup></th><th>Error : |f(ε<sub>n</sub><sup>(0)</sup>)|</th></tr></thead><tbody><tr><td>0</td><td>1.000000000000000e+00</td><td>9.999999999999e-01</td><td>1.000000000000000e+00</td><td>9.999999999999e-01</td></tr><tr><td>1</td><td>7.214333614238075e-01</td><td>1.709019659994e-01</td><td>//</td><td>//</td></tr><tr><td>2</td><td>6.426297075549519e-01</td><td>2.122526111795e-02</td><td>6.115427878626955e-01</td><td>2.603889264412e-02</td></tr><tr><td>3</td><td>6.295897535420759e-01</td><td>5.877144614143e-04</td><td>//</td><td>//</td></tr><tr><td>4</td><td>6.292075245956608e-01</td><td>5.042359407791e-07</td><td>6.300570086829122e-01</td><td>1.306918594083e-03</td></tr><tr><td>5</td><td>6.292071960945523e-01</td><td>3.723688024592e-13</td><td>//</td><td>//</td></tr><tr><td>6</td><td>6.292071960943098e-01</td><td>1.110223024625e-16</td><td>6.292071823737325e-01</td><td>2.106050034367e-08</td></tr></tbody></table>
<figure><img src="image_19.png" /></figure>

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

The $\epsilon$-algorithm seems to converge to the solution. However, this convergence is slower than the original sequence. In this case, we can see that the $\epsilon$-algorithm does not speed up the convergence of the sequence generated by Newton's method. During the convergence phase, Newton's method converges in a quadratic way when the solution is a simple real root. It seems possible to say that the $\epsilon$-algorithm fails to accelerate the convergence of sequences of order greater than or equal to 2. Experimentally, we cannot do better than order 2 with the $\epsilon$-algorithm. However, it is possible to speed up the convergence of a Newton's method when the solution is a multiple real root of the function $f$. In this context, Newton's method does not converge in a quadratic way.

Let us take the function $f(x) = x^3 - 5x^2 + 8x - 4 = (x - 2)^2(x - 1)$ where 2 is a multiple real root. In taking the initial point $S_0 = 4$, we get:

<table><caption>Convergence table for the initial sequence and the accelerated sequence</caption><thead><tr><th>n</th><th>S<sub>n</sub></th><th>Error : |f(S<sub>n</sub>)|</th><th>ε<sub>n</sub><sup>(0)</sup></th><th>Error : |f(ε<sub>n</sub><sup>(0)</sup>)|</th></tr></thead><tbody><tr><td>0</td><td>4.000000000000000e+00</td><td>1.2000000000000e+01</td><td>4.000000000000000e+00</td><td>1.2000000000000e+01</td></tr><tr><td>2</td><td>2.760869565217391e+00</td><td>1.019407002547e+00</td><td>1.843749999999998e+00</td><td>2.059936523437e-02</td></tr><tr><td>4</td><td>2.254024574811046e+00</td><td>8.092030547366e-02</td><td>2.952364869803970e+00</td><td>1.770792682392e+00</td></tr><tr><td>6</td><td>2.073327528611925e+00</td><td>5.771203180609e-03</td><td>1.931371920890656e+00</td><td>4.386587806457e-03</td></tr><tr><td>8</td><td>2.019276742640324e+00</td><td>3.787559057233e-04</td><td>2.009569946753165e+00</td><td>9.246033372356e-05</td></tr><tr><td>10</td><td>2.004887650796044e+00</td><td>2.400589202800e-05</td><td>1.999731999853756e+00</td><td>7.180482519686e-08</td></tr><tr><td>12</td><td>2.001226371937800e+00</td><td>1.505832578629e-06</td><td>2.000001918540625e+00</td><td>3.677058657558e-12</td></tr><tr><td>14</td><td>2.000306874665135e+00</td><td>9.420095992140e-08</td><td>1.999999996690187e+00</td><td>1.776356839400e-15</td></tr><tr><td>16</td><td>2.000076736318657e+00</td><td>5.888912468776e-09</td><td>2.000000000039581e+00</td><td>1.776356839400e-15</td></tr><tr><td>18</td><td>2.000019185222217e+00</td><td>3.680771243352e-10</td><td>2.000000000174146e+00</td><td>3.552713678800e-15</td></tr><tr><td>20</td><td>2.000004796398688e+00</td><td>2.300737378391e-11</td><td>2.000000000244601e+00</td><td>0.0000000000000e+00</td></tr></tbody></table>

The Newton's method isn't of order 2 and the $\epsilon$-algorithm converges faster.

In the general case, the $\epsilon$-algorithm allows us to accelerate the convergence of sequences which converge linearly. However, when the order is superior or equal to two, we have difficulties in accelerating the convergence.
## 2.2 The vectorial ε-algorithm

In this section, we'll work on vector sequences. These sequences will be denoted by $(\mathbf{S}_n)_n$, where $n \in \mathbb{N}$. For each $n$, $\mathbf{S}_n \in \mathbb{R}^d$ with $d \ge 2$. If $(\mathbf{S}_n)_n$ converges, we'll denote by $\mathbf{S}$ its finite limit. Now, let us take another sequence $(\mathbf{V}_n)_n$ which also converges to $\mathbf{S}$. We say that $(\mathbf{V}_n)_n$ converges faster than $(\mathbf{S}_n)_n$ if

$$ \lim_{n \to \infty} \frac{\|\mathbf{V}_n - \mathbf{S}\|}{\|\mathbf{S}_n - \mathbf{S}\|} = 0. $$

By analogy with the section 2.2 on comparison between two sequences, in denoting by $(v_n)_n$ the sequence defined by $v_n = \|\mathbf{V}_n - \mathbf{S}\|$, $\forall n \in \mathbb{N}$, and $(u_n)_n$ the sequence defined by $u_n = \|\mathbf{S}_n - \mathbf{S}\|$, $\forall n \in \mathbb{N}$, we have $v_n = o(u_n)$.

### 2.2.1 Presentation

The vectorial form of the ε-algorithm has also been studied by P. Wynn. The formulation of this algorithm is based on the expression (1) of the scalar ε-algorithm and is adapted to work with vectors.
Let $(\mathbf{S}_n)_n$ a vector sequence with $\mathbf{S}_n \in \mathbb{R}^d$, for $n \in \mathbb{N}$.

The vectorial ε-algorithm looks like :

$$ \begin{gathered} \varepsilon_{-1}^{(n)} = \mathbf{0} \in \mathbb{R}^d \qquad \varepsilon_0^{(n)} = \mathbf{S}_n \in \mathbb{R}^d \qquad n = 0, 1, \dots \\ \varepsilon_{k+1}^{(n)} = \varepsilon_{k-1}^{(n+1)} + \left[ \varepsilon_k^{(n+1)} - \varepsilon_k^{(n)} \right]^{-1} \qquad n, k = 0, 1, \dots \end{gathered} \tag{5} $$

where, for $n$ and $k$ fixed, $\varepsilon_k^{(n)} \in \mathbb{R}^d$.

To be able to apply this algorithm, we have to define the inverse of a vector. Indeed, $\varepsilon_k^{(n+1)} - \varepsilon_k^{(n)}$ is a vector in $\mathbb{R}^d$. For that, P. Wynn defined the inverse of a vector $\mathbf{v} = (v_i)_{1 \le i \le d} \in \mathbb{R}^d$ by the relation

$$ \mathbf{v}^{-1} = \frac{\mathbf{v}}{(\mathbf{v}; \mathbf{v})} $$

with the dot product $(\mathbf{v}; \mathbf{v})$ defined by

$$ (\mathbf{v}; \mathbf{v}) = \sum_{i=1}^{d} v_i^2 = \|\mathbf{v}\|_2^2. $$

As this algorithm is based on the formulation of the scalar ε-algorithm, the vectorial ε-algorithm is presented in the same way.

Indeed, we can construct the same double-entry table called the ε table where the elements of the sequence $(\mathbf{S}_n)_n$ are located in the second column. However, the quantities $\varepsilon_k^{(n)}$ will be vectors. We progress from left to right and from top to bottom in the ε table. Moreover, for each new element of the initial sequence $(\mathbf{S}_n)_n$, we can construct an ascending diagonal on the ε table.
Still based on the results of the scalar ε-algorithm, only quantities with an even lower index can be integrated in the new sequence created by the ε-algorithm. We can also find a rule of the cross for the vectorial ε-algorithm which is :

$$ \begin{gathered} \varepsilon_{-2}^{(n)} = \infty \qquad \varepsilon_{0}^{(n)} = \mathbf{S}_{n} \qquad n = 0, 1, \dots \\ \frac{\varepsilon_{2k+2}^{(n-1)} - \varepsilon_{2k}^{(n)}}{\|\varepsilon_{2k+2}^{(n-1)} - \varepsilon_{2k}^{(n)}\|_{2}^{2}} + \frac{\varepsilon_{2k-2}^{(n+1)} - \varepsilon_{2k}^{(n)}}{\|\varepsilon_{2k-2}^{(n+1)} - \varepsilon_{2k}^{(n)}\|_{2}^{2}} = \frac{\varepsilon_{2k}^{(n+1)} - \varepsilon_{2k}^{(n)}}{\|\varepsilon_{2k}^{(n+1)} - \varepsilon_{2k}^{(n)}\|_{2}^{2}} + \frac{\varepsilon_{2k}^{(n-1)} - \varepsilon_{2k}^{(n)}}{\|\varepsilon_{2k}^{(n-1)} - \varepsilon_{2k}^{(n)}\|_{2}^{2}} \quad n, k = 0, 1, \dots \end{gathered} \tag{6} $$

### 2.2.2 Theory

In contrast to the scalar ε-algorithm, which is based on the Shanks transformation and on which a mathematical theory can be built, the vectorial ε-algorithm is poor in theoretical results. It was created from the formulation of the scalar ε-algorithm but adapted to vectors.

However, the following result, given in [Bre77] page 145, is closed to the theorem 2 :

**Theorem 12.** Let $\mathbf{S} \in \mathbb{R}$ and $(\mathbf{S}_n)_n$ a sequence of vectors of $\mathbb{R}^d$ such that :

$$ \sum_{i=0}^{k} a_i (\mathbf{S}_{n+i} - \mathbf{S}) = 0 \qquad n = 0, 1, \dots $$

where the $a_i$ are real numbers with $a_k \neq 0$.

If we apply the vectorial ε-algorithm to this sequence $(\mathbf{S}_n)_n$, then :

$$ \varepsilon_{2k}^{(n)} = \mathbf{S} \quad \text{for} \quad n = 0, 1, \dots \quad \text{if} \quad \sum_{i=0}^{k} a_i \neq 0 $$

and

$$ \varepsilon_{2k}^{(n)} = 0 \quad \text{for} \quad n = 0, 1, \dots \quad \text{if} \quad \sum_{i=0}^{k} a_i = 0. $$

It is not possible to express the vectors $\varepsilon_k^{(n)}$ of the expression (5) as generalised determinants. Then, many of the properties for the scalar ε-algorithm have no equivalent in the vectorial ε-algorithm ([Bre78], page 111). However, we still have the following result :

**Theorem 13.** Let a be an arbitrary non-zero real constant and **b** an arbitrary vector of $\mathbb{R}^d$. If the application of the vectorial ε-algorithm to the sequences $(\mathbf{S}_n)_n$ and $(a\mathbf{S}_n + \mathbf{b})_n$ respectively provides vectors noted $\varepsilon_k^{(n)}$ and $\tilde{\varepsilon}_k^{(n)}$, then :

$$ \tilde{\varepsilon}_{2k}^{(n)} = a\varepsilon_{2k}^{(n)} + \mathbf{b} \qquad n, k = 0, 1, \dots $$

$$ \tilde{\varepsilon}_{2k+1}^{(n)} = \frac{\varepsilon_{2k}^{(n)}}{a} \qquad n, k = 0, 1, \dots $$

*Proof.* Available in [Bre77], page 146. □
### 2.2.3 Implementation of the algorithm

We proceed in the same way as the scalar $\epsilon$-algorithm. Our objective is to construct the sequence $(\epsilon_{2k}^{(0)})_k$ for $k = 0, 1, ...$

The algorithm we have implemented is based on the relation (5) of the vectorial $\epsilon$-algorithm. When a new element $S_k$ ($k \in \mathbb{N}$) of the initial sequence is known, an ascending diagonal of the $\epsilon$ table can be computed until the element $\epsilon_k^{(0)}$ is reached. If $k$ is an even number, $\epsilon_k^{(0)}$ is an element of the accelerated sequence. Then, we check if this element is close enough to the limit or not with a stopping criterion and a tolerance on it. The computation of the error is made with the euclidean norm. If the limit $S$ is known, the stopping criterion is

$$ \|\epsilon_k^{(0)} - S\|_2 < tol. $$

If the limit $S$ isn't known, the stopping criterion will be the absolute error between two consecutive terms of the accelerated sequence

$$ \|\epsilon_k^{(0)} - \epsilon_{k-2}^{(0)}\|_2 < tol $$

or the relative error

$$ \frac{\|\epsilon_k^{(0)} - \epsilon_{k-2}^{(0)}\|_2}{\|\epsilon_{k-2}^{(0)}\|_2} < tol. $$

During the computation of a new element $\epsilon_{k+1}^{(n)}$ with $n$ and $k$ given, we first compute the dot product

$$ (\epsilon_k^{(n+1)} - \epsilon_k^{(n)}; \epsilon_k^{(n+1)} - \epsilon_k^{(n)}). $$

Then, we make the appropriate additions and subtractions component by component. Let us remark these computations can be done in parallel.

All these computations are done in quadruple precision.

**Remark 13.** Another method for accelerating the convergence of vector sequences would be to use the scalar $\epsilon$-algorithm component per component.

### 2.2.4 Application to the fixed point method

Based on the results obtained with the scalar $\epsilon$-algorithm, the vectorial $\epsilon$-algorithm will be only applied on the fixed point method.

The goal of a fixed point problem is to find a $d$-tuple $X \in \mathbb{R}^d$ ($d \ge 2$) such that, for $F : \mathbb{R}^d \to \mathbb{R}^d$,

$$ F(X) = X. $$

To find it, we construct a sequence of elements in $\mathbb{R}^d$ with the instructions
$\mathbf{S}_0$ given,

$$ \mathbf{S}_{n+1} = F(\mathbf{S}_n), \quad \forall n \in \mathbb{N}, $$

until the sequence $(\mathbf{S}_n)_n$ converges. The stopping criterion is $\|F(\mathbf{S}_n) - \mathbf{S}_n\|_2 \le \delta$ with $\delta$ a small positive real number.

We applied the vectorial $\varepsilon$-algorithm on two examples. One with a function $F: \mathbb{R}^2 \to \mathbb{R}^2$, and another one with $F: \mathbb{R}^3 \to \mathbb{R}^3$. For both, we worked with $\delta = 5.0 \times 10^{-16}$.

The two examples are taken from [Bur] (page 637).

**First example:** Let us take the function $F: \mathbb{R}^2 \to \mathbb{R}^2$ defined by

$$ F(\mathbf{X}) = \left( \begin{array}{c} \frac{x_1^2 + x_2^2 + 8}{10} \\ \frac{x_1 x_2^2 + x_1 + 8}{10} \end{array} \right). $$

The initial point is $\mathbf{S}_0 = (0.5; 0.25)$ and the sequence must converge to $\mathbf{S} = (1.0; 1.0)$.

We obtained:

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

<table><caption>Convergence table for the initial and accelerated sequences</caption><thead><tr><th>n</th><th>Error : |S<sub>n</sub> - F(S<sub>n</sub>)|</th><th>Error : |ε<sub>n</sub><sup>(0)</sup> - F(ε<sub>n</sub><sup>(0)</sup>)|</th></tr></thead><tbody><tr><td>0</td><td>6.8810342836306226e-01</td><td>6.8810342836306226e-01</td></tr><tr><td>2</td><td>4.9725203875172062e-02</td><td>3.6056763919160670e-02</td></tr><tr><td>4</td><td>7.4229251363374903e-03</td><td>5.1841593910959645e-04</td></tr><tr><td>6</td><td>1.1750746342499222e-03</td><td>1.7725331298734953e-06</td></tr><tr><td>8</td><td>1.8769352951098791e-04</td><td>3.5256169223706905e-09</td></tr><tr><td>10</td><td>3.0022829460351418e-05</td><td>1.3047359312753205e-12</td></tr><tr><td>12</td><td>4.8034445167189884e-06</td><td>7.0216669371534024e-16</td></tr><tr><td>14</td><td>7.6854579310693708e-07</td><td>0.0000000000000000e+00</td></tr></tbody></table>

The vectorial $\varepsilon$-algorithm accelerates the convergence of the sequence generated by the fixed point method. Indeed, it needs 37 iterations for the initial sequence to reach the tolerance $5.0 \times 10^{-16}$. But, it requires only 14 iterations with the vectorial $\varepsilon$-algorithm.
Second example : Let us take the function $F : \mathbb{R}^3 \to \mathbb{R}^3$ defined by

$$F(\mathbf{X}) = \left( \begin{array}{c} \frac{13 - x_2^2 + 4x_3}{15} \\ \frac{11 + x_3 - x_1^2}{10} \\ \frac{22 + x_2^3}{25} \end{array} \right).$$

The initial point is $S_0 = (9.0; 5.0; -2.0)$.

We obtained :

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

<table><caption>Convergence table for the initial and accelerated sequences</caption><thead><tr><th>n</th><th>Error : |S<sub>n</sub> - F(S<sub>n</sub>)|</th><th>Error : |ε<sub>n</sub><sup>(0)</sup> - F(ε<sub>n</sub><sup>(0)</sup>)|</th></tr></thead><tbody><tr><td>0</td><td>1.7824482538850258e+01</td><td>1.7824482538850258e+01</td></tr><tr><td>2</td><td>1.5321978712938131e+01</td><td>3.4789461935093007e+00</td></tr><tr><td>4</td><td>7.7892752286272426e-01</td><td>4.3538247150481340e+00</td></tr><tr><td>...</td><td>...</td><td>...</td></tr><tr><td>16</td><td>3.5819957529488642e-08</td><td>4.0418217762533309e-14</td></tr><tr><td>18</td><td>2.6143085098005391e-09</td><td>5.5511151231257827e-16</td></tr><tr><td>20</td><td>1.9091823617961324e-10</td><td>2.4825341532472731e-16</td></tr></tbody></table>

The same conclusion can be established for this example. We need 31 iterations to reach the tolerance with the initial sequence. With the vectorial ε-algorithm, it only needs 20 iterations.

**Remark 14.** In all the previous examples for the scalar and vectorial ε-algorithm, we only compared the number of iterations between the initial and accelerated sequences. Indeed, these are simple test cases. But, we also compare the execution time in the following section with large non-linear systems.

## 3 Application to concrete problems

All the problems which are described in this section can be resolved with codes available on the open source solver called *pugs*. It stands for Parallel Unstructured Grid Solvers. This project has
been launched in 2018. It is mainly, but not only, designed to deal with finite volume methods. *pugs* is a parallel software that uses *MPI* and multi-threading for parallelism. Multi-threading is achieved through an encapsulation of some *Kokkos* mechanisms. The philosophy of *pugs* is to provide "simple" numerical tools that are assembled together through a high-level language (a DSL close to the mathematics) to build more complex solvers. This approach is inspired by the success of FreeFEM, which uses a similar approach. Modules in *pugs* are written in C++, and then called via a .pgs file itself written in a more user friendly DSL language ([Ple], page 124).

Recall that all the numerical computations have been made on the following laptop [A].

## 3.1 Stationary diffusion problem

The study problem is a stationary diffusion problem [DDFV] written :

$$ \left\{ \begin{array}{ll} -\nabla \cdot (\kappa \nabla u) + \lambda u = f & \text{in } \Omega, \\ u = g & \text{on } \Gamma_D, \\ \kappa \nabla u \cdot \mathbf{n} = g & \text{on } \Gamma_N, \end{array} \right. \qquad (7) $$

where $\Omega$ is a bounded open domain of $\mathbb{R}^2$, with $\partial\Omega = \Gamma_D \cup \Gamma_N$ ($\Gamma_D \cap \Gamma_N = \emptyset$) and $\mathbf{n} \in \mathbb{R}^2$ the outgoing unit normal vector.

We have $f, \lambda \in L^2(\Omega)$, with $\lambda \ge 0$, and $g \in H^{1/2}(\partial\Omega)$.

The tensor-valued diffusion coefficient $\kappa$ is supposed to be bounded and to satisfy the following uniform ellipticity condition,

$$ \forall \mathbf{x} \in \Omega, \forall \mathbf{y} \in \mathbb{R}^2, \quad \alpha_{min} \|\mathbf{y}\|^2 \le \mathbf{y}^t \kappa(\mathbf{x}) \mathbf{y} \le \alpha_{max} \|\mathbf{y}\|^2, $$

where $\alpha_{min}$ and $\alpha_{max}$ are positive coefficients.

Under the above conditions, and if $\lambda > 0$ or $\Gamma_D$ is of positive length, system (7) has a unique solution in $H^1(\Omega)$. Such a solution satisfies a positiveness principle. It means if $f \ge 0$ and $g \ge 0$, then $u \ge 0$.

### 3.1.1 Problem 1 : $\kappa$ constant

Taking the domain $\Omega = ]0, 1]^2$, with the diffusion coefficient $\kappa$ a constant scalar function defined by $\kappa = 1$ and functions $f$ and $g$ defined by $f : \mathbb{R}^2 \to \mathbb{R} : (x, y) \mapsto 2\pi^2 \sin(\pi x) \sin(\pi y)$ and $g : \mathbb{R}^2 \to \mathbb{R} : (x, y) \mapsto 0$, we want to solve the preceding problem with $\lambda = 0$ and $\Gamma_N = \emptyset$:

$$ \left\{ \begin{array}{ll} -\nabla \cdot (\kappa \nabla u) = f & \text{in } \Omega, \\ u = g & \text{on } \Gamma_D = \partial\Omega. \end{array} \right. \qquad (8) $$

The exact solution of this problem is

$$ u_{ex} : \mathbb{R}^2 \to \mathbb{R} : (x, y) \mapsto \sin(\pi x) \sin(\pi y). $$

The method used to solve this problem is a finite volume method called DDFV (Discrete Duality Finite Volume), which is used for deformed meshes and which was already implemented
in pugs. The method is described in the article [DDFV]. This method is based on the construction of a primal and dual mesh. It returns a positive solution when the positiveness principle is verified in our problem.

We will resolve the previous problem with the DDFV method on two different meshes. The first one will be a deformed cartesian mesh and the second one will be a mesh generated randomly. Let us see how the problem is solved.

After discretizing the problem on the primal and dual mesh, we have to solve the following non-linear system

$$A(\mathbf{u})\mathbf{u} = \begin{pmatrix} A^{\text{primal}}(\mathbf{u}^{\text{primal}}, \mathbf{u}^{\text{dual}}) & \mathbf{0} \\ \mathbf{0} & A^{\text{dual}}(\mathbf{u}^{\text{primal}}, \mathbf{u}^{\text{dual}}) \end{pmatrix} \begin{pmatrix} \mathbf{u}^{\text{primal}} \\ \mathbf{u}^{\text{dual}} \end{pmatrix} = \begin{pmatrix} \mathbf{b}^{\text{primal}} \\ \mathbf{b}^{\text{dual}} \end{pmatrix} \quad (9)$$

with,

$$\mathbf{u}^{\text{primal}} = (u_i)_{1 \le i \le n}, \quad \mathbf{u}^{\text{dual}} = (u_r)_{1 \le r \le m}, \quad \mathbf{u} = (\mathbf{u}^{\text{primal}}, \mathbf{u}^{\text{dual}}),$$

$$\mathbf{b}^{\text{primal}} = (b_i)_{1 \le i \le n}, \quad \mathbf{b}^{\text{dual}} = (b_r)_{1 \le r \le m}, \quad \mathbf{b} = (\mathbf{b}^{\text{primal}}, \mathbf{b}^{\text{dual}}),$$

where $\mathbf{u}^{\text{primal}}$ is the solution on the primal mesh computed with the DDFV method, $n$ is the number of cells on the primal mesh, $\mathbf{u}^{\text{dual}}$ is the solution on the dual mesh computed by DDFV and $m$ is the number of cells on the dual mesh.

To solve this non-linear system, we use Picard's successive iteration method. Starting with a strictly positive initial function $u^0$ discretized on meshes, we construct the matrix $A(\mathbf{u}^0)$ and solve the system $A(\mathbf{u}^0)\mathbf{u}^1 = \mathbf{b}$. By repeating this process, we construct a vector sequence $(\mathbf{u}^k)_k$ which, if it converges to a positive vector, tends to the physical solution of the scheme. We stop the algorithm when the relative error between two iterates $\mathbf{u}^{k+1} - \mathbf{u}^k$ is sufficiently small. The algorithm can be summarised as follows:

<caption>**Algorithm 1** Algorithm solving (8)</caption>

Require: tol, $\mathbf{u}^0$
$k \leftarrow 0$
$\mathbf{u}^1 \leftarrow A^{-1}(\mathbf{u}^0)\mathbf{b}$
while ( tol $\|\mathbf{u}^k\|_2 < \|\mathbf{u}^{k+1} - \mathbf{u}^k\|_2$) do
$\mathbf{u}^{k+1} \leftarrow A^{-1}(\mathbf{u}^k)\mathbf{b}$
$k \leftarrow k + 1$
end while

Let us denote by $\| \cdot \|_2$ the $L^2$-norm weighted by the area of the cells and defined by $\|\mathbf{u}\|_2 = \sqrt{\sum_i V_i u_i^2}$ where $V_i$ is the cell area for the cell $i$ (for the primal or dual mesh).

The objective is to accelerate the convergence of the sequence $(\mathbf{u}^k)_k$ with the $\varepsilon$-algorithm.
The strategy adopted is to accelerate separately the primal sequence ((u<sup>primal</sup>)<sup>k</sup>)<sub>k</sub> and the dual sequence ((u<sup>dual</sup>)<sup>k</sup>)<sub>k</sub> obtained by solving the following non-linear systems

$$
\mathbf{A}^{\text{primal}} ((\mathbf{u}^{\text{primal}})^k, (\mathbf{u}^{\text{dual}})^k) (\mathbf{u}^{\text{primal}})^{k+1} = \mathbf{b}^{\text{primal}}
$$

$$
\mathbf{A}^{\text{dual}} ((\mathbf{u}^{\text{primal}})^k, (\mathbf{u}^{\text{dual}})^k) (\mathbf{u}^{\text{dual}})^{k+1} = \mathbf{b}^{\text{dual}}
$$

and we stop the algorithm when we have

$$
\|(\mathbf{u}^{\text{primal}})^{k+1} - (\mathbf{u}^{\text{primal}})^k\|_2 < \text{tol} \|(\mathbf{u}^{\text{primal}})^k\|_2
$$

and,

$$
\|(\mathbf{u}^{\text{dual}})^{k+1} - (\mathbf{u}^{\text{dual}})^k\|_2 < \text{tol} \|(\mathbf{u}^{\text{dual}})^k\|_2.
$$

Sequences ((u<sup>primal</sup>)<sup>k</sup>)<sub>k</sub> and ((u<sup>dual</sup>)<sup>k</sup>)<sub>k</sub> are accelerated with the vectorial ε-algorithm.

Denote that all the following computations have been done with the first iterate defined by

$$
u^0 : \mathbb{R}^2 \to \mathbb{R} : (x, y) \mapsto xy.
$$

First case : The primal mesh is a deformed cartesian mesh. The dual mesh is obtained from
the primal one. The deformation on the primal mesh is given by

$$
(x, y) \rightarrow (x + 0.1 \sin(2\pi x) \sin(2\pi y), y + 0.1 \sin(2\pi x) \sin(2\pi y)).
$$

<figure><img src="image_23.png" /><figcaption>Figure 1: Primal mesh 32 × 32 used</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Figure 2: Dual mesh 32 × 32 used</figcaption></figure>

Here are the results obtained by accelerating the convergence of the initial sequence using
the vectorial ε-algorithm, working in quadruple precision on the computation of the accelerated
sequence terms and in imposing a tolerance of 10⁻¹² on the computation of the error :
<figure><img src="image_25.png" /><figcaption>Figure 3: $u^0$ on the primal mesh 32 × 32</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Figure 4: $u_{ex}$ on the dual mesh 32 × 32</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>Figure 5: $u^k$ on the primal mesh 32 × 32</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>Figure 6: Error $|u_{ex} - u^k|$ on the primal mesh 32 × 32</figcaption></figure>

<table><thead><tr><th colspan="3">Number of iterations to reach the imposed tolerance</th></tr><tr><th>Mesh</th><th>Original sequence</th><th>Accelerated sequence</th></tr></thead><tbody><tr><td>32 × 32</td><td>75</td><td>48</td></tr><tr><td>64 × 64</td><td>84</td><td>44</td></tr><tr><td>128 × 128</td><td>89</td><td>44</td></tr><tr><td>256 × 256</td><td>92</td><td>46</td></tr><tr><td>512 × 512</td><td>96</td><td>54</td></tr></tbody></table>

<table><thead><tr><th colspan="4">Execution time of the while loop</th></tr><tr><th>Mesh</th><th>Original sequence</th><th>Accelerated sequence</th><th>Gain</th></tr></thead><tbody><tr><td>32 × 32</td><td>5.86 sec</td><td>3.96 sec</td><td>32.42%</td></tr><tr><td>64 × 64</td><td>27.60 sec</td><td>14.53 sec</td><td>47.35%</td></tr><tr><td>128 × 128</td><td>311 sec = 5 min 11 sec</td><td>158 sec = 2 min 38 sec</td><td>49.19%</td></tr><tr><td>256 × 256</td><td>3136 sec = 52 min 16 sec</td><td>1499 sec = 24 min 59 sec</td><td>52.20%</td></tr><tr><td>512 × 512</td><td>35 665 sec ≈ 9 h 54 min</td><td>19 224 sec ≈ 5 h 20 min</td><td>46.09%</td></tr></tbody></table>

The execution times in the table come from an average after running the code 20 times for
meshes 32, 64 and 128; 3 times for the mesh 256 and 1 time for the mesh 512.
For every meshes, the number of iterations needed to reach tolerance is reduced. Then, it has an impact on execution time. Indeed, the execution time is reduced for all meshes, from 32.42% to 52.20% for the best one.

Let us look at the time taken by the $\varepsilon$-algorithm in the computation of the accelerated sequence.

<table><thead><tr><th colspan="4">Execution time of the ε-algorithm in the while loop</th></tr><tr><th>Mesh</th><th>Accelerated sequence</th><th>ε-algorithm</th><th>Rate</th></tr></thead><tbody><tr><td>32 × 32</td><td>3.96 sec</td><td>0.61 sec</td><td>15.40%</td></tr><tr><td>64 × 64</td><td>14.53 sec</td><td>1.62 sec</td><td>11.15%</td></tr><tr><td>128 × 128</td><td>158 sec = 2 min 38 sec</td><td>13.67 sec</td><td>8.65%</td></tr><tr><td>256 × 256</td><td>1499 sec = 24 min 59 sec</td><td>35.90 sec</td><td>2.39%</td></tr><tr><td>512 × 512</td><td>19 224 sec ≈ 5 h 20 min</td><td>153.65 sec ≈ 2 min 33 sec</td><td>0.80%</td></tr></tbody></table>

Despite the use of the quadruple precision, the execution time of the $\varepsilon$-algorithm counts for nothing in all the computation of the accelerated sequence compared to the fixed point method. The most costly in terms of computation time must be the resolution of the non-linear system (9).

We show below the evolution of the relative error for primal and dual sequences, original and accelerated, with the mesh 512 × 512:

<table><thead><tr><th colspan="5">Evolution of the relative error <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding-bottom: 1px;">||u<sup>k+1</sup> - u<sup>k</sup>||<sub>2</sub></div><div style="padding-top: 1px;">||u<sup>k</sup>||<sub>2</sub></div></div> on the deformed mesh 512 × 512</th></tr><tr><th>k</th><th>Original primal seq.</th><th>Accelerated primal seq.</th><th>Original dual seq.</th><th>Accelerated dual seq.</th></tr></thead><tbody><tr><td>2</td><td>1.452996e-01</td><td>1.157419e+00</td><td>1.452990e-01</td><td>1.157416e+00</td></tr><tr><td>4</td><td>1.884952e-02</td><td>3.641065e-02</td><td>1.886333e-02</td><td>3.630004e-02</td></tr><tr><td>6</td><td>5.010165e-03</td><td>3.636998e-03</td><td>5.013034e-03</td><td>3.657097e-03</td></tr><tr><td>8</td><td>1.830138e-03</td><td>4.136477e-04</td><td>1.830770e-03</td><td>4.101603e-04</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>50</td><td>1.187475e-08</td><td>1.835707e-12</td><td>1.187429e-08</td><td>2.008249e-12</td></tr><tr><td>52</td><td>7.522371e-09</td><td>1.365785e-12</td><td>7.522082e-09</td><td>4.036462e-13</td></tr><tr><td>54</td><td>4.781693e-09</td><td>2.613116e-13</td><td>4.781498e-09</td><td>5.449388e-13</td></tr></tbody></table>

The graph associated with this table is available in [C].

The accelerated primal and dual sequences converge to a solution of the problem. Moreover, during the computation, we check that this numerical solution obtained with the $\varepsilon$-algorithm is positive.

Let us see if the accelerated primal and dual sequences converge to the exact solution $u_{ex}$. For that, let us compute the error between $(u^{primal})^k$ and $u_{ex,h}$, and between $(u^{dual})^k$ and $u_{ex,h}$. We denote by $u_{ex,h}$ the exact solution $u_{ex}$ discretized on the mesh. We obtained:
<table><caption>Evolution of the error $||u^k - u_{ex,h}||_2$ on the mesh 32 × 32</caption><thead><tr><th>k</th><th>Original primal seq.</th><th>Accelerated primal seq.</th><th>Original dual seq.</th><th>Accelerated dual seq.</th></tr></thead><tbody><tr><td>2</td><td>1.546730e-02</td><td>1.938168e-02</td><td>1.496902e-02</td><td>1.701920e-02</td></tr><tr><td>4</td><td>2.800747e-03</td><td>1.804444e-03</td><td>5.912512e-03</td><td>5.616156e-03</td></tr><tr><td>6</td><td>1.487380e-03</td><td>1.062689e-03</td><td>5.829869e-03</td><td>5.866729e-03</td></tr><tr><td>8</td><td>1.129700e-03</td><td>1.022527e-03</td><td>5.844580e-03</td><td>5.851490e-03</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>44</td><td>9.823983e-04</td><td>9.823985e-04</td><td>5.865792e-03</td><td>5.865792e-03</td></tr><tr><td>46</td><td>9.823984e-04</td><td>9.823985e-04</td><td>5.865792e-03</td><td>5.865792e-03</td></tr><tr><td>48</td><td>9.823984e-04</td><td>9.823985e-04</td><td>5.865792e-03</td><td>5.865792e-03</td></tr></tbody></table>

<table><caption>Evolution of the error $||u^k - u_{ex,h}||_2$ on the mesh 128 × 128</caption><thead><tr><th>k</th><th>Original primal seq.</th><th>Accelerated primal seq.</th><th>Original dual seq.</th><th>Accelerated dual seq.</th></tr></thead><tbody><tr><td>2</td><td>1.520472e-02</td><td>1.934173e-02</td><td>1.510301e-02</td><td>1.901317e-02</td></tr><tr><td>4</td><td>2.703203e-03</td><td>1.867681e-03</td><td>2.680005e-03</td><td>1.933023e-03</td></tr><tr><td>6</td><td>9.620485e-04</td><td>3.085660e-04</td><td>1.122782e-03</td><td>7.508588e-04</td></tr><tr><td>8</td><td>3.882245e-04</td><td>1.475394e-04</td><td>7.635105e-04</td><td>6.975487e-04</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>40</td><td>6.171682e-05</td><td>6.171819e-05</td><td>6.976548e-04</td><td>6.976549e-04</td></tr><tr><td>42</td><td>6.171720e-05</td><td>6.171819e-05</td><td>6.976549e-04</td><td>6.976549e-04</td></tr><tr><td>44</td><td>6.171751e-05</td><td>6.171819e-05</td><td>6.976549e-04</td><td>6.976549e-04</td></tr></tbody></table>

<table><caption>Evolution of the error $||u^k - u_{ex,h}||_2$ on the mesh 512 × 512</caption><thead><tr><th>k</th><th>Original primal seq.</th><th>Accelerated primal seq.</th><th>Original dual seq.</th><th>Accelerated dual seq.</th></tr></thead><tbody><tr><td>2</td><td>1.511144e-02</td><td>1.937977e-02</td><td>1.510168e-02</td><td>1.932448e-02</td></tr><tr><td>4</td><td>2.788547e-03</td><td>1.919393e-03</td><td>2.776707e-03</td><td>1.924043e-03</td></tr><tr><td>6</td><td>9.839635e-04</td><td>2.930143e-04</td><td>9.812493e-04</td><td>3.023657e-04</td></tr><tr><td>8</td><td>3.844514e-04</td><td>1.125941e-04</td><td>3.891549e-04</td><td>1.374690e-04</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>50</td><td>3.858532e-06</td><td>3.858930e-06</td><td>8.596670e-05</td><td>8.596669e-05</td></tr><tr><td>52</td><td>3.858665e-06</td><td>3.858930e-06</td><td>8.596670e-05</td><td>8.596669e-05</td></tr><tr><td>54</td><td>3.858754e-06</td><td>3.858930e-06</td><td>8.596670e-05</td><td>8.596669e-05</td></tr></tbody></table>

When the size of our mesh decreases, the numerical solution converges to the exact solution, for the original and accelerated sequences. In $L^2$-norm, the DDFV scheme is of order 2. This order of convergence is preserved with the accelerated sequence:

**Second case:** The primal mesh is a random mesh. To be able to compare the execution times between the convergence of the initial and accelerated sequences, we need to perform the computations on the same randomly generated mesh. To do it, we fix the seed as below:

The interest in working with a random mesh is that solving the non-linear system becomes more complex and takes longer time.

Here are the results obtained by accelerating the convergence of the initial sequence using the vectorial $\epsilon$-algorithm, working in quadruple precision on the computation of the accelerated sequence terms and in imposing a tolerance of $10^{-12}$ on the computation of the error:
<figure><img src="image_29.png" /><figcaption>L² errors for the accelerated primal sequence on the deformed mesh with the problem 1</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>Figure 7: Primal mesh 32 × 32 used</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>Figure 8: Dual mesh 32 × 32 used</figcaption></figure>

<table><caption>Number of iterations to reach the imposed tolerance</caption><thead><tr><td>Mesh</td><td>Original sequence</td><td>Accelerated sequence</td></tr></thead><tbody><tr><td>32 × 32</td><td>66</td><td>38</td></tr><tr><td>64 × 64</td><td>118</td><td>54</td></tr><tr><td>128 × 128</td><td>120</td><td>62</td></tr><tr><td>256 × 256</td><td>135</td><td>92</td></tr><tr><td>512 × 512</td><td>199</td><td>116</td></tr></tbody></table>

<table><caption>Execution time of the while loop</caption><thead><tr><td>Mesh</td><td>Original sequence</td><td>Accelerated sequence</td><td>Gain</td></tr></thead><tbody><tr><td>32 × 32</td><td>4.28 sec</td><td>3.30 sec</td><td>22.90%</td></tr><tr><td>64 × 64</td><td>39.03 sec</td><td>18.67 sec</td><td>52.17%</td></tr><tr><td>128 × 128</td><td>385 sec = 6 min 25 sec</td><td>215 sec = 3 min 35 sec</td><td>44.16%</td></tr><tr><td>256 × 256</td><td>4430 sec = 1 h 13 min 50 sec</td><td>3125 sec = 52 min 5 sec</td><td>29.45%</td></tr><tr><td>512 × 512</td><td>76 331 sec ≈ 21 h 12 min</td><td>46 505 sec ≈ 12 h 55 min</td><td>39.07%</td></tr></tbody></table>
The execution times in the table comes from an average after running the code 20 times for meshes 32 and 64; 10 times for the mesh 128; 2 times for the mesh 256 and 1 time for the mesh 512.

For every meshes, the number of iterations needed to reach tolerance is reduced by the $\varepsilon$-algorithm. Then, it has an impact on execution time. Indeed, the execution time is reduced for all meshes, from 22.90% to 52.17% for the best one.

Let us look at the time taken by the $\varepsilon$-algorithm in the computation of the accelerated sequence.

<table><thead><tr><th colspan="4">Execution time of the ε-algorithm in the while loop</th></tr><tr><th>Mesh</th><th>Accelerated sequence</th><th>ε-algorithm</th><th>Rate</th></tr></thead><tbody><tr><td>32 × 32</td><td>3.30 sec</td><td>0.42 sec</td><td>12.73%</td></tr><tr><td>64 × 64</td><td>18.67 sec</td><td>2.42 sec</td><td>12.96%</td></tr><tr><td>128 × 128</td><td>215 sec = 3 min 35 sec</td><td>23.47 sec</td><td>10.92%</td></tr><tr><td>256 × 256</td><td>3125 sec = 52 min 5 sec</td><td>124 sec ≈ 2 min 4 sec</td><td>3.97%</td></tr><tr><td>512 × 512</td><td>46 505 sec ≈ 12 h 55 min</td><td>681 sec ≈ 11 min 21 sec</td><td>1.46%</td></tr></tbody></table>

We have the same conclusion than the previous case. Despite the use of the quadruple precision, the execution time of the $\varepsilon$-algorithm counts for nothing in all the computation of the accelerated sequence compared to the fixed point method.

We show below the evolution of the relative error for primal and dual sequences, original and accelerated, with the mesh 512 × 512:

<table><thead><tr><th colspan="5">Evolution of the relative error $\frac{\|u^{k+1} - u^k\|_2}{\|u^k\|_2}$ on the random mesh 512 × 512</th></tr><tr><th>k</th><th>Original primal seq.</th><th>Accelerated primal seq.</th><th>Original dual seq.</th><th>Accelerated dual seq.</th></tr></thead><tbody><tr><td>2</td><td>3.166227e-02</td><td>1.151446e+00</td><td>3.398303e-02</td><td>1.177142e+00</td></tr><tr><td>4</td><td>1.639473e-03</td><td>3.351557e-02</td><td>8.492141e-04</td><td>1.466143e-02</td></tr><tr><td>6</td><td>4.812553e-04</td><td>3.469474e-03</td><td>2.893217e-04</td><td>2.442970e-04</td></tr><tr><td>8</td><td>1.829692e-04</td><td>6.438433e-04</td><td>9.822530e-05</td><td>4.297215e-04</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>112</td><td>4.963164e-10</td><td>2.453042e-12</td><td>3.832522e-10</td><td>6.858652e-13</td></tr><tr><td>114</td><td>4.297076e-10</td><td>1.382528e-12</td><td>3.319570e-10</td><td>9.053758e-13</td></tr><tr><td>116</td><td>3.720193e-10</td><td>8.139259e-13</td><td>2.875861e-10</td><td>7.316865e-13</td></tr></tbody></table>

The graph associated with this table is available in [C].

The accelerated primal and dual sequences converge to a solution of the problem. Moreover, during the computation, we check that this numerical solution obtained with the $\varepsilon$-algorithm is positive.

Let us see if the accelerated primal and dual sequences converge to the exact solution $u_{ex}$. For that, let us compute the error between $(u^{primal})^k$ and $u_{ex,h}$, and between $(u^{dual})^k$ and $u_{ex,h}$. We recall that we denote by $u_{ex,h}$ the exact solution $u_{ex}$ discretized on the mesh. We obtained:
<table><caption>Evolution of the error $||u^k - u_{ex,h}||_2$ on the mesh $32 \times 32$</caption><thead><tr><th>k</th><th>Original primal seq.</th><th>Accelerated primal seq.</th><th>Original dual seq.</th><th>Accelerated dual seq.</th></tr></thead><tbody><tr><td>2</td><td>1.722941e-02</td><td>1.711647e-02</td><td>1.214998e-02</td><td>1.197910e-02</td></tr><tr><td>4</td><td>1.778901e-03</td><td>1.500327e-03</td><td>1.067636e-02</td><td>1.066647e-02</td></tr><tr><td>6</td><td>5.633050e-04</td><td>5.967955e-04</td><td>1.076340e-02</td><td>1.074696e-02</td></tr><tr><td>8</td><td>6.968387e-04</td><td>7.505946e-04</td><td>1.078295e-02</td><td>1.077943e-02</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>34</td><td>7.971614e-04</td><td>7.971638e-04</td><td>1.079182e-02</td><td>1.079182e-02</td></tr><tr><td>36</td><td>7.971627e-04</td><td>7.971638e-04</td><td>1.079182e-02</td><td>1.079182e-02</td></tr><tr><td>38</td><td>7.971633e-04</td><td>7.971638e-04</td><td>1.079182e-02</td><td>1.079182e-02</td></tr></tbody></table>

<table><caption>Evolution of the error $||u^k - u_{ex,h}||_2$ on the mesh $128 \times 128$</caption><thead><tr><th>k</th><th>Original primal seq.</th><th>Accelerated primal seq.</th><th>Original dual seq.</th><th>Accelerated dual seq.</th></tr></thead><tbody><tr><td>2</td><td>1.843409e-02</td><td>1.821347e-02</td><td>8.056016e-03</td><td>7.610406e-03</td></tr><tr><td>4</td><td>2.475368e-03</td><td>2.053853e-03</td><td>2.677584e-03</td><td>2.677001e-03</td></tr><tr><td>6</td><td>5.972997e-04</td><td>3.666791e-04</td><td>2.701141e-03</td><td>2.694457e-03</td></tr><tr><td>8</td><td>1.786508e-04</td><td>7.273517e-05</td><td>2.705395e-03</td><td>2.693655e-03</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>58</td><td>4.915487e-05</td><td>4.916213e-05</td><td>2.708191e-03</td><td>2.708191e-03</td></tr><tr><td>60</td><td>4.915672e-05</td><td>4.916213e-05</td><td>2.708191e-03</td><td>2.708191e-03</td></tr><tr><td>62</td><td>4.915809e-05</td><td>4.916213e-05</td><td>2.708191e-03</td><td>2.708191e-03</td></tr></tbody></table>

<table><caption>Evolution of the error $||u^k - u_{ex,h}||_2$ on the mesh $512 \times 512$</caption><thead><tr><th>k</th><th>Original primal seq.</th><th>Accelerated primal seq.</th><th>Original dual seq.</th><th>Accelerated dual seq.</th></tr></thead><tbody><tr><td>2</td><td>1.856981e-02</td><td>1.822785e-02</td><td>7.486792e-03</td><td>6.996415e-03</td></tr><tr><td>4</td><td>2.574905e-03</td><td>2.124638e-03</td><td>6.727851e-04</td><td>7.315175e-04</td></tr><tr><td>6</td><td>6.711839e-04</td><td>4.024026e-04</td><td>6.856359e-04</td><td>6.987853e-04</td></tr><tr><td>8</td><td>2.292313e-04</td><td>8.245509e-05</td><td>6.738612e-04</td><td>6.705496e-04</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>112</td><td>3.094151e-06</td><td>3.094262e-06</td><td>6.703425e-04</td><td>6.703425e-04</td></tr><tr><td>114</td><td>3.094168e-06</td><td>3.094262e-06</td><td>6.703425e-04</td><td>6.703425e-04</td></tr><tr><td>116</td><td>3.094182e-06</td><td>3.094262e-06</td><td>6.703425e-04</td><td>6.703425e-04</td></tr></tbody></table>

When the size of our mesh decreases, the numerical solution converges to the exact solution, for the original and accelerated sequences. In $L^2$-norm, the DDFV scheme is of order 2. This order of convergence is preserved with the accelerated sequence.

### 3.1.2 Problem 2: $\kappa$ discontinuous

We take again the equation (8) but working with a discontinuous diffusion coefficient $\kappa$. This discontinuity only occurs along the faces of the cells in the primal mesh. We have,

$$ \kappa(x, y) = \begin{cases} 1 & \text{if } x \le \frac{1}{2}, \\ 2 & \text{if } x > \frac{1}{2}. \end{cases} $$

The parameters that differ from problem 1 are

$f(x, y) = 2\pi^2 \cos(\pi x)\cos(\pi y) + 20$ and $g(x, y)$ equals to the exact solution on $\Gamma_D$.
The exact solution is

$$u_{ex}(x, y) = \begin{cases} \cos(\pi x) \cos(\pi y) - 10x^2 + 12 & \text{if } x \le \frac{1}{2}, \\ \frac{1}{2} \cos(\pi x) \cos(\pi y) - 5x^2 + \frac{43}{4} & \text{if } x > \frac{1}{2}. \end{cases}$$

After the discretisation of the problem on the primal and dual meshes, we have to resolve a non-linear system like (9). To solve this non-linear system, we use Picard's successive iteration method, starting with the same initial function $u^0$ as the previous problem.

**First case** : The primal mesh is the deformed cartesian mesh seen in the problem 1.

<figure><img src="image_32.png" /><figcaption>Figure 9: $u^0$ on the primal mesh 32 × 32</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>Figure 10: $u_{ex}$ on the primal mesh 32 × 32</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>Figure 11: $u^k$ on the primal mesh 32 × 32</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>Figure 12: Error $|u_{ex} - u^k|$ on the primal mesh 32 × 32</figcaption></figure>

Here are the results obtained by accelerating the convergence of the initial sequence, generated by Picard's method, in using the vectorial $\epsilon$-algorithm, working with quadruple precision on the computation of the accelerated sequence terms and imposing a tolerance of $10^{-12}$ on the error :
<table><caption>Number of iterations to reach the imposed tolerance</caption><thead><tr><td>Mesh</td><td>Original sequence</td><td>Accelerated sequence</td></tr></thead><tbody><tr><td>32 × 32</td><td>94</td><td>50</td></tr><tr><td>64 × 64</td><td>98</td><td>48</td></tr><tr><td>128 × 128</td><td>99</td><td>48</td></tr><tr><td>256 × 256</td><td>100</td><td>52</td></tr><tr><td>512 × 512</td><td>101</td><td>52</td></tr></tbody></table>

<table><caption>Execution time of the while loop</caption><thead><tr><td>Mesh</td><td>Original sequence</td><td>Accelerated sequence</td><td>Gain</td></tr></thead><tbody><tr><td>32 × 32</td><td>6.46 sec</td><td>3.85 sec</td><td>40.40%</td></tr><tr><td>64 × 64</td><td>32.01 sec</td><td>16.51 sec</td><td>48.42%</td></tr><tr><td>128 × 128</td><td>313 sec = 5 min 13 sec</td><td>155 sec = 2 min 35 sec</td><td>50.48%</td></tr><tr><td>256 × 256</td><td>3303 sec = 55 min 03 sec</td><td>1725 sec = 28 min 45 sec</td><td>47.77%</td></tr><tr><td>512 × 512</td><td>38 635 sec ≈ 10 h 43 min</td><td>19 951 sec ≈ 5 h 32 min</td><td>48.36%</td></tr></tbody></table>

The execution times in the table comes from an average after running the code 20 times for meshes 32, 64 and 128; 5 times for the mesh 256 with the accelerated sequence and 1 time with the original one; and 3 times for the mesh 512 with the accelerated sequence and 1 time for the original sequence.

For every meshes, the number of iterations needed to reach tolerance is reduced. Then, it has an impact on execution time. Indeed, the execution time is reduced for all meshes, from 40.40% to 50.48% for the best one.

Let us look at the time taken by the $\varepsilon$-algorithm in the computation of the accelerated sequence.

<table><caption>Execution time of the $\varepsilon$-algorithm in the while loop</caption><thead><tr><td>Mesh</td><td>Accelerated sequence</td><td>$\varepsilon$-algorithm</td><td>Rate</td></tr></thead><tbody><tr><td>32 × 32</td><td>3.85 sec</td><td>0.65 sec</td><td>16.88%</td></tr><tr><td>64 × 64</td><td>16.51 sec</td><td>1.90 sec</td><td>11.51%</td></tr><tr><td>128 × 128</td><td>155 sec = 2 min 35 sec</td><td>15.59 sec</td><td>10.06%</td></tr><tr><td>256 × 256</td><td>1725 sec = 28 min 45 sec</td><td>44.30 sec</td><td>2.55%</td></tr><tr><td>512 × 512</td><td>19 951 sec ≈ 5 h 32 min</td><td>145 sec ≈ 2 min 25 sec</td><td>0.73%</td></tr></tbody></table>

We have the same conclusion than the previous problem. The execution time of the $\varepsilon$-algorithm counts for nothing in all the computation of the accelerated sequence compared to the fixed point method.

We show below the evolution of the relative error for primal and dual sequences, original and accelerated, with the mesh 512 × 512 :
<table><caption>Evolution of the relative error $\frac{\|\mathbf{u}^{k+1} - \mathbf{u}^k\|_2}{\|\mathbf{u}^k\|_2}$ on the deformed mesh 512 × 512</caption><thead><tr><th>k</th><th>Original primal seq.</th><th>Accelerated primal seq.</th><th>Original dual seq.</th><th>Accelerated dual seq.</th></tr></thead><tbody><tr><td>2</td><td>2.941657e-01</td><td>2.965763e+01</td><td>2.947488e-01</td><td>2.967602e+01</td></tr><tr><td>4</td><td>3.916448e-02</td><td>7.445862e-02</td><td>3.919937e-02</td><td>7.547522e-02</td></tr><tr><td>6</td><td>1.191368e-02</td><td>5.647665e-02</td><td>1.191478e-02</td><td>5.645852e-02</td></tr><tr><td>8</td><td>4.582234e-03</td><td>6.531037e-04</td><td>4.581756e-03</td><td>6.486434e-04</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>48</td><td>5.927514e-08</td><td>5.867576e-12</td><td>5.927519e-08</td><td>1.575381e-12</td></tr><tr><td>50</td><td>3.752934e-08</td><td>4.345110e-13</td><td>3.752930e-08</td><td>1.024804e-12</td></tr><tr><td>52</td><td>2.385586e-08</td><td>8.821833e-13</td><td>2.385584e-08</td><td>7.426335e-13</td></tr></tbody></table>

The graph associated with this table is available in [D].

The accelerated primal and dual sequences converge to the solution of the problem. Moreover, during the computation, we check that this numerical solution obtained with the $\varepsilon$-algorithm is positive.

We have also checked that the order 2 in $L^2$-norm of the DDFV scheme is preserved with the accelerated sequence.

**Second case** : The primal mesh is the random mesh seen in the problem 1. However, as we work with a discontinuity which occurs at $x = 1/2$, we should redo the computations with a random mesh where there is a clear demarcation along the straight line $x = 1/2$. So, in this case, we don't have the order 2 in $L^2$-norm.

Here are the results obtained by accelerating the convergence of the initial sequence, generated by Picard's method, in using the vectorial $\varepsilon$-algorithm, working with quadruple precision on the computation of the accelerated sequence terms and imposing a tolerance of $10^{-12}$ on the error :

<table><caption>Number of iterations to reach the imposed tolerance</caption><thead><tr><th>Mesh</th><th>Original sequence</th><th>Accelerated sequence</th></tr></thead><tbody><tr><td>32 × 32</td><td>62</td><td>34</td></tr><tr><td>64 × 64</td><td>117</td><td>54</td></tr><tr><td>128 × 128</td><td>134</td><td>66</td></tr><tr><td>256 × 256</td><td>136</td><td>80</td></tr><tr><td>512 × 512</td><td>212</td><td>154</td></tr></tbody></table>

<table><caption>Execution time of the while loop</caption><thead><tr><th>Mesh</th><th>Original sequence</th><th>Accelerated sequence</th><th>Gain</th></tr></thead><tbody><tr><td>32 × 32</td><td>4.09 sec</td><td>2.71 sec</td><td>33.74%</td></tr><tr><td>64 × 64</td><td>38.75 sec</td><td>19.25 sec</td><td>50.32%</td></tr><tr><td>128 × 128</td><td>424 sec = 7 min 04 sec</td><td>224 sec = 3 min 44 sec</td><td>47.17%</td></tr><tr><td>256 × 256</td><td>4066 sec = 1 h 08 min</td><td>2463 sec = 41 min 03 sec</td><td>39.42%</td></tr><tr><td>512 × 512</td><td>77 588 sec ≈ 21 h 33 min</td><td>58 983 sec ≈ 16 h 23 min</td><td>23.98%</td></tr></tbody></table>

The execution times in the table comes from an average after running the code 20 times for meshes 32 and 64; 10 times for the mesh 128; 2 times with the accelerated sequence and 1 time with the original one for the mesh 256; and 1 time for the mesh 512.

The execution time is reduced for all meshes, from 23.98% to 50.32% for the best one.

Let us look at the time taken by the $\varepsilon$-algorithm in the computation of the accelerated sequence.
<table><caption>Execution time of the ε-algorithm in the while loop</caption><thead><tr><th>Mesh</th><th>Accelerated sequence</th><th>ε-algorithm</th><th>Rate</th></tr></thead><tbody><tr><td>32 × 32</td><td>2.71 sec</td><td>0.31 sec</td><td>11.44%</td></tr><tr><td>64 × 64</td><td>19.25 sec</td><td>2.38 sec</td><td>12.36%</td></tr><tr><td>128 × 128</td><td>244 sec = 3 min 44 sec</td><td>25.71 sec</td><td>10.54%</td></tr><tr><td>256 × 256</td><td>2463 sec = 41 min 03 sec</td><td>95 sec ≈ 1 min 35 sec</td><td>3.86%</td></tr><tr><td>512 × 512</td><td>55 983 sec ≈ 16 h 23 min</td><td>1 181 sec ≈ 19 min 41 sec</td><td>2.11%</td></tr></tbody></table>

We have the same conclusion than the previous case. Despite the use of the quadruple precision, the execution time of the ε-algorithm counts for nothing in all the computation of the accelerated sequence compared to the fixed point method.

We show below the evolution of the relative error for primal and dual sequences, original and accelerated, with the mesh 512 × 512 :

<table><caption>Evolution of the relative error $\frac{\|\mathbf{u}^{k+1} - \mathbf{u}^k\|_2}{\|\mathbf{u}^k\|_2}$ on the deformed mesh 512 × 512</caption><thead><tr><th>k</th><th>Original primal seq.</th><th>Accelerated primal seq.</th><th>Original dual seq.</th><th>Accelerated dual seq.</th></tr></thead><tbody><tr><td>2</td><td>2.941657e-01</td><td>2.965763e+01</td><td>2.947488e-01</td><td>2.967602e+01</td></tr><tr><td>4</td><td>3.916448e-02</td><td>7.445862e-02</td><td>3.919937e-02</td><td>7.547522e-02</td></tr><tr><td>6</td><td>1.191368e-02</td><td>5.647665e-02</td><td>1.191478e-02</td><td>5.645852e-02</td></tr><tr><td>8</td><td>4.582234e-03</td><td>6.531037e-04</td><td>4.581756e-03</td><td>6.486434e-04</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>48</td><td>5.927514e-08</td><td>5.867576e-12</td><td>5.927519e-08</td><td>1.575381e-12</td></tr><tr><td>50</td><td>3.752934e-08</td><td>4.345110e-13</td><td>3.752930e-08</td><td>1.024804e-12</td></tr><tr><td>52</td><td>2.385586e-08</td><td>8.821833e-13</td><td>2.385584e-08</td><td>7.426335e-13</td></tr></tbody></table>

The graph associated with this table is available in [D].

The accelerated primal and dual sequences converge to the solution of the problem. Moreover, during the computation, we check that this numerical solution obtained with the ε-algorithm is positive.

## 3.2 Non-stationary diffusion problem

The problem studied is the following one :

$$ \left\{ \begin{array}{ll} \partial_t u - \nabla \cdot (\kappa \nabla u) + \lambda u = f & \text{in } \Omega \times (0, T], \\ u = g & \text{on } \Gamma_D \times (0, T], \\ \kappa \nabla u \cdot \mathbf{n} = g & \text{on } \Gamma_N \times (0, T], \\ u(\cdot, 0) = u^0 & \text{in } \Omega, \end{array} \right. \qquad (10) $$

where $\Omega$ is a bounded open domain of $\mathbb{R}^2$, with $\partial\Omega = \Gamma_D \cup \Gamma_N (\Gamma_D \cap \Gamma_N = \emptyset)$ and $\mathbf{n} \in \mathbb{R}^2$ the outgoing unit normal vector. The time interval is defined by $[0,T]$, for a given simulation time $T > 0$.

We have $u^0, f, \lambda \in L^2(\Omega)$, with $\lambda \ge 0$, and $g \in H^{1/2}(\partial\Omega)$.

The tensor-valued diffusion coefficient $\kappa$ is supposed to be bounded and to satisfy the following uniform ellipticity condition,
$$ \forall \mathbf{x} \in \Omega, \forall \mathbf{y} \in \mathbb{R}^2, \alpha_{min} \|\mathbf{y}\|^2 \leq \mathbf{y}^t \boldsymbol{\kappa}(\mathbf{x}) \mathbf{y} \leq \alpha_{max} \|\mathbf{y}\|^2, $$

where $\alpha_{min}$ and $\alpha_{max}$ are positive coefficients.

### 3.2.1 Problem 3: $\kappa$ constant

Let us take the domain $\Omega =]0, 1]^2$ over the time interval $[0,1]$.

The diffusion coefficient $\kappa$ is a constant scalar function defined by $\kappa = 1$.

The functions $f$ and $g$ are defined by $f: \mathbb{R}^2 \times \mathbb{R}^+ \to \mathbb{R}: (x, y, t) \mapsto 0$ and $g: \mathbb{R}^2 \times \mathbb{R}^+ \to \mathbb{R}: (x, y, t) \mapsto 0$. The initial function $u^0$ is defined by $u^0: \mathbb{R}^2 \to \mathbb{R}: (x, y) \mapsto \sin(\pi x)\sin(\pi y)$.

We want to solve the preceding problem (10) with $\lambda = 0$ and $\Gamma_N = \emptyset$:

$$ \left\{ \begin{array}{ll} \partial_t u - \nabla \cdot (\kappa \nabla u) = f & \text{in } \Omega \times (0, 1], \\ u = g & \text{on } \Gamma_D \times (0, 1], \\ u(\cdot, 0) = u^0 & \text{in } \Omega. \end{array} \right. \qquad (11) $$

The exact solution of this problem is

$$ u_{ex}: \mathbb{R}^2 \times \mathbb{R}^+ \to \mathbb{R}: (x, y, t) \mapsto e^{-2\pi^2 t} \sin(\pi x) \sin(\pi y). $$

The method used to solve this problem is also the finite volume method DDFV. Indeed, in discretizing in time (11) with the time step $\Delta t = \frac{1}{N}$ (with $N \in \mathbb{N}^*$ fixed), we obtain

$$ \frac{u^{n+1} - u^n}{\Delta t} - \nabla \cdot (\kappa \nabla u^{n+1}) = 0 $$

which is equivalent to

$$ u^{n+1} - \Delta t \nabla \cdot (\kappa \nabla u^{n+1}) = u^n. $$

This last equation is the same than (7) with $\lambda = 1$ and $f = u^n$.

It means that at each time step, we have to solve a non-linear system (9).

To solve this problem, we work on the deformed mesh defined in section 4.1.1. At each time step, we start the Picard's method with the vector $\mathbf{u}_0 = 1.0$.

Here are the results obtained by accelerating the convergence of the initial sequence generated by Picard's method at each time step, in using the vectorial $\varepsilon$-algorithm, in working with quadruple precision on the computation of the accelerated sequence terms and in imposing a tolerance of $10^{-12}$ on the error at each time step:

<table><thead><tr><th colspan="5">Execution time</th></tr><tr><th>Mesh</th><th>Time step</th><th>without acceleration process</th><th>with acceleration process</th><th>Gain</th></tr></thead><tbody><tr><td>32 × 32</td><td>1/128</td><td>10 min 43 sec</td><td>7 min 45 sec</td><td>27.68 %</td></tr><tr><td>64 × 64</td><td>1/256</td><td>103 min 34 sec ≈ 1 h 43 min</td><td>60 min 43 sec ≈ 1 h</td><td>41.37 %</td></tr><tr><td>128 × 128</td><td>1/512</td><td>1556 min 34 sec ≈ 25 h 56 min</td><td>734 min ≈ 12 h 15 min</td><td>52.84 %</td></tr></tbody></table>
The execution times in the previous table comes from an average after running the code 2 times for the mesh 32; once time for the mesh 64; 2 times with the acceleration process and 1 time without the acceleration process for the mesh 128.

The execution time is reduced for all meshes. We gain from 27.68% to 52.84% on the initial computation time.

Let us have a look at the time taken by the ε-algorithm in the computation of the acceleration process.

<table><thead><tr><th colspan="5">Execution time of the ε-algorithm in the accelerated process</th></tr><tr><th>Mesh</th><th>Time step</th><th>Accelerated process</th><th>ε-algorithm</th><th>Rate</th></tr></thead><tbody><tr><td>32 × 32</td><td>1/128</td><td>7 min 45 sec</td><td>65.73 sec ≈ 1 min 5 sec</td><td>14.14%</td></tr><tr><td>64 × 64</td><td>1/256</td><td>60 min 43 sec ≈ 1 h</td><td>366 sec = 6 min 6 sec</td><td>10.04%</td></tr><tr><td>128 × 128</td><td>1/512</td><td>734 min ≈ 12 h 15 min</td><td>4856 sec ≈ 1 h 20 min 56 sec</td><td>11.03%</td></tr></tbody></table>

For every meshes, the ε-algorithm represents around 10% of the total computing time of the acceleration process. It costs almost nothing. Finally, the ε-algorithm seems to be a good way to reduce the computation time of a non-stationary diffusion problem resolved with a Picard's method, in the case of a κ constant.

We find below the evolution in time of a cross-section for the computed function with the acceleration process:

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

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

### 3.2.2 Problem 4: κ discontinuous

Let us take the domain $\Omega = ]0, 1]^2$ over the time interval $[0,3]$.

We are still working on the equation (11) but with a discontinuous diffusion coefficient κ defined by

$$ \kappa(x, y) = \begin{cases} 1/5 & \text{if } x \le \frac{1}{2}, \\ 4/5 & \text{if } x > \frac{1}{2}. \end{cases} $$

This discontinuity only occurs along the faces of the cells from the primal mesh.
<figure><img src="image_38.png" /><figcaption>Figure 14: Evolution in time of the solution computed along the coordinate $y = 0.55$ with the mesh 128</figcaption></figure>

The parameters that differ from problem 3 are

$$f(x, y, t) = e^{-t} \cos(\pi x) \sin(\pi y) \left( \frac{8\pi^2}{5} - 1 \right) \quad \text{in } \Omega$$

and

$g(x, y, t)$ equals to the exact solution on $\Gamma_D$.

The initial condition is

$$u^0(x, y) = \begin{cases} 4 \cos(\pi x) \sin(\pi y) + 15 & \text{if } x \le \frac{1}{2}, \\ \cos(\pi x) \sin(\pi y) + 15 & \text{if } x > \frac{1}{2}. \end{cases}$$

and the exact solution is

$$u_{ex}(x, y, t) = \begin{cases} 4e^{-t} \cos(\pi x) \sin(\pi y) + 15 & \text{if } x \le \frac{1}{2}, \\ e^{-t} \cos(\pi x) \sin(\pi y) + 15 & \text{if } x > \frac{1}{2}. \end{cases}$$

To solve this problem, we proceed in the same way as the problem 3. At each time step, we have to solve a non-linear system (9) and we start the Picard's method with the vector $\mathbf{u}_0 = 1.0$. However, we work on a random mesh as follows:

There is a clear demarcation along the straight line $x = 1/2$ on the primal mesh. We preserved the order 2 in $L^2$-norm of the DDFV scheme.

Here are the results obtained by accelerating the convergence of the initial sequence generated by Picard's method at each time step, in using the vectorial $\varepsilon$-algorithm, in working with quadruple precision on the computation of the accelerated sequence terms and in imposing a tolerance of $10^{-12}$ on the error at each time step:
<figure><img src="image_39.png" /><figcaption>Figure 15: Primal mesh 32 × 32 used</figcaption></figure>

<figure><img src="image_40.png" /><figcaption>Figure 16: Dual mesh 32 × 32 used</figcaption></figure>

<table><thead><tr><th colspan="5">Execution time</th></tr><tr><th>Mesh</th><th>Time step</th><th>without acceleration process</th><th>with acceleration process</th><th>Gain</th></tr></thead><tbody><tr><td>32 × 32</td><td>1/32</td><td>43,53 sec</td><td>37.00 sec</td><td>15.02 %</td></tr><tr><td>64 × 64</td><td>1/64</td><td>354.5 sec ≈ 5 min 54 sec</td><td>322.3 sec ≈ 5 min 22 sec</td><td>9.08 %</td></tr><tr><td>128 × 128</td><td>1/128</td><td>95 min ≈ 1 h 35 min</td><td>84 min ≈ 1 h 24 min</td><td>10.98 %</td></tr><tr><td>256 × 256</td><td>1/256</td><td>1717 min ≈ 28 h 37 min</td><td>1456 min ≈ 24 h 16 min</td><td>15.20 %</td></tr></tbody></table>

The execution times in the previous table comes from an average after running the code 10 times for the mesh 32; 6 time for the mesh 64; 1 time for the mesh 128 and 256. We gain around 10% for the meshes 64 and 128 on the initial computation time, and we can gain around 15% for the meshes 32 and 256.

Let us have a look at the time taken by the ε-algorithm in the computation of the acceleration process.

<table><thead><tr><th colspan="5">Execution time of the ε-algorithm in the accelerated process</th></tr><tr><th>Mesh</th><th>Time step</th><th>Accelerated process</th><th>ε-algorithm</th><th>Rate</th></tr></thead><tbody><tr><td>32 × 32</td><td>1/32</td><td>37.00 sec</td><td>2.78 sec</td><td>7.51%</td></tr><tr><td>64 × 64</td><td>1/64</td><td>322.3 sec ≈ 5 min 22 sec</td><td>15.74 sec</td><td>4.89%</td></tr><tr><td>128 × 128</td><td>1/128</td><td>84 min ≈ 1 h 24 min</td><td>331 sec ≈ 5 min 31 sec</td><td>6.57%</td></tr><tr><td>256 × 256</td><td>1/256</td><td>1456 min ≈ 24 h 16 min</td><td>1853 sec ≈ 30 min 53 sec</td><td>2.12%</td></tr></tbody></table>

For every meshes, the ε-algorithm represents less than 8% of the total computing time of the acceleration process. It costs only 2.12% for the mesh 256. Finally, the ε-algorithm seems to be a good way to reduce the computation time of a non-stationary diffusion problem resolved with a Picard's method, in the case of a discontinuous κ.
<figure><img src="image_41.png" /><figcaption>Figure 17: Evolution in time of the solution computed along the coordinate $y = 0.51$ with the mesh $32 \times 32$</figcaption></figure>

# 4 Conclusion

The aim of this study was to accelerate the convergence of numerical sequences. To do so, there are several ways to accelerate the convergence of a sequence. At first, we decided to study the $\varepsilon$-algorithm because while being efficient, its implementation is relatively easy.

To begin with, we discovered how this algorithm works for scalar and vectors sequences. Then, some theoretical results were presented and the implementation of the algorithm was explained. We have seen that it is interesting to work with quadruple precision when we compute the terms of the accelerated sequence. It allows us to avoid the truncation errors as far as possible, with very little impact on the computation time of the $\varepsilon$-algorithm compared to the overall computation time. With some test cases, we were able to identify some practical results for this algorithm.

Indeed, the $\varepsilon$-algorithm allows us to accelerate linearly converging sequences, but, accelerating sequences which converge in a quadratic way seems to be extremely complicated. Applied to stationary and non-stationary diffusion problems where the equations are solved by using a fixed point method, the $\varepsilon$-algorithm seems to be a good way to accelerate the convergence of these sequences.

Based on the results obtained, it may be interesting to often use a convergence acceleration method like the $\varepsilon$-algorithm when the iterative fixed point method is used. Indeed, we reduce the computation time for a low cost use of the $\varepsilon$-algorithm.

However, we have seen experimentally that it is difficult to accelerate the convergence of sequence which converges in a quadratic way. Moreover, we haven't mentioned it, but we have also difficulties in accelerating the convergence of sequences which converge in a logarithmic way. So, it can be interesting to study other accelerated processes which could be able to accelerate these types of sequences. For that, we can mention the $\rho$-algorithm and the $\theta$-algorithm, both available in [Bre78].
# Appendices

## Appendix A Laptop information

**Architecture**: x86_64

Processor operating mode(s): 32-bit; 64-bit

Address sizes: 39 bits physical, 48 bits virtual

Byte order: Little Endian

**Processors**: 8

**Manufacturer identifier**: Genuine Intel

Model name: Intel (R) Core (TM) i5-8250U CPU @ 1.60GHz

Processor family: 6

Model: 142

Thread(s) per cores: 2

Core(s) per socket: 4

Socket(s): 1

Max processor speed in MHz: 3 400,00

Min processor speed in MHz: 400,00

Bogo MIPS: 3600,00

**Virtualization features**:

Virtualisation: VT-x

**Caches** (sum of all):

L1d: 128 KiB ( 4 instances )

L1i: 128 KiB ( 4 instances )

L2: 1 MiB ( 4 instances )

L3: 6 MiB ( 1 instance )
Appendix B Particular rule

Let us denote $W := \varepsilon_{k-3}^{(n+2)}$, $P := \varepsilon_{k-2}^{(n+1)}$, $R := \varepsilon_{k-2}^{(n+2)}$, $N := \varepsilon_{k-1}^{(n)}$, $C := \varepsilon_{k-1}^{(n+1)}$, $S := \varepsilon_{k-1}^{(n+2)}$, $T := \varepsilon_k^{(n)}$, $U := \varepsilon_k^{(n+1)}$ and $E := \varepsilon_{k+1}^{(n)}$.

$N$

The situation is the following one :

$P$

$T$

$W$

$C$

$E$

$R$

$U$

$S$

The particular rule is rewriting

$$E = \frac{a}{\left(1 + \frac{a}{C}\right)}$$

with

$$a = \frac{S}{\left(1 - \frac{S}{C}\right)} + \frac{N}{\left(1 - \frac{N}{C}\right)} - \frac{W}{\left(1 - \frac{W}{C}\right)}.$$

By the recurrence (1), we have :

$$C = W + \frac{1}{R - P};$$

$$T = P + \frac{1}{C - N};$$

$$U = R + \frac{1}{S - C};$$

$$E = C + \frac{1}{U - T}.$$

By replacing the values of $U$ and $T$ in the last equality, we obtain the rule of the cross

$$E = C + \frac{1}{R + \frac{1}{S-C} - P - \frac{1}{C-N}} = C + \frac{1}{\frac{1}{N-C} + \frac{1}{S-C} - \frac{1}{W-C}}.$$

From this rule of the cross, we get

$$E = C + \frac{1}{\frac{1}{N-C} + \frac{1}{S-C} - \frac{1}{W-C}} = C + \frac{C}{\frac{C}{N-C} + \frac{C}{S-C} - \frac{C}{W-C}}$$

iff

$$E = C - \frac{C}{\frac{C}{C-N} + \frac{C}{C-S} - \frac{C}{C-W}} = C - \frac{C}{\frac{C-N+N}{C-N} + \frac{C-S+S}{C-S} - \frac{C-W+W}{C-W}}$$
⇔

$$E = C - \frac{C}{1 + \frac{N}{C-N} + \frac{S}{C-S} - \frac{W}{C-W}} = C - \frac{C}{1 + \frac{1}{C} \left( \frac{N}{C-N} + \frac{S}{C-S} - \frac{W}{C-W} \right)}$$

⇔

$$E = C - \frac{C}{1 + \frac{a}{C}} = \frac{C\left(1 + \frac{a}{C}\right) - C}{\left(1 + \frac{a}{C}\right)} = \frac{a}{1 + \frac{a}{C}}.$$
Appendix C Graphs for the problem 1 on the mesh 512 × 512

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

<figure><img src="image_43.png" /></figure>
Appendix D Graphs for the problem 2 on the mesh 512 × 512

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

<figure><img src="image_45.png" /></figure>
References

[Bel] BELHANNACHE Farida. *Introduction aux Équations aux Différences*, Maître de conférences, Département de Mathématiques, Université Mohammed Seddik Ben Yahia-Jijel.

[DDFV] BLANC X. , HERMELINE F. , LABOURASSE E. , PATELA J. . *Monotonic dia- mond and DDFV type finite-volume schemes for 2D elliptic problems*, To be published (CICP).

[Bre77] BREZINSKI Claude. *Accélération de la convergence en analyse numérique*, Lecture Notes in Mathematics, Éditions Springer, 1977.

[Bre78] BREZINSKI Claude. *Algorithmes d'accélération de la convergence étude numérique*, Collection langages et algorithmes de l'informatique, Éditions Technip, 1978.

[Bur] BURDEN Richard L., FAIRES J. Douglas. *Numerical Analysis*, Ninth Edition.

[Cuv] CUVELIER François. *Analyse numérique élémentaire*, Notes de cours, Sup GAlilée, Ingénieurs MACS 1ere année, Université Paris XIII/Institut Galilée.

[Hig] HIGHAM Nicholas J. . *Accuracy and Stability of Numerical Algorithms*, Seconde Edition, SIAM, 2002.

[Ple] PLESSIER Alexiane. *Implicit semi-Lagrangian schemes for compressible gas dynamics*. Numerical Analysis [math.NA]. Sorbonne Université, 2023. English. ffNNT : 2023SORUS110ff. fftel-04135003

[Vig] VIGNES Jean. *Algorithmes numériques analyse et mise en œuvre 2 : équations et systèmes non linéaires*, Collection langages et algorithmes de l'informatique, Éditions Technip, 1980.

[Quad] https://en.wikipedia.org/wiki/Quadruple-precision_floating-point_format.