Phuc Ngo
email: hoai-diem-phuc.ngo@loria.fr

Digital line segment detection for table reconstruction in document images

Keywords: Digital line / segment detection, Blurred segment, Adaptive directional scan, Materialized table extraction, Digital geometry

published or not. The documents may come

Digital line segment detection for table reconstruction in document images

Introduction

The process of table analysis involves table detection and table structure recognition. The task of table detection aims to locate regions of image that contain tables using the bounding boxes, while table structure recognition involves the identification of row and column layout information for individual table cells. It is an important task in many document analysis applications. Thank to their compact form, tables are frequently used to summarized information in a structured and relational way. They are present in a large variety of documents such as reports, scientific papers, business documents, invoice, . . . This paper addresses the table detection in document images. It is a difficult problem due to the variety of layouts and formats of tables in images (see Fig. 1). Furthermore, when tables are correctly extracted from images, it allows to perform the recognition task more efficiently, and thus improve the table information extraction such as the text recognition within each cell and the interpretation of tabular data in the document.

In the context of table detection, several approaches have been proposed in the literature. One of the earliest works on identifying tabular regions in document images is of Watanabe et al. [START_REF]Structure analysis of table-form documents on the basis of the recognition of vertical and horizontal line segments[END_REF] in which the authors proposed to identify blocks enclosed by vertical and horizontal line segments. More precisely, line segments are first detected and corner points are subsequently determined. The method then exploits the connective relationships between the extracted corner points for the individual blocks using global and local tree structures. Following this idea, the methods in [START_REF] Cesarini | Trainable table location in document images[END_REF][START_REF] Gatos | Automatic table detection in document images[END_REF][START_REF] Green | Model-based analysis of printed tables[END_REF][START_REF] Héloïse | Straight line reconstruction for fully materialized table extraction in degraded document images[END_REF] use parallel lines, the horizontal and vertical space as features to extract the table regions. Other conventional techniques for table analysis have been studied based on the textual [START_REF] Göbel | Icdar 2013 table competition[END_REF] (a-b) and from UNLV Table Dataset [START_REF] Shahab | An open approach towards the benchmarking of table structure recognition systems[END_REF] (c-d).

layout analysis of documents. We can mention the work of T. Kieninger [START_REF] Kieninger | Table structure recognition based on robust block segmentation[END_REF] which relies on the detection of text block units for extracting the table cells. More precisely, the method starts with an arbitrary word as block seed, and recursively extends this block to all the words that interlace with their vertical neighbors. Then, it uses the gap information to prevent the table columns. Similarly, S. Mandal et al. in [START_REF] Mandal | Simple and effective table detection system from document images[END_REF] supposed that tables are present in distinct columns and the spaces between the fields are larger than those between words in text lines. Then, the table extraction relies on the formation of word blob text lines and on finding the set of text lines that may form a table. In [START_REF] Shafait | Table detection in heterogeneous documents[END_REF], the method is based on the layout analysis module of Tesseract OCR (Optical Character Recognition) to detect tables. More specifically, it uses the tab-stops to inform text aligns and to extract blocks composing the layout of a document. These blocks are then evaluated as table candidates using a specific strategy. Recently, the deep-learning approaches based on convolutional neural network (CNN) have been significantly investigated for table analysis. T. Kieninger and A. Dengel [START_REF] Thomas | The t-recs table recognition and analysis system[END_REF] can be considered as the pioneers to apply unsupervised learning method to the table detection using a bottomup clustering of given word segments. After that, we have the supervised methods such as in [START_REF] Arif | Table detection in document images using foreground and background features[END_REF][START_REF] Gilani | Table detection using deep learning[END_REF] that used Faster R-CNN based model for table detection, or [START_REF] Minghao | Tablebank: A benchmark dataset for table detection and recognition[END_REF][START_REF] Prasad | Cascadetabnet: An approach for end to end table detection and structure recognition from image-based documents[END_REF] used CNN to performed table detection and table structure recognition at the same time.

The present work aims at designing a method of table detection using image processing techniques. In particular, we are interested in materialized tables in documents. Such tables have rectangular form and are composed by straight line segments intersecting at right angle. They may not be closed but must contain at least two surrounded horizontal lines and inside vertical lines (see Fig. 1). For detecting these tables, the method relies on a recent line detector called Fast Blurred Segment Detector (FBSD) which is developed in the context of digital geometry. This detector is fast, accurate and robust to noise. It can be run in pure automatic mode, with quite few parameter settings. The extracted lines are then used as candidates to reconstruct the tables in document images. The method has been evaluated and compared to other state-of-the-art methods on different public datasets and shown a very competitive result. The paper is organized as follows. Sec. 2 recalls the main theoretical notions used in this work. Sec. 3 describes the FBSD detector for extracting straight line segments in images. The proposed method of table extraction is presented in Sec. 4. Then, Sec. 5 shows the experimental results, conclusion and future works are provided in Sec. 6.

Background Notions

We recall hereafter several notions of digital geometry [START_REF] Klette | Digital geometry -Geometric methods for digital picture analysis[END_REF] used in this work. We refer the reader to the given references for more details.

0 ≤ ax + by + c < ω (1)
Hereafter, we note V(D) = (-b, a) the director vector of D, w(D) = ω its arithmetical width, h(D) = -c its shift to origin, and p(D) = max(|a|, |b|) its period. The thickness µ = ω-1 p(D) of D is the minimum of the vertical and horizontal distances between lines ax + by + c = 0 and ax + by + c = ω -1. When ω = max(|a|, |b|) then D is the narrowest 8-connected line and called naif digital line, and ω = |a| + |b| then D is 4-connected line and called standard digital line.

A digital straight segment is a finite subset of D. Example of blurred segment is given in Fig. 2 (c). In [START_REF] Isabelle | Optimal blurred segments decomposition of noisy shapes in linear time[END_REF], it is shown that blurred segments can be detected in linear time by a recognition algorithm based on an incremental growth of the convex hull of added points.

Adaptive Directional Scan

From the notion of blurred segment, two notions of directional scan are proposed in [START_REF] Even | Thick line segment detection with fast directional tracking[END_REF][START_REF] Kerautret | Blurred segments in gray level images for interactive line extraction[END_REF] for thick line detection in gray-level images.

Definition 3. [17]

A static directional scan is an ordered partition into scans S i restricted to the grid domain G ⊂ Z 2 of a thick digital straight line D, called scan strip.

Each scan S i is a segment of a naive line N i , called scan line, orthogonal to D. The directional scan is defined as:

DS = S i = D ∩ N i ∩ G V(N i ) • V(D) = 0 h(N i ) = h(N i-1 ) + p(D) (2)
In this definition, V(N i ) • V(D) = 0 expresses the orthogonality between the scan lines N i and the scan strip D. The shift p(D) between successive scans N i-1 and N i guarantees that all points of D are traversed only one time. The scans S i can be iteratively parsed from the start scan S 0 to both ends (see Fig. 3 (a)).

An adaptive directional scan is a dynamical version of the directional scan with an on-line registration to a moving search direction. Compared to static directional scans where the scan strip remains fixed to the initial line D 0 , here the scan strip D i follows a curve C to track while scan lines N i remain fixed (see Fig. 3 (b)). Definition 4. [START_REF] Even | Thick line segment detection with fast directional tracking[END_REF] An adaptive directional scan is defined by:

ADS =            S i = D i ∩ N i ∩ G V(N i ) • V(D 0 ) = 0 h(N i ) = h(N i-1 ) + p(D 0 ) D i = D( C i , w(D 0 )), i > 0            (3)
where C i is a triplet composed of the director vector (a i , b i ) and the shift to origin c i of an estimate at position i of the tangent to the curve C to track.

The obtained thick digital line is used to update the scan strips and lines. The last clause expresses the scan bounds update at iteration i.

FBSD: Fast Blurred Segment Detector

In this section, we describe briefly the efficient tool FBSD [START_REF] Even | Thick line segment detection with fast directional tracking[END_REF] to detect blurred segments from gray-level images based on discrete geometry notions presented previously. More precisely, the method examines the gradient image of the input for edge detection. Firstly, it uses gradient magnitude to select image points for the edge candidates. Then, it relies on gradient orientation to track edges and to detect the blurred segments.

The obtained segments encode straight edges, and its thickness parameter characterizes the edge sharpness. In particular, a linear-time blurred segment recognition algorithm is applied to incrementally extend the segment and the adaptive directional scan allows to orientate the edge candidates in appropriate directions in gradient image.

The process for a single detection, given an input seed segment AB crossing the edges to be detected, is summarized in Fig. 4. It is composed of two steps. The initial detection consists in building and extending a rough blurred segment B of assigned thickness ε 0 , based on points with highest gradient magnitude found in each scan of the directional scan defined by AB with a large tolerance to gradient direction. Validity tests, for small and sparse segments, are then applied to decide the detection pursuit.

In the second step of fine tracking, the final blurred segment B is built and extended with points that correspond to local maxima of the gradient, ranked by magnitude order, and with gradient direction close to start point gradient direction. At this refinement step, the adaptive directional scan, defined by the found center position C and the rough blurred segment direction D, is used in order to extend the segment in the appropriate direction. As soon as the thickness of the expanding blurred segment gets stable, a procedure is applied to focalize the search direction and to avoid further insertion of ouliers. Final length and sparsity thresholds of the output segment B can be set accordingly. They are the only parameters of this local detector, together with the input assigned thickness ε 0 .

By considering the prominent local maxima of the gradient magnitude under AB, the previous process also allows the detection of all the edges crossed by the seed segment AB, called multi-detection. Furthermore, in order to avoid multiple detections of the same edge, an occupancy mask is used to collect the dilated points of all the detected blurred segments, so that these points can not be used any more in the detection.

In [START_REF] Even | Thick line segment detection with fast directional tracking[END_REF], a full automatic detection of all straight lines in image is also proposed. More precisely, a line segment, that crosses the whole image, is swept in both directions, vertical then horizontal, from the center to the borders of the image. At each position, the multi-detection algorithm is run to collect all the blurred segments found under the sweeping segment. Then, small segments are rejected in order to avoid the formation of misaligned segments when the sweeping segment crosses an image edge near one of its ends. For more details of the FBSD detector, we refer the reader to the paper [START_REF] Even | Thick line segment detection with fast directional tracking[END_REF].

Removing segments belonging to text

Table cell extraction

Input (gray-level) image

Main Table Extraction Framework

This section presents the proposed method for extracting tables from document images.

In [START_REF] Héloïse | Straight line reconstruction for fully materialized table extraction in degraded document images[END_REF], a table is said fully materialized if it is composed and surrounded by straight line segments intersecting at right angle. In our work, we concern with tables of rectangular form which are also composed by straight line segments intersecting at right angle. However, they may not be closed, compared to fully materialized tables, but must contain at least two surrounded horizontal lines and inside vertical lines. We call them materialized tables. Such tables are quite frequent in a wide variety of documents such as invoices, scientific papers, commercial / administrative documents, question forms, . . . . Some examples of materialized tables are given in Fig. 1.

The proposed method relies on the edge line detector FBSD described in the previous section. The lines in document images are not all table edges, they can be part of graphics, images or separators. It should be mentioned that the acquisition noise during the sequential scanning / printing process may disrupt, damage or misalign the straight line segments (see Fig. 1 (c-d)) and it makes the task of table reconstruction and detection more difficult. This leads us to a proposal of line segment recovery before the reconstruction of tables from document images.

The automatic table extraction process is composed of 6 steps. These steps are detailed in the following. The overview of the whole process is given in Fig. 5 and illustrated in Fig. 6.

Line segment detection: This first step consists in extracting line segments in document images (see Fig. 6 (b)). For this, we use the FSBS detector [START_REF] Even | Thick line segment detection with fast directional tracking[END_REF]. These lines consist candidate segments of tables to be extracted and also other segments belonging to graphics, images, separators, texts, . . .

Horizontal and vertical segment extraction:

As we address only materialized tables which are composed of horizontal and vertical segments intersecting at right angle, we thus remove all segments which are not horizontal and vertical (see Fig. 6 (c)). A tolerance of α = ± 5 • angle is applied for such segment selection due to impact of scanning procedure of documents.

Line segment recovery: In order to deal with broken lines, especially in degraded conditions and due to the acquisition noise, and to ensure a better stability of table detection, the extracted horizontal (resp. vertical) segments are prolonged and connected if they are close to each other (see Fig. 6 (d)). A distance threshold δ = 30 pixels is used to test the closeness between two horizontal (resp. vertical) segments.

Suppression of segments belonging to text:

The extracted segments from the previous step do not always belong to borders of cell tables. In particular, most of them are associated to text in document image. In order to eliminate such segments, we consider the original gray-level image and verify the line model of border cell tables. A local analysis around the extracted segment is performed. It can be observed that the pattern of the 1D gray intensity profile the segment of a cell border is typical with a bottom (of very low intensity) in the middle and high intensity on both sides1 , and it is not always the case for segments associated to text (see Fig. 7). This profile pattern is repetitive along the segment associated to the border cell tables, but not in case of text. Therefore, by simply computing the ratio of the number of profile patterns corresponding to the border cell segment over the segment length, we can distinguish the segments associated to text. In our experiments, the ratio is set at 75% for eliminating text segments, and w = 7 pixels around the extracted segment is used for the local analysis of intensity profile. The result of this step is illustrated in Fig. 6 (e).

Table cell extraction:

This step aims at detecting table cells before the final step of table reconstruction. The cell extraction is based on the vertical and horizontal segments limiting the table cells. Indeed, for the materialized tables, these cells are limited by four segments defining a rectangle. In case of open tables, i.e., they are not bounded by vertical segments (see Fig. 1 (c-d)), the cells at extremities are surrounded by only three segments (two horizontal and one vertical segments), then we can close the cell by using the end points of the horizontal segments. By that way, we can eliminate, from the extracted segments, the separator lines which do not form a close rectangle.

Table reconstruction: Finally, a table is built as unions of all its cells (see Fig. 6 (f)). For this step, we compute the connected components from the extracted cells, and the corresponding table is obtained as the bounding box of each connected component. A verification of table size is performed at the end to remove small bounding boxes related to noise or other elements which may not be the tables. For this, we keep only the bounding boxes of size bigger than 1% of image size.

Experimental Results

This section presents results of the proposed method for table extraction from document images. In the experimental stage, the proposed approach is validated through comparisons with other recent methods in the literature: S. Faisal and S. Ray [START_REF] Shafait | Table detection in heterogeneous documents[END_REF], H. Alhéritière et al. [START_REF] Héloïse | Straight line reconstruction for fully materialized table extraction in degraded document images[END_REF], D. Prasad et al. [START_REF] Prasad | Cascadetabnet: An approach for end to end table detection and structure recognition from image-based documents[END_REF] and M. Li et al. [START_REF] Minghao | Tablebank: A benchmark dataset for table detection and recognition[END_REF]. In [START_REF] Shafait | Table detection in heterogeneous documents[END_REF], the authors proposed a method based on components of the layout analysis module of Tesseract to locate tables in document images. In [START_REF] Héloïse | Straight line reconstruction for fully materialized table extraction in degraded document images[END_REF], the method relies on the local Radon Transform for line extraction and then the pattern recognition techniques for reconstructing tables from the extracted lines. [START_REF] Prasad | Cascadetabnet: An approach for end to end table detection and structure recognition from image-based documents[END_REF] and [START_REF] Minghao | Tablebank: A benchmark dataset for table detection and recognition[END_REF] are deep-learning approaches based on convolution networks which allow to detect the regions of tables and recognize the structural body cells from the detected tables at the same time. For more details about these methods, we refer the reader to the given references.

Comparison results are run on three publicly available datasets for table detection: ICDAR 2013 Table Competition [START_REF] Göbel | Icdar 2013 table competition[END_REF], UNLV Table Dataset [START_REF] Shahab | An open approach towards the benchmarking of table structure recognition systems[END_REF] and Marmot Dataset [START_REF] Fang | Dataset, ground-truth and performance metrics for table detection evaluation[END_REF]. More precisely, the ICDAR 2013 [START_REF]Icdar 2013 table competition dataset[END_REF] contains 128 documents with a total of 150 tables: 75 tables in 27 excerpts from the EU and 75 tables in 40 excerpts from the US Government. The UNLV Table Dataset [START_REF]Unlv table dataset[END_REF] contains 427 examples in scanned image format from a variety of sources ranging from technical reports and business letters to newspapers and magazines. The Marmot Dataset [START_REF]Marmot dataset[END_REF] has 2000 documents mainly from research papers composed of both Chinese and English pages, in which 1000 pages containing at least one table, while the other 1000 pages do not contain tables, but have complex layout where some page components may be mistakenly recognized as tables, such as matrices and figures. All these datasets are provided with information, in XML format, of table structure ground-truth data such as bounding box coordinates, rows, columns, cells, . . .

To evaluate the performance of the table detection algorithms, we use the standard evaluation metrics in the literature [START_REF] Hu | Medium-independent table detection[END_REF][START_REF] Kieninger | An approach towards benchmarking of table structure recognition results[END_REF]: precision (P), recall (R) and F 1 -measure (F 1 ) which is an harmonic mean of P and R. The metrics are computed by summing up the overlapping area of the obtained result and the ground-truth. More precisely, let D be the set of pixels of the extracted tables, and G of the ground-truth. Then, the considered measure are computed as:

P = D ∩ G D , R = D ∩ G G , F 1 = 2 * R * P R + P (4)
The results are reported in Tab. 1 on the overall average score of each dataset. We observe that our method archives well-scored and competitive results comparing to the others on the considered datasets. More precisely, for the ICDAR 2013 dataset, [START_REF] Prasad | Cascadetabnet: An approach for end to end table detection and structure recognition from image-based documents[END_REF] 2 Results extracted from [START_REF] Héloïse | Straight line reconstruction for fully materialized table extraction in degraded document images[END_REF]. 3 Results on a subset of ICDAR 2013 containing 88 tables in [START_REF] Héloïse | Straight line reconstruction for fully materialized table extraction in degraded document images[END_REF]. Fig. 8: Limit cases of our proposal. Images from UNLV [START_REF] Shahab | An open approach towards the benchmarking of table structure recognition systems[END_REF] and Marmot [START_REF] Fang | Dataset, ground-truth and performance metrics for table detection evaluation[END_REF] datasets. performed best at the absolute score of 100% on the evaluation measures, and our proposal is as good as [START_REF] Minghao | Tablebank: A benchmark dataset for table detection and recognition[END_REF]. In case of the UNLV Table dataset, [START_REF] Minghao | Tablebank: A benchmark dataset for table detection and recognition[END_REF] provides the best score and our method comes in second place. This can be explained by the very low quality of scanned document images, and also this dataset contains no materialized tables (see Fig. 8). Finally, on the Marmot dataset, the proposed method has the best precision and F 1 measures, but less good than [START_REF] Minghao | Tablebank: A benchmark dataset for table detection and recognition[END_REF] on recall. This lower recall can be explained by the non detection of materialized tables and the false detection of graphics in this dataset. Though, it should be mentioned that this dataset contains high quality images.

Further comparisons with other methods in the literature are also available in the referenced papers [START_REF] Fang | Dataset, ground-truth and performance metrics for table detection evaluation[END_REF][START_REF] Minghao | Tablebank: A benchmark dataset for table detection and recognition[END_REF][START_REF] Prasad | Cascadetabnet: An approach for end to end table detection and structure recognition from image-based documents[END_REF][START_REF] Shafait | Table detection in heterogeneous documents[END_REF].

Conclusion

In this paper, we presented a method for table extraction from document images based on notions of digital geometry. The method relies on the FBSD detector [START_REF] Even | Thick line segment detection with fast directional tracking[END_REF] to detect the straight line segments in image and then reconstruct the tables from the extracted lines. This detector is fast, accurate and robust to noise. In particular, it allows to handle noisy document images. Indeed, the sequential printing and scanning of a document may deteriorate its contents, including the table structures such as disrupting table lines. The proposed approach has been evaluated and compared with other methods on different table datasets. The results are encouraging and competitive to the existing methods. The source code of the proposed method, based on OpenCV library [START_REF]OpenCV: Open computer vision[END_REF], is available at the GitHub repository: https://github.com/ngophuc/TableExtraction.

In the field of document analysis and interpretation of tabular data in document images, the table structure recognition [START_REF] Farrukh | Interpreting data from scanned tables[END_REF][START_REF] Thomas | The t-recs table recognition and analysis system[END_REF] is also an active research topic. It involves identifying the rows, columns, individual table cells. In an advantage context, it is to retrieve the textual content of the cells. It is widely used in image captioning, video description, and many other applications. In this context, the proposed method can be helpful as a first step of table detection before extracting other information. On the other hand, we can use the method to address the need for generating data for the machine learning approaches which need many data for the training process.

As mentioned previously, we address the materialized tables in this work. However, in the considered datasets, there are tables not for this case study, and also graphics or matrices that our algorithm mistakenly detects as tables, see Fig. 8 for some examples. In future, we would like to extend the method for more general cases of tables and to conduct more comparisons with the state-of-the-art methods on more data.

Fig. 1 :

 1 Fig. 1: Examples of materialized tables in document images. Images extracted from ICDAR 2013 Table Competition [13] (a-b) and from UNLV Table Dataset [26] (c-d).

Fig. 2 :

 2 Fig. 2: Examples of (a) a standard segment belongs to the 4-connected digital line D(2, -3, 3, 5), (b) a naif segment belongs to the 8-connected digital line D(2, -3, 0, 3) and (c) a blurred segment of thickness ε = 1.5 belongs to the digital line D(3, -4, 3, 7).

2. 1

 1 Digital Straight Line and Blurred Segment Definition 1. [24] A digital straight line of integer characteristics (a, b, c, ω) ∈ Z 4 with a and b relatively prime, denoted by D(a, b, c, ω), is the set of digital points p = (x, y) ∈ Z 2 satisfying the inequalities:

Figure 2 (

 2 a-b) show examples of digital straight lines and segments. Definition 2. [16] A blurred segment B of assigned thickness ε is a sequence of points that are all covered by a digital straight line D of thickness µ ≤ ε. The covering digital line with minimal thickness is called the optimal line of B.

Fig. 3 :

 3 Fig. 3: (a) A static directional scan with the start scan S 0 in red, odd scans in green, even scans in blue, bounds of scan lines N i with dashed lines and bounds of scan strip D with bold dashed lines. (b) An adaptive directional scan in which the scan strip is dynamically fit to the curve C, and the scans are accordingly shifted to cover C.

Fig. 4 :

 4 Fig.4: Flowchart of the thick line segment detection process by FBSD[START_REF] Even | Thick line segment detection with fast directional tracking[END_REF].

Fig. 6 :

 6 Fig. 6: Different steps of the table extraction process: (a) input image, (b) straight segments detected by FBSD (green), (c) candidate horizontal (blue) and vertical (red) segments, (d) recovering nearby segments by extending and connecting horizontal (blue) and vertical (red) segments, (e) removing segments associated to text and (f) final result of table extraction (black box).

Fig. 7 :

 7 Fig. 7: Gray intensity profile (red segment) analysis along extracted segment (blue segment): (a) border cell table, and (b) text.

  Flowchart of the proposed method for table extraction from document images.

	Gradient image	FBSD detector	Straight segments	Filtering segments	Horizontal/vertical segments	Segment recovery	Long segments
					Tables	Table reconstruction
	Fig. 5:

Table 1 :

 1 Evaluation results of table detection in document images on different datasets.

	Dataset	Method	P	R	F 1
		Faisal et al. [25] 2	0.25	0.33	0.28
	ICDAR 2013	Alhéritière et al. [14] 3	0.98	0.83	0.90
	Dataset	Prasad et al. [23]	1	1	1
		M. Li et al. [22]	0.97	0.80	0.88
		Our proposal	0.91	0.85	0.88
		Faisal et al. [25]	0.86	0.79	0.82
	UNLV Table	Alhéritière et al. [14]	0.73	0.66	0.70
	Dataset	M. Li et al. [22]	0.92	0.96	0.94
		Our proposal	0.83	0.82	0.83
	Marmot	M. Li et al. [22]	0.77	0.98	0.86
	Dataset	Our proposal	0.93	0.81	0.87

A similar line model has been used in[START_REF] Héloïse | Straight line reconstruction for fully materialized table extraction in degraded document images[END_REF] for line verification and table extraction. Their model consists in separating the 1D intensity profile, from left to right, into 3 zones: the intensity should begin to increase then (potentially) stabilize and decrease.


Ahmad Deeb

Aziz Hamdouni

Dina Razafindralandy

Performance of Borel-Padé-Laplace integrator for the solution of stiff and non-stiff problems

Keywords: Borel-Laplace summation, divergent series, time integrator, stiff equations

A stability analysis of the Borel-Padé-Laplace series summation technique, used as explicit time integrator, is carried out. Its numerical performance on stiff and non-stiff problems is analyzed. Applications to ordinary and partial differential equations are presented. The results are compared with those of many popular schemes designed for stiff and non-stiff equations.

Introduction

Stiff problems occur in many areas of engineering science, such as mechanics, electrical and chemical engineering (see for instance [START_REF] Hairer | Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems, Second Revised[END_REF][START_REF] Iserles | A first course in the numerical analysis of differential equations[END_REF][START_REF] Butcher | Numerical Methods for Ordinary Differential Equations[END_REF][START_REF] Lambert | Numerical Methods for Ordinary Differential Systems: The Initial Value Problem[END_REF][START_REF] Willoughby | Stiff differential systems[END_REF]). However, their resolution has remained a challenge for numerical analysts. The reason is that many numerical methods designed for general ordinary differential equations exhibit a high instability when solving stiff problems, unless an excessively small time step is used. As a consequence, numerical schemes with better stability properties have been developed especially for stiff problems.

One method used to estimate the largest time step allowed by a given numerical scheme without breaking its stability is the analysis of the linear stability domain. The scheme is called A-stable if this domain contains the half complex plane with negative real part, meaning that the scheme is stable in some sense however large is the time step, for the numerical solution of a 1D linear equation. The notion of a linear stability domain will be recalled later. See also [START_REF] Hairer | Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems, Second Revised[END_REF][START_REF] Lambert | Numerical Methods for Ordinary Differential Systems: The Initial Value Problem[END_REF][START_REF] Hackbusch | The Concept of Stability in Numerical Mathematics[END_REF] for different notions of stability. Of course, even if a scheme is A-stable, the time step is limited in pratice due to precision requirements.

Among the most widely used numerical schemes for stiff equations, one can mention implicit linear multistep methods based on backward difference formu-when authors showed that the heat equation and many equations in mechanics (Burgers and Navier-Stokes equations, . . . ), quantum physics or astronomy have divergent but Gevrey Taylor series [START_REF] Lutz | On the Borel summability of divergent solutions of the heat equation[END_REF][START_REF] Lysik | Borel summable solutions of the Burgers equation[END_REF][START_REF] Costin | Borel summability of Navier-Stokes equation in R 3 and small time existence[END_REF][START_REF] Dyson | Divergence of perturbation theory in quantum electrodynamics[END_REF][START_REF] Suslov | Divergent perturbation series[END_REF][START_REF] Kontopoulos | Order and chaos in dynamical astronomy[END_REF]. The Borel-Laplace summation method has been transformed into numerical algorithm [START_REF] Thomann | Formal and numerical summation of formal power series solutions of ODE's[END_REF] and used for the first time as a time integrator by Razafindralandy and Hamdouni [START_REF] Razafindralandy | Time integration algorithm based on divergent series resummation, for ordinary and partial differential equations[END_REF]. Since then, a number of features of the Borel-Laplace integrator was studied. For example, it generally allows much larger time steps than other explicit methods for the resolution of many problems [START_REF] Razafindralandy | Time integration algorithm based on divergent series resummation, for ordinary and partial differential equations[END_REF]. Its ability to cross some types of singularities, its high-order symplecticity, or its high-order iso-spectrality in solving a Lax pair problem have been stated in [START_REF] Deeb | Borel-Laplace summation method used as time integration scheme[END_REF]. Another advantage of BL is that decreasing or increasing the approximation order is as simple as changing the value of a parameter in the code. However, nowhere in the cited works on the Borel-Laplace integrator stiffness has been addressed. One aim of the present article is, as mentioned, to fill this gap.

The Borel-Laplace algorithm that will be discussed here results from the representation of the Borel sum as a Laplace integral. It makes use of a Padé approximation, as will be seen later, and is accordingly named Borel-Padé-Laplace algorithm (BPL). A representation of the Borel sum as an inverse factorial series also leads to an efficient algorithm [START_REF] Deeb | Comparison between Borel-Padé summation and factorial series, as time integration methods[END_REF] but will not be used.

This paper is organized as follows. In section 2, the Borel-Padé-Laplace algorithm is briefly recalled. In section 3, a linear stability analysis is carried out. The stability regions, corresponding to different values of parameters, are plotted. In section 4, numerical performance on stiff and non-stiff ODE problems as well as on a PDE is analyzed.

Borel-Padé-Laplace integrator

Consider an ordinary differential equation or a semi-discretized partial differential equation :

   du dt = F (t, u), u(0) = u 0 , (1)
where

u : [0, T ] -→ R n t -→ u(t)
is the unknown, n ∈ N * is the dimension of the system and F is a non-linear operator

F : R × R n -→ R n (t, v) -→ F (t, v).
Borel-Laplace integrator is based on approximating the solution to (1) via a (convergent or divergent) time series

u(t) = ∞ k=0 u k t k ∈ (C[[t]]) n (2)
and performing a Borel-Laplace summation procedure on this series. In equation (2), C[[t]] stands for the ring of formal power series in t, with complex coefficients. To simplify, assume that n = 1. The terms u k are obtained by inserting directly the series expansion (2) into equation [START_REF] Hairer | Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems, Second Revised[END_REF]. This leads to explicit relations of the form

u k+1 = 1 k + 1 F k (u 0 , • • • , u k ) (3)
where

F k is the k-th Taylor coefficient of F (t, u(t)) at t = 0. It is generally a non-linear function of u 0 , u 1 , • • • , u k .
The expression of F k will be explicitely given for each equation we will be dealing with.

Regarding the validity domain of series (2), there are two possible scenarios. The first one is that the (exact or numerical) radius of convergence of series is zero. In this case, a summation procedure is required. A summation procedure consists in finding an analytic representation (as an integral, a rational function, a continued fraction, . . . ) of the exact solution from the series. The one chosen here is the Borel-Laplace summation in which the solution is represented by a Laplace transform of a rational function (see section 2.1).

The other possibility is that the series is convergent with a finite or theoretically infinite radius of convergence. However, if the series converges slowly, the partial sum

u K (t) = K k=0 u k t k (4)
may give an acceptable approximation only for values of t much smaller than the radius of convergence. In this case, an acceleration of convergence is optional but advisable.

In the present paper, the Borel-Laplace summation will systematically be applied to the series. If the series is convergent, it will act as a convergence acceleration technique. As will be seen, it not only extends the stability region but also increases the speed of the method.

Theoretical setting of Borel-Padé-Laplace summation

The theory behind Borel-Laplace summation can be found in many papers [START_REF] Borel | Leçons sur les séries divergentes[END_REF][START_REF] Ramis | Poincaré et les développements asymptotiques (Première partie)[END_REF][START_REF] Ramis | Les développements asymptotiques après Poincaré : continuité et... divergences[END_REF][START_REF] Costin | Asymptotics and Borel Summability[END_REF] and shall not be reproduced here. Only the computational aspects are presented. Let us assume that series (2) is a p-Gevrey series in a neighborhood of the origin, that is,

|u k | ≤ CA k (k!) p , ∀k ≥ 0 (5)
for some positive real numbers A and C. In fact, it is known that most of series arising in engineering problems are p-Gevrey series for some positive rational number p. In the sequel, we consider only the case p = 1. The numerical summation is done in three stages. First, the Borel transform

B u(ξ) = ∞ k=0 B k ξ k ∈ C[[ξ]] (6)
of series ( 2) is considered. In this expression,

B k = u k+1 k! , k ≥ 0. (7)
Series ( 6) is convergent at the origin. Next, B u(ξ) is prolonged analytically into a function P (ξ) in the vicinity of a semi-line of the complex plane, linking 0 to ∞. Lastly, the Laplace transform (at 1/t), which is the formal inverse of the Borel transform, is applied to the prolonged function. At the end of the procedure, one gets an integral representation S u(t) = u 0 + P (ξ)e -ξ/t d ξ of the solution. The function S u(t) is called the Borel sum of series [START_REF] Iserles | A first course in the numerical analysis of differential equations[END_REF]. The Borel-Laplace summation procedure is summarized in Table 1. If the initial series ( 2) is convergent at the origin, then the Borel sum S u(t) takes the same value as the original sum (2) for t inside the disc of convergence. However, the domain of definition of S u is generally larger than the disc of convergence of the series u. If the initial series is a divergent but Gevrey series, S u is a sectorially analytical function, having the series u as Gevrey asymptotics.

u(t) = +∞ k=0 u k t k ∼ S u(t) = u0 + P (ξ)e -ξ/t d ξ Borel           Laplace B u(ξ) = +∞ k=0 B k ξ k --------------→ Prolongation P (ξ)

Algorithm

Numerically, only a finite number of terms u k can be computed. The series is then represented by the degree K polynomial u(t) u K (t)

u(t) u K (t) = K k=0 u k t k , (8)
that is the partial sum already been defined in [START_REF] Lambert | Numerical Methods for Ordinary Differential Systems: The Initial Value Problem[END_REF]. The Borel transform of u K (t) is a degree (K -1) polynomial. The prolongation is carried out via a Padé approximation [START_REF] Brezinski | Rationnal approximation to formal power serie[END_REF][START_REF] Brezinski | Padé approximations[END_REF]. Other prolongation techniques exist but none of them has been used as part of a Borel-Laplace based time integrator, to the best of our knowledge. A usual Gauss-Laguerre quadrature permits to compute the Laplace transform [START_REF] Stroud | Gaussian quadrature formulas (without numerical tables)[END_REF]. In simulations, the semi-line is the positive real axis. Due to the realization of the prolongation via a Padé aprroximation, the algorithm is called Borel-Padé-Laplace (BPL).

The function Su K (t) obtained with BPL provides an approximate solution to the equation as long as an accuracy criterion is met. When this is no longer the case, the algorithm (computation of u k 's and Borel summation) is restarted using the last acceptable value Su K (t f ) as initial condition. BPL is then a step-by-step method over time. One way of evaluating the accuracy of the approximate solution is to calculate the residue of the equation. This strategy is rather expensive but, as will be seen, is fast enough to compete with all the other numerical schemes under consideration.

The Borel-Padé-Laplace algorithm can be summarized as follows, for a onedimensional problem, with the residue as quality criteria:

1. Start with t 0 = 0 and u 0 = u(t 0 ). 2. Compute the first K coefficients of the series u:

u k+1 = 1 k + 1 F k (u 0 , . . . , u k ), k = 0, • • • , K -1
where F k is the k-th Taylor coefficient of F (t, u(t)) at t = t 0 . 3. Apply the Borel transformation. In other words, compute the first K coefficients of B u:

B k = u k+1 k! , k = 0, • • • , K -1.
B 0 + B 1 ξ + • • • + B K-1 ξ K-1 + O(ξ K ) = a 0 + a 1 ξ + • • • + a Ka ξ Ka 1 + b 1 ξ + • • • + b K b ξ K b =: P (ξ). (9)
5. Obtain the approximate Borel sum by computing the Laplace transform with a N G -point Gauss-Laguerre quadrature formula 1 :

Su K (t) = u 0 + t N G i=1 P (tξ i )w i . (10)
6. Find t f such that the relative residue norm is smaller than a tolerance parameter for all t ≤ t f :

dSu K dt -F t, Su K < Su K (11)
Take Su K (t) as the approximation of u(t) for t ∈ (t 0 , t f ].

1 Note that +∞ 0 P (ξ) e -ξ/t dξ = t +∞ 0 P (tξ) e -ξ dξ 7. Return to step 2 with t 0 = t f , u 0 = Su K (t f ).

In this algorithm, K a and K b are any positive integers such that K a +K b = K-1.

Their influence on the stability region will be analyzed in section 3. The reals ξ i are the roots of the N G -th Gauss-Laguerre polynomials and the w i are the corresponding weights. Note also that a singular value decomposition will be carried out to improve the robustness of the Padé approximation in numerical tests, following an algorithm discussed in [START_REF] Gonnet | Robust Padé approximation via SVD[END_REF].

In stage 6, the final time t f can be determined as follows. First, we evaluate the numerical of convergence τ of the series (2) from its K first terms using a simple criteria developped in [START_REF] Cochelin | A path-following technique via an asymptotic-numerical method[END_REF]. It writes

τ = δ u 1 u K 1/(K-1) (12)
where δ is a small tolerance parameter such that

u K (t) -u K-1 (t) u K (t) -u 0 u K t K u 1 t ≤ δ.
Next, we check if [START_REF] Butcher | Numerical Methods for Ordinary Differential Equations[END_REF] is verified by t = τ . If it is, we try further with 2τ , 4τ , . . . until condition [START_REF] Butcher | Numerical Methods for Ordinary Differential Equations[END_REF] is not verified any longer. If, on the contrary, the first try t = τ is not successful, we try further with τ /2, τ /4 and so on. The last successful value will be taken as t f . The quantity t f -t 0 is considered as the time step of the algorithm. The cut-off order K can be thought as the order of the scheme. Note that one advantage of Borel-Laplace integrator is that, in contrast to many schemes such as BDF or Runge-Kutta, increasing the order is very simple. Increasing K is enough; the algorithm does not need to be modified, no coefficient has to be changed.

In the next section, the linear stability of BPL is analysed. We study in particular the influence of the summation procedure.

Stability analysis

First, we recall the notion of stability domain for a discrete scheme. Consider the scalar linear model problem

   du dt = λu u(t 0 ) = u 0 ( 13
)
where λ is a complex number with a negative real part. The solution of this equation decreases exponentially to zero when t grows. Consider an iterative scheme with a constant time step h, providing discrete approximate solutions v n h u(t n ) of ( 13) at discrete times t n = t 0 + nh as follows:

v n+1 h = R(λ, h)v n h ( 14
)
for some function R of λ and h. The stability domain of this method is defined as the following subset of the complex plane [START_REF] Hairer | Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems, Second Revised[END_REF][START_REF] Iserles | A first course in the numerical analysis of differential equations[END_REF]:

D = { (λh) ∈ C : |R(λ, h)| < 1 }. ( 15
)
When the time step h is such that λh lies in the stability region, the approximate solution decreases to zero, like the exact one, when n grows. BPL is not a discrete scheme, in the sense that it does not provide a discrete approximation of the solution but a continuous one. It is however relatively easy to adapt to it the notion of stability domain. For simplicity, assume that t 0 = 0. For equation ( 13), we have:

F (t, u) = λu and F k (u 0 , . . . , u k ) = λu k . (16)
Equation ( 3) permits to compute the coefficients of the time series:

u k+1 = λu k k + 1 . (17)
When inserted into series (2), these coefficients lead of course to the Taylor expansion of the exact solution e λt . It is a convergent series. We carry out two different linear stability investigations. The first one is when the solution is approximated by the Taylor series truncated at order K, without the Borel summation procedure, and the second one is when the solution is approximated with the BPL scheme. When the summation procedure is not applied, the method will be called time Asymptotic Numerical Method (ANM) as in computational solid and fluid mechanics [START_REF] Cochelin | Méthode asymptotique numérique, Methodes numériques[END_REF]. With ANM, we have:

u(h) = K k=0 (hλ) k k! u 0 . (18)
Comparing this relation to ( 14), we define the domain of linear stability, for a given truncation order K, as

D K AN M = z ∈ C such that K k=0 z k k! ≤ 1 . ( 19
)
This domain is plotted in Figure 1 for K from 2 to 10. As can be observed, D K AN M grows with K. The growth is however rather slow. Let us use the positive number |D K AN M | defined as follows as a quantification of the size of This quantity increases almost linearly as can be seen in Figure 1. The slope of the curve is about 0.375. In fact, the region D K AN M coincides with the stability region of an explicit K-th order Runge-Kutta method. The reason to this is that the stability function of an explicit Runge-Kutta method is a truncated Taylor expansion of the exponential function.

D K AN M : |D K AN M | = sup d ≥ 0 such that [-d, 0] ∈ D K AN M . (20)
We now apply the summation procedure to the series. To make it more concrete, let us set u 0 = 1 and K = 4. Like previously, the truncated series solution is

4 k=0 (λt) k k! = 1 + λt 1 + (λt) 2 2 + (λt) 3 6 + (λt) 4 24 . (21)
Its Borel transform reads

3 k=0 λ k+1 (k + 1)! ξ k k! = λ 1 + λξ 2 + (λξ) 2 12 + (λξ) 3 144 (22)
We take k a = 1 and k b = 2. The [1/2] Padé approximant of ( 22) is

P (ξ) = λ 48 + 14λξ 48 -10λξ + (λξ) 2 . ( 23
)
We then have the following approximate solution of the linear equation ( 13) with the Borel-Padé-Laplace scheme

u(t) 1 + t N G i=1 P (tξ i )ω i = 1 + λt N G i=1 Q(λtξ i )ω i ( 24
)
where Q is the rational function

Q(z) = 48 + 14z 48 -10z + z 2 . ( 25
)
As already mentioned, ξ i is the i-th root of the N G -th Laguerre polynomial and ω i is the corresponding weight in Gauss-Laguerre quadrature. Note that the Padé approximant [START_REF] Costin | Borel summability of Navier-Stokes equation in R 3 and small time existence[END_REF] has no pole on the integration domain (the real positif axis) of the Laplace integral.

The stability region of Borel-Padé-Laplace integrator, for K = 4, is

D 4 BP L = z ∈ C such that 1 + z N G i=0 Q(zξ i )ω i ≤ 1 . ( 26
)
This region is plotted in Figure 3, with N G = 100 Gauss points, along with the stability regions for other values of K, ranging from 2 to 10. In this Figure, the Padé approximant in Borel space is chosen as closed as possible to the diagonal, that is

K b = K a or K b = K a + 1 (27)
depending on the parity of K. This choice will be discussed later.

As can be seen in Figure 3, the regions do not include the half complex plane with negative real part. Indeed, as an explicit scheme, BPL is not Astable. However, this figure clearly shows that, for a fixed K ≥ 4, the stability region of BPL is much wider than that of the simple truncated series scheme or that of an explicit Runge-Kutta scheme. Note that when K = 2, the Borel summation has no effect, and the stability region is the same as in Figure 1a.

Another striking point is that the growth of the stability region with K is not as regular as in the case where the Borel summation is not applied. This is due to the choice of (almost) diagonal Padé approximant. However, if we consider either only odd K or only even K, the growth is regular again. In all cases, the overall growth rate is higher than in Figure (1a). Indeed, let us quantify the size of D K BP L as previously with

|D K BP L | = sup d ≥ 0 such that [-d, 0] ∈ D K BP L (28)
The dependence of |D K BP L | on K is plotted in Figure 4 (along with the evolution of |D K AN M | for comparison), for K ranging from 2 to 12. The average slope of the curve of |D K BP L | is about 0.543. Since we are in the particular situation where the Taylor coefficients u k of the solution decrease very rapidly with k, only few terms of the series are numerically meaningful. More precisely, the coefficients B k of the Borel transformed series B u are below our machine precision (about 2 • 10 -16 ) for k > 11. So, taking K ≥ 12 does not bring any substantial improvement.

These observations indicate the importance of the Borel summation procedure, even in a situation where it has not been developed for. Indeed, the summation procedure enlarges very significantly the stability region, even when the Taylor series of the solution is convergent, with an infinite radius of convergence. This stability region can even be larger if we play with the parameters of the Padé approximants in Borel space. For example, we plot in Figure 5 the evolution of D K BP L with K when K a is fixed to 1. As can be observed, the stability domain grows with K and is almost always far larger compared to Figures 1 and3. The growth rate is also far higher.

To end up, we would like to analyse graphically the influence of K a (or K b ) when K is fixed. The stability regions corresponding to K = 10 and different Padé degrees are plotted in Figure 6. It can be observed that the stability region grows with the degree K b of the Padé denominator. When K a = 0 (Figure 6j), BPL tends to be A(α)-stable for some angle α. A theoretical study on the optimal choice of K a and K b , which take into account the stability and the precision, would be very interesting but has not been carried out yet. In the sequel, a (almost-) diagonal Padé approximant satisfying relation [START_REF] Thomann | Formal and numerical summation of formal power series solutions of ODE's[END_REF] will be chosen. This is motivated by some good properties of diagonal Padé (convergence [START_REF] Baker | An investigation of the applicability of the Padé approximant method[END_REF][START_REF] Nuttall | The convergence of Padé approximants of meromorphic functions[END_REF][START_REF] Baker | Defects and the convergence of Padé approximants[END_REF], invariance under linear fractional transformation [START_REF] Baker | Essentials of Padé Approximants[END_REF], A-stability of diagonal Padé approximants to the exponential function [START_REF] Hairer | Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems, Second Revised[END_REF], . . . ). As seen, it may not correspond to an optimal choice but it will be shown that it is good enough to obtain a very competitive performance in terms of computation time.

Note that some of the plotted stability regions are not complete. Indeed, the whole stability regions may contain other parts in the complex plane, but these parts have been excluded from the graphics.

In the next section, the performance of BPL in solving stiff and non-stiff equations is analysed. As mentioned, all the previous figures were plotted with N G = 100 Gauss points. For N G = 20 and for N G = 200, only small changes

Numerical performance

Unless otherwise stated, the order K of BPL is set to 10. The degrees of the numerator and the denominator of the Padé approximant ( 9) are K a = 4 and K b = 5.

We compare BPL with some classical numerical schemes. Some of them are popular choices for solving stiff equations. These schemes have either a forth or a tenth consistency order.

Classical schemes

The following schemes are considered.

• The 4-stage 4-th order explicit Runge-Kutta algorithm with a Fehlberg adaptive time step [START_REF] Fehlberg | Klassische Runge-Kutta-Formeln vierter und niedrigerer Ordnung mit Schrittweiten-Kontrolle und ihre Anwendung auf Wärmeleitungsprobleme[END_REF][START_REF] Hairer | Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems, Second Revised[END_REF], called RK4 hereinafter.

• The 5-stage 10-th order implicit adaptive Gauss-Legendre method which is a Runge-Kutta scheme combined with a Gauss-Legendre quadrature [START_REF] Hairer | Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems, Second Revised[END_REF], refered as GAU.

• The 4-step 4-th order implicit backward differentiation formula [START_REF] Curtiss | Integration of stiff equations[END_REF][START_REF] Hairer | Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems, Second Revised[END_REF], initialized with RK4, and named BDF in this article.

• The 4-th order exponential time differencing method combined with the adaptive Runge-Kutta-Fehlberg method [START_REF] Cox | Exponential time differencing for stiff systems[END_REF]. This method is generally called ETDRK4, but will simply be shortened to ETD.

RK4 has been chosen for its popularity and speed in solving non-stiff equations, GAU for its order 10 (the same order as that set for BPL), BDF for its popularity in solving stiff equations and ETD because it is a relatively recent integrator for stiff equations. All of these methods are adaptive. The step size is updated at each time iteration with a formula

h n+1 = 0.9 h n τ e n+1 1 k+1
where k is the order of the scheme and τ is a small tolerance parameter which can be choosen to adjust the accuracy of the method. e n+1 is an estimation of the local error. Indications on how it is computed are given below for each scheme.

For RK4 and GAU, which are multi-stage one-step methods, the approximate solution at t = t n+1 can be written as follows:

u n+1 = u n + h n (b 1 k 1 + • • • + b s k s )
where s is the number of stages. The intermediate values k i and the coefficients b i are defined in equation (1.8) and in Table 5.1 of [START_REF] Hairer | Solving Ordinary Differential Equations I: Nonstiff Problems, Second Revised[END_REF] for RK4, and in equation (7.7) of [START_REF] Hairer | Solving Ordinary Differential Equations I: Nonstiff Problems, Second Revised[END_REF] and in page 71 of [START_REF] Hairer | Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems, Second Revised[END_REF] for GAU. The estimated error e n+1 is obtained from the difference between u n+1 and a second estimation

u * n+1 = u n + h n (b * 1 k 1 + • • • + b * s k s )
of the approximate solution, having at least an order k of constistency. The coefficients b * i are provided in Table 5.1 of [START_REF] Hairer | Solving Ordinary Differential Equations I: Nonstiff Problems, Second Revised[END_REF] for RK4. They are defined in [START_REF] Hairer | Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems, Second Revised[END_REF], equation (8.16), and in [START_REF] De Swart | On the construction of error estimators for implicit Runge-Kutta methods[END_REF], section 2, for GAU.

For BDF, the approximate solution is determined from a relation

a n+1 1 u n+1 + • • • + a n+1 s u n+1-s = h n f (t n+1 , u n+1 ) (29)
where s = 4 is the step number. The coefficients a n+1 i , in the adaptive case, are described in appendix G of [START_REF] Co | Methods of Applied Mathematics for Engineers Scientists[END_REF], subsection G4. The estimation of the local error is explained in pages 372-373, in Theorem 6.2 and in Table 6.2 of [START_REF] Hairer | Solving Ordinary Differential Equations I: Nonstiff Problems, Second Revised[END_REF].

Lastly, if equation ( 1) is a scalar ODE then the ETD approximate solution is defined in [START_REF] Cox | Exponential time differencing for stiff systems[END_REF], equation [START_REF] Deeb | Borel-Laplace summation method used as time integration scheme[END_REF]. When equation ( 1) is not scalar, a pseudoinversion and an exponentiation of a matrix is needed. They are carried out respectively with a singular value decomposition and a matrix Padé approximation. The local error is obtained from a perturbation of equation ( 29) of [START_REF] Cox | Exponential time differencing for stiff systems[END_REF].

All the schemes are implemented entirely in python with a fairly equal effort in optimization. The computations are run on a single processor. We focus on accuracy, the size of time step and computation time.

We now apply these schemes to some classes of differential problems.

Lotka-Volterra equations

Consider a prey-predator system, dynamically governed by the Lotka-Volterra equations [START_REF] Hofbauer | The Theory of Evolution and Dynamical Systems: Mathematical Aspects of Selection[END_REF]:

       du dt = α u -β uv, dv dt = -δ v + γ uv, (30)
where u and v are respectively the number of preys and predators in the population, and α, β, δ, γ are real positive constants. The reproduction parameter α is the natural (exponential) growth rate of preys in absence of predators whereas δ is the natural decline rate of predators in absence of preys. βv is the mortality rate of prey depending on the the number v of predators and δu is the birth rate of predators depending on the number of prey eaten. It is straight forward to show that system [START_REF] Deeb | Comparison between Borel-Padé summation and factorial series, as time integration methods[END_REF] possesses the first integral:

I(u, v) = βv + γu -α ln v -δ ln u. ( 31
)
We first choose a set of coefficients for which the problem is not stiff.

Non-stiff case

Take an initial population which consists of two preys and one predator, that is u 0 = 2 and v 0 = 1. The reproduction/decline parameters are set to α = 2/3 and δ = 2 and the predation parameters to β = 4/3 and γ = 2.

For BPL, the function F k which operates in the recurrence relation ( 3) is defined by

F k (u 0 , . . . , u k , v 0 , . . . , v k ) =       αu k + β k l=0 u l v k-l -δv k + γ k l=0 u l v k-l       . ( 32
)
The parameter of BPL which is used in the accuracy criterion ( 11) is set such that the mean error on the first integral ( 31) is about 1.35•10 -7 over a simulation time T = 1000. This mean or overall error is defined as an approximation of

1 T T 0 I u(t), v(t) -I u(0), v(0) dt. ( 33
)
The approximate solution over 40 seconds is visualized in Figure 7. It necessitated 254 iterations. To obtain the smooth plots in Figure 7, not only the value of the solution at the discrete times (t i ) i=0,...,254 but also at some intermediate times t ∈]t i , t i+1 [ are plotted. In contrast to many other schemes, no interpolation method is needed for this. Formula ( 10) directly provides the approximate solution within each interval ]t i , t i+1 [. The accuracy parameter τ are set for each step (BDF, ETD, GAU and RK4) such that their a posteriori accuracies are comparable to that of BPL. The mean errors are reported in Table 2. They are around 4 • 10 -7 .

Figure 8a shows the evolution of the time steps of the different methods. The solution being periodic, only the evolution over the last 100 seconds are plotted. As can be seen, it is with BPL that the time step is the largest. Since we are in a non-stiff case, the classical Runge-Kutta method has a good performance and competes with the 10-th order Gauss scheme in terms of time step. Figure 8b represents the same data as Figure 8a but with a logarithmic scale in ordinate. It shows that the time step of BPL is about 60 times larger than that of BDF and about 80 times as large as that of ETD. It is confirmed in Table 2 which compares the mean values. As for CPU time, the two explicit integrators, BPL and RK4, have a comparable performance (see Table 2). They need about 10 times less computation time than GAU and at least 100 times less than BDF and ETD.

In a second test, each scheme is run with multiple values of the (residue or estimated error) tolerance. The mean time step is plotted in Figure 9 against the overall accuracy. This figure clearly shows that, amongst the considered schemes, BPL has always the largest mean time step, whatever the precision. This mean time step is about 6.6 times as large as that of the Gauss scheme with the same order, for an error around 3 • 10 -9 . This large time step results in a faster computation. Indeed, as can be noticed in Figure 10, BPL requires much less CPU time than GAU, for comparable precisions. Only RK4 is faster than BPL for a medium or a low precision. But when a high precision is required, BPL tends to be more interesting.

Increasing the stiffness ratio

We now examine the behaviour of the schemes when the stiffness ratio varies. The stiffness ratio r is defined as the spectral condition number of the linear part of equations [START_REF] Deeb | Comparison between Borel-Padé summation and factorial series, as time integration methods[END_REF], that is

r = max(α, δ) min(α, δ) (34)
since α and δ are positive real numbers. For the sake of simplicity, and since it will be the case in the numerical experiments, assume that δ > α, such that r = δ/α. In fact, r appears naturally when equations [START_REF] Deeb | Comparison between Borel-Padé summation and factorial series, as time integration methods[END_REF] are adimensionalized with the variables

u * = γu δ , v * = βv α , t * = αt. (35)
Indeed, equations (30) can be written as follows:

       du * dt * = u * (1 -v * ), dv * dt * = rv * (-1 + u * ). (36)
All the parameters of the equations are kept at the same value as before, except δ which is increased. As before, simulations are run with multiple values of the (residue or estimated error) tolerance until 1000 seconds. The error on the first integral and the CPU time are recorded and plotted hereafter.

For a moderate stiffness ratio r = 8, Figure 11 shows that RK4 and BPL compete in terms of CPU time, even if BPL indicates a slight advantage for high precision simulations. It can also be stated in this figure that GAU can provide a very accurate solution, due to its high order, but with a higher cost than BPL. BDF and ETD are much more expensive than the other schemes.

For r = 16, we have approximately the same picture, except that BPL becomes more interesting than RK4 even for moderate precisions. This can be observed in Figure 12. When the stiffness ratio is set to a high value r = 32, the situation changes significantly. First, as can be seen in Figure 13, RK4 cannot reach very high precision any longer, compared to BPL. The precision that can be achieved with GAU is still very high but not as high as with r = 16. ETD also looses precision.

Only BPL is able to maintain the same precision as previously.

Concerning the numerical cost, the increase of CPU time needed by BPL and ETD is very small compared to that of GAU.

Lastly, BDF does not appear in Figure 13. Indeed, although it is a popular method for stiff equations, it fails with r = 32. It diverges as soon as t reaches few seconds. This behaviour has also been observed with the optimized BDF solver of the python scipy package, with the optimized BDF solver of Scilab, and with the option CVODE BDF of the package Sundials of Julia language.

For r = 64, ETD also fails. As remarked in Figure 14, RK4 gives moderately accurate solutions, and even wrong solutions for some values of the predicted error tolerance. Indeed, even for very small value of the tolerance, the overall error may be larger than one. Figure 14 also shows that the precision of GAU seems to stagnate around 2.6 • 10 -4 . Only BPL can provide highly accurate solutions. Moreover, its CPU cost is very small compared to that of GAU.

At last, with r = 128, the Gauss method also fails. This behaviour has as well been observed with the Gauss solver of the Matlab package numeric::odesolve. For this value of the stiffness ratio, RK4 cannot give an accurate solution any longer, whereas with BPL, the error can be as small as 7.4•10 -10 (see Figure 15).

These numerical experiments shows that BPL is an interesting alternative method for stiff problems. First, its arbitrary high order allows to get highly accurate solutions. With Lotka-Volterra equations, it never fails for values of r up to 128. Moreover, its cost is generally much smaller than that of the other methods, due to its explicit property.

The previous tests show the performance of BPL for the resolution of stiff and non-stiff ODE's. In the next subsection, we examine its efficiency in solving partial differential equations.

Korteweg-de-Vries equation

In this subsection, we consider the Korteweg-de-Vries equation (KdV)

∂u ∂t + c 0 ∂u ∂x + β ∂ 3 u ∂x 3 + α 2 ∂u 2 ∂x = 0 (37)
which models waves on shallow water surfaces [START_REF] Korteweg | On the change of form of long waves advancing in a rectangular canal, and on a new type of long stationary waves[END_REF]. In this equation, the linear propagation velocity c 0 , the non-linear coefficient α and the dispersion coefficient β are positive constants, linked to the gravity acceleration g and the mean depth d of the water by:

c 0 = gd, α = 3 2 g d , β = d 2 c 0 6 . (38)
In order to focus on the performance of the time integrators, we choose a high order scheme, namely a spectral method, for the space discretization. The solution is assumed to be periodic with period X in space, and integrable. It is approximated by its truncated Fourier series:

u(x, t) |m|≤M ûm (t) e imωx , (39)
where M ∈ N and ω = 2π X . The substitution of equation ( 39) into (37) leads to a (2M + 1)-dimensional ODE

dû dt = Aû + N (û) (40)
where the array û contains the unknowns ûm , A is a diagonal matrix with diagonal entries

A m m = -c 0 iωm + iβω 3 m 3 (41)
and N (û) is a non-linear array containing convolution terms:

N (û) = - 1 2 iαmω û * û. ( 42
)
Convolution operations are performed in physical space and the standard dealiasing 3/2 rule is applied. With BPL, each component ûm (t) of the Fourier coefficient array û(t) is decomposed into its Taylor series

ûm (t) = K k=0 ûm k t k . ( 43
)
The series coefficients are computed explicitely as follows:

ûk+1 = 1 k + 1 (-c 0 iωm + iβω 3 m 3 )û k - 1 2 iαmω k l=0 ûn * ûk-l . (44)
For the simulations, the initial condition is the periodic prolongation of the function

u 0 (x) = U sech 2 (κx), x ∈ - X 2 , X 2 , (45)
U being a constant and κ = 3U 4d 3 . The corresponding exact solution is the traveling wave u(x, t) = u 0 (x -ct).

with c = c 0 1 + U 2d . We take X = 24π, d = 2, g = 10 and U = 1 2 . The solution is periodic in time, with a period T 14.986s.

We use D = 2M to indicate the size of the system, instead of the dimension 2M + 1 of equation [START_REF]Automatic differentiation: applications, theory, and implementations[END_REF]. The simulations are run over one period, for some values of D between 64 and 512. It is hard to calibrate the tolerance parameter τ of all the schemes to have the same (a posteriori) overall error at each value of D. So, this calibration has not been done. Instead, we require more accuracy to BPL than to the other methods, in order not to overestimate the performance of BPL. The overall error are recorded in Table 3. The error reported in this table is an approximation of

T 0 u computed (t) -u exact (t) u exact (t) dt. ( 47
)
The evolution of the computation time of each scheme is plotted in Figure 16. It can be seen there that BDF requires a very high cost for D = 128, despite the low precision (see second column of Table 3). As a consequence, it has not been used for higher values of D.

Figure 16 also shows that, among the considered schemes, RK4 is the fastest for (non-stiff) small-sized problems, for the given precisions. But for high degrees of freedom, BPL becomes the most interesting in terms of computational time. BPL also has the smallest slope.

Figure 17 indicates that BPL has a very large mean time step compared to the other schemes, whatever the size of the problem is. It is also striking that the mean time step does not vary very much with the size of the problem. However, BPL is not the only scheme which presents this characteristics since ETD exhibits the same behavior, but with much smaller time steps.

The large time step of BPL is of a great importance in its performance. Indeed, the CPU time spent at each time step is very high with BPL in comparison to the other schemes, as can be stated in Figure 18. One reason for this is the evaluation of the residue in step 6 of the algorithm presented in section 2.2. This evaluation is done multiple times at each time step to decide if the solution is still accurate enough. Another precision evaluation is desirable, but not available yet. Fortunately, this expensive precision evaluation is largely counter-balanced by large time steps.

In all of the previous simulations, the order K of the time series in BPL was set to 10. In our last test, the effect of K on the performance of BPL is analysed. For this, the size of the problem is set to D = 128. A residue tolerance = 1 • 10 -4 is chosen. Figure 19 shows the L 1 relative error (defined in equation ( 47)) over one period. This figure reveals a fluctuation of the error according to the parity of K. Note that such fluctuation is not uncommon when manipulating truncated series. Moreover, the parity of K intervenes in the choice of the Padé approximants in Borel space. Indeed, when K is odd, the numerator and the denominator of the Padé approximant have the same degree; and when K is even, the denominator has a higher degree than the numerator (see choice in equation ( 27)). Figure 19 however tells us that globally, the accuracy increases with the order K of the series, for a fixed value of the residue tolerance.

The mean time step has also a globally increasing tendency with K as can be seen in Figure 20, passing from ∆t mean = 0.0256 for K = 4 to ∆t mean = 0.156 when K = 14. As a consequence, the CPU time decreases with K, as can be noted in Figure 21. These results tend to indicate that high values of K accelerate the computation.

Conclusion

In this article, we studied the linear stability of the Borel-Padé-Laplace integrator. It has been shown that if the summation procedure is not applied, the scheme has the same linear stability domain as an explicit Runge-Kutta integrator. But when the summation is carried out, the linear stability domain enlarges very significantly, even when the Taylor series of the solution is convergent. It has been obsereved that the size of this domain increases with the truncation order of the series. We also saw that the choice of Padé approximants in Borel space has a substantial impact on the size of the linear stability domain.

Even if BPL is not A-stable, it has been shown that this scheme is more efficient than many explicit and implicit ones, in solving stiff problems. It runs without any particular difficulty for a wide range of values of the stiffness ratio. Due to its high order, it can reach very high precisions even when the stiffness number is high. Moreover, its explicit property makes it very fast compared to the other integrators.

Numerical tests on non-stiff Lotka-Volterra and on Korteweg-de-Vries equations showed that for small-size systems, the popular 4-th order Runge-Kutta method has a comparable speed than BPL when only a moderate precision is needed. But when high precision is required BPL becomes more interesting. It is even more true when the size of the system is large.

It is worth to notice that increasing the approximation order of BPL does not require any programming effort. One has simply to raise the cut-off parameter K of the series, without changing anything else in the algorithm. As could be observed in the last part of the article, the higher this value is, the faster BPL is.

Despite its speed, one optimization should be brought to the algorithm of BPL. Indeed, a numerical test with Korteweg-de-Vries equation showed that a BPL time step is rather expensive, due among others to many evaluations of the residue. A more efficient accuracy estimation should be developed. This should increase the speed of the scheme.

To obtain the previous results, the computation was done on a single processor. Note however that BPL also presents some advantage regarding parallelization. Indeed, the summation algorithm can be done component-wise, letting the computation to be shared between many processors.

In this paper, only the computational aspects of BPL are discussed. The founding theory was skipped. Yet, some optimizations may be brought to the algorithm with help of theoretical considerations. For instance, a theoretical study of the equation may be helpful to determine the actual Gevrey order (which was set to one in this article). However, it is conceivable to evaluate numerically this Gevrey order from the coefficients of the series. A theoretical study of the equation may also help to find a better (than the real positive semi-line) integration direction in the Laplace transform.

4 .

 4 Compute a [K a /K b ] Padé approximant P (ξ) of the polynomial with coefficients B k , i.e. determine a k and b k such that

Figure 1 :Figure 2 :

 12 Figure 1: Linear stability regions of the truncated Taylor series approximation (without Borel summation)

  K = 2, Padé [0/1]

  K = 3, Padé [1/1]

  K = 10, Padé[4/5]

Figure 3 :

 3 Figure 3: Linear stability regions of Borel-Padé-Laplace integrator with increasing K and close-to-diagonal Padé approximants.

Figure 4 :

 4 Figure 4: Evolution of |D K AN M | and |D K BP L | with K

  K = 2, Padé [1/0]

  K = 3, Padé [1/1]

Figure 5 :

 5 Figure 5: Linear stability regions of Borel-Padé-Laplace integrator with increasing K and fixed Ka = 1

Figure 6 :

 6 Figure 6: Linear stability regions of Borel-Padé-Laplace integrator for fixed K = 10 and decreasing Ka

  Trajectory in (u, v) plane

Figure 7 :

 7 Figure 7: Approximate solution with BPL

Figure 8 :Figure 9 :

 89 Figure 8: Non-stiff Lotka-Volterra. Evolution of time step

Figure 10 :

 10 Figure 10: Non-stiff Lotka-Volterra. Evolution of CPU with the mean error

Figure 11 : 8 Figure 12 :

 11812 Figure 11: Lotka-Volterra. Stiffness ratio r = 8

Figure 13 :

 13 Figure 13: Lotka-Volterra. Stiffness ratio r = 32

Figure 14 :

 14 Figure 14: Lotka-Volterra. Stiffness ratio r = 64

Figure 15 :

 15 Figure 15: Lotka-Volterra. Stiffness ratio r = 128

  • 10 -1 3.71 • 10 -4 2.92 • 10 -3 5.11 • 10 -3 1.83 • 10 -3 128 1.08 • 10 -2 3.54 • 10 -4 3.27 • 10 -3 5.81 • 10 -3 1.69 • 10 -3 256 -3.61 • 10 -4 3.66 • 10 -3 4.00 • 10 -3 1.23 • 10 -3 512 -3.17 • 10 -4 4.11 • 10 -3 2.65 • 10 -3 6.50 • 10 -4

Figure 16 :

 16 Figure 16: KdV. Evolution of the computation time with the size D of the problem

Figure 17 :

 17 Figure 17: KdV. Evolution of the mean time step with the size D of the problem

Figure 18 :

 18 Figure 18: KdV. Mean CPU time per time step

Figure 19 :Figure 20 :Figure 21 :

 192021 Figure 19: KdV. Evolution of the error with K

Table 1 :

 1 Borel-Laplace summation

Table 2 :

 2 Mean time step 2.42 • 10 -3 1.65 • 10 -1 2.07 • 10 -4 3.70 • 10 -2 3.10 • 10 -2 Error on the first integral and CPU time

		BDF	BPL	ETD	GAU	RK4
	Mean error	5.56 • 10 -7 1.35 • 10 -7 7.22 • 10 -7	4.6 • 10 -7	2.38 • 10 -7
	CPU	5.96 • 10 2	4.32	9.34 • 10 2	4.50 • 10 1	4.96

Table 3 :

 3 KdV. Overall error

		10 5 10 6		BDF BPL ETD GAU RK4
		10 4
	CPU	10 3
		10 2
		10 1
		100	200	300 D	400	500