Mathias J'ai

Tout d'abord, je remercie le bon Dieu, tout puissant, de m'avoir donné la force pour survivre, ainsi que l'audace pour dépasser toutes les difficultés et mener à bien ce travail. Ce travail n'aurait jamais vu le jour sans l'appui d'un certain nombre de personnes que j'aimerais remercier ici.

Je tiens en tout premier lieu à remercier mon directeur de thèse Monsieur Jean-Marc Cadou et mon co-encadrant Monsieur Yann Guevel pour la confiance qu'ils m'ont accordée en me proposant ce sujet. Même s'il y a eu de nombreux moments difficiles durant les trois ans, je n'aurais jamais pu y aller jusqu'au bout, sans leur soutien, leur appui scientifique, leurs outils pédagogiques, et surtout leur côté humain. Je n'oublierai jamais le bénéfice qu'ils ont apporté sur le développement et l'amélioration de ma réflexion en mathématiques pures vers le monde numérique. Je ne trouve pas les mots pour définir à quel point je leur suis reconnaissante. Pendant ces trois années, j'ai vraiment eu l'impression de faire partie non seulement d'une équipe, mais aussi d'une famille. Je les remercie de m'avoir permis de réaliser cette thèse dans des conditions optimales.

Je remercie particulièrement Monsieur Aziz Hamdouni qui m'a fait l'honneur de présider mon jury de thèse. J'apprécie ses remarques et ses conseils sur la méthode de resommation. Je remercie également Monsieur Bruno Cochelin, Monsieur Mejdi Azaiez pour leur rôle de rapporteurs et Madame Marianne Beringhier pour son rôle d'examinatrice. Je les remercie pour l'intérêt qu'ils ont porté à mon travail de thèse et pour leurs remarques constructives qui ont été bien prises en compte. Je suis fière d'avoir partagé mon travail avec eux. Je remercie également Monsieur Nicolas Jacques et Monsieur Cyrille Allery d'avoir accepté d'être les membres du CSI (comité de suivi individuel) pour ces trois années de thèse.

Je tiens à remercier tous les autres membres du groupe, en particulier Laetitia, Grégory et Khaoula. Je n'oublierai pas les moments partagés au bureau avec Nejia, Diane

n'a pas été facile pour moi. Je remercie également les autres doctorants : Xavier, Adrien, Alemeh, Nora, Malik, Delphine, Ghita, Thomas, Miguel, Maxime et Julien pour tous les bons moments partagés durant ces trois années. Je n'oublierai pas Amine, mon partenaire de déjeuner tous les jours au RU. Je le remercie pour toutes les pauses-café, les différents cadeaux, les chocolats, les belles surprises, et les différents conseils administratifs et moraux. Je remercie également Asmahan, avec elle, nous nous sommes soutenus dans les bons et dans les mauvais moments (notamment le COVID) et je suis immensément fière du chemin que nous avons parcouru ensemble.

Mes pensées vont également à Bertrand qui a sauvé mes fichiers et mon ordinateur plusieurs fois. Merci également à Béatrice, Virginie, Noluenn, et Roseline pour leurs aides. Je remercie aussi Ahmad pour nos différents échanges scientifiques à la rochelle sur les méthodes de resommation. Je tiens à remercier mes autres amis libanais en France et à l'étranger qui m'ont permis de vivre des moments de bonheur durant ces 3 années : Freddy, Marwan, Hassan, Marianne, Grace, Joe, Ghewa, Zeinab, Jihane, Maria, Abbas, Patrick, Marie-Noël, Jad, Machhour, Marinelle, Elie, Joseph, Mohamad, Tatiana et Romaine. J'ai une pensée particulière pour Marie-Noël, avec elle, j'ai partagé des jours d'aventure en Bretagne que je n'oublierai jamais. Je tiens à remercier Kassem pour les différentes soirées à la résidence qui m'ont donné l'impression d'avoir une seconde famille en Bretagne. Je sais que cela a été dur pour lui de me supporter dans mes moments stressants et de m'accompagner durant ces années. Je n'ai pas de mots pour lui dire combien je suis reconnaissante, 1000 merci Kassem. Merci aussi pour le chef Hamza pour tous les plats et les délicieux desserts qu'il a préparé dans la résidence. Merci également à Assi pour toutes les promenades en voiture, les dîners ensemble, l'aide apportée pour faire les courses, et pour avoir été là quand je me sentais seule.

Et enfin, mes remerciements les plus sincères vont à mes parents qui m'ont toujours encouragé à poursuivre mes études à l'étranger. Je leur remercie pour leur compréhension et leur soutien pendant les moments difficiles à l'étranger. Je ne serai jamais assez reconnaissante pour tout ce que vous avez fait pour moi depuis le début. Sans vous, je n'en serais pas là ! Je n'oublierai jamais à remercier ma soeur Nancy et mon frère Elie pour les appels et les visios chaque jour pour me soutenir. Je réalise vraiment ce que signifie avoir un frère et une soeur malgré la distance de 3155 km qui nous sépare.

À mon père et ma mère,

TABLE OF CONTENTS INTRODUCTION

This thesis is related to the numerical solution of time-dependent problems. In this context, consider an ordinary differential or a partial derivative equation in the following form:

du dt (t, x) = F (t, u) u, x ∈ R n ( 1
)
where u is the unknown and F is a nonlinear operator.

The possible resolution of such a problem can be analytical or numerical. The analytical solution does not exist in some cases, and it is difficult to find in others. This is especially true for problems with multiple degrees of freedom. As a result, a numerical method is employed. The usual numerical schemes for evolution problems are based on time discretization and, if the problem is modeled by a partial differential equation, a spatial discretization is then used. Established approaches for spatial discretization include the finite element method (FEM), the finite difference method (FDM), and the finite volume method (FVM). These latter aim to find approximate solutions of partial differential equations, in a set of points rather than in the whole domain. Concerning time discretization, the first numerical method for approximating the first-order derivative in time, which is du dt (t, x) , was introduced in 1768 by Leonhard Euler who treated it in his book [START_REF] Euler | Institutionum Calculi Integralis[END_REF]. Since then, various techniques have been developed for different orders, which are based on the temporal discretization of the problem to be solved, by advancing the model by successive discrete-time steps denoted by ∆t. Depending on the type of formula used to approximate the solutions, a distinction is made between single and multi-step schemes, and between explicit and implicit schemes. Explicit methods calculate the state of a system at a later time from the current state, whereas implicit methods find a solution by solving an equation involving both the current and later states of the system. For example, Eq.( 1) can be discretized using the simplest explicit and implicit methods, which are the forward and backward Euler methods, in the following forms:

u(t + ∆t) = u(t) + ∆t F (u(t)) Forward Euler scheme (2)
u(t + ∆t) = u(t) + ∆t F (u(t + ∆t)) Backward Euler scheme

The explicit method is simpler to program and takes less time to compute. However, it has a low range of stability. Therefore, a small enough step size is required to avoid divergence. On the other hand, the implicit method has high stability and converges when the appro-priate parameters are set. They are used because many problems encountered in practice are stiff, requiring impractically small time steps to keep the error bounded. However, the main drawback of implicit time integrators is their high computational cost compared to explicit methods. This is due to the linear system resolution. Nevertheless, it is important to note that explicit integrators are more suitable for computing unsteady phenomena such as shocks and acoustic waves. One may cite the common time integration methods used in the literature which are Euler, Crank-Nicholson, Newmark, Runge-Kutta, Adams-Bashforth, Adams-Moulton, Backward differentiation formulas, Heun, etc... A review of several time integrators of different type can be found in [START_REF] Muscat | Coupling of time integration schemes for compressible unsteady flows[END_REF]. The appropriate use of a numerical scheme for a given equation is evaluated using three criteria proposed in numerical analysis. The first criterion is the convergence which guarantees that the approximated solution is close to the exact solution. The second one is the order which quantifies the quality of the approximation. The last one is stability, which shows the behavior of the error. In many cases, when solving evolving problems, these criteria or one of them are no longer respected. For example, the classical integration schemes might be unstable and the third criterion is no longer satisfied. Thus, the choice of their parameters is an essential and delicate issue. For explicit schemes, beyond the critical time step related to the Courant-Friedrichs-Lewy (CFL) condition [START_REF] Courant | Über die partiellen Differenzengleichungen der mathematischen Physik[END_REF], solving the problem by classical discretization techniques becomes imprecise. This leads to a destruction in the invariance properties of the physical phenomena, thus the conservation of energy and momentum is no longer respected. With a time step smaller than the critical one, it may be possible to respect the physical phenomenon, but on the other hand, the number of iterations and the computation time may increase. Another example to highlight the limitations of classical time integration schemes is the numerical resolution of coupled problems, such as those of Fluid-Structure Interaction (FSI). In these cases, both the structure and fluid domains require the use of a distinct critical time step [START_REF] Van Zuijlen | Implicit and explicit higher order time integration schemes for structural dynamics and fluid-structure interaction computations[END_REF]. It is then necessary to consider the smaller of these. Furthermore, depending on the numerical techniques for tracking the free surface, managing the breakup of the boundary, or activating two-phase zones in the fluid domain, many critical time steps must be considered before running a simulation.

It is admitted that it is essential to use numerical tools that are both robust in terms of numerical accuracy and efficient in terms of computation time. Thus, the present study is dedicated to developing alternatives schemes to these methods by presenting numerical schemes that and are not a priori subject to a restrictive stability condition. The latter may have the potential to overcome the previously listed constraints of classical time integration schemes and are able to make the calculations more efficient with a reduced computation time. These chosen alternative methods are based on time pertur-bation approach. We note that several attempts based on perturbation theory and series development of the solution exist in the literature. Perturbation methods have been used to solve nonlinear equations for many years. However, the studies based on time perturbation theory to overcome the critical time step and thus to obtain continuous solutions in time, are not very frequent. Some researchers who have performed some problems using the temporal perturbation are cited. Hamada and Aboanber [START_REF] Hamada | Confirmation of accuracy of generalized power series method for the solution of point kinetics equations with feedback[END_REF], [START_REF] Aboanber | Power series solution (PWS) of nuclear reactor dynamics with newtonian temperature feedback[END_REF] developed the Power Series Solution denoted in their work by "PWS" to obtain approximate transient solutions to point kinetics equations with feedback. Lozi et al. [START_REF] Lozi | A New Reliable Numerical Method for Computing Chaotic Solutions of Dynamical Systems: The Chen Attractor Case[END_REF], [START_REF] Lozi | A new accurate numerical method of approximation of chaotic solutions of dynamical model equations with quadratic nonlinearities[END_REF] described the method of construction of approximate chaotic solutions to dynamical model equations with quadratic nonlinearities, using an accurate numerical method based on time series. Fafard et al. [START_REF] Fafard | Application of an asymptotic method to transient dynamic problems[END_REF], Berrahma-Chekroun et al. [START_REF] Berrahma-Chekroun | Regular Article[END_REF], [START_REF] Berrahma-Chekroun | Resolution of the transient dynamic problem with arbitrary loading using the asymptotic method[END_REF] applied a new approach called the Dynamic Asymptotic method ("DAM") for the solution of transient dynamic problems. Braikat et al. [START_REF] Braikat | A High Order Continuation Based On Time Power Series Expansion And Time Rational Representation For Solving Nonlinear Structural Dynamic Problems[END_REF] performed some tests also on transient nonlinear mechanical problems. Fellah and Fafard [START_REF] Fellah | An approach based on time expansion, Padé approximants and Shanks transformation for linear structural dynamic problems[END_REF] also applied an approach based on time expansion for linear structural dynamic problems.

The perturbation theory used in the last five cited works is related to a developing area of research that deals with nonlinear problems called the Asymptotic Numerical Method and is denoted by "ANM". Thus, the application of ANM to temporal problems was initiated by the last five researchers mentioned previously, and then applied by D. Razafindralandy [START_REF] Razafindralandy | Contribution à l'étude mathématique et numérique de la simulation des grandes échelles[END_REF]. ANM has just been used as a temporal perturbation method in [START_REF] Fafard | Application of an asymptotic method to transient dynamic problems[END_REF]- [START_REF] Razafindralandy | Contribution à l'étude mathématique et numérique de la simulation des grandes échelles[END_REF]. Note that the latter has been extensively studied and successfully applied, as high order perturbation method during the past 30 years to many problems as solid mechanics in [START_REF] Daya | A numerical method for nonlinear eigenvalue problems application to vibrations of viscoelastic structures[END_REF]- [START_REF] Boumediene | Nonlinear forced vibration of damped plates by an asymptotic numerical method[END_REF], as well as fluid mechanics in [START_REF] Cadou | ANM for stationary Navier-Stokes equations and with Petrov-Galerkin formulation[END_REF]- [START_REF] Medale | High performance computations of steady-state bifurcations in 3D incompressible fluid flows by Asymptotic Numerical Method[END_REF]. It is used, in particular, for bifurcation analysis both in the field of structural design and fluid mechanics.

In this work, ANM based on the time perturbation theory is used. Since the objective is to solve evolution problems and develop alternative methods to time integration schemes, it is very reasonable to take the time as the perturbation parameter for ANM. The main idea is to represent the unknowns of the problem given in Eq. [START_REF] Euler | Institutionum Calculi Integralis[END_REF], in the form of a polynomial by using a perturbation parameter denoted by t:

u( t) = u 0 + u 1 t + u 2 t2 + • • • + u N tN + . . . ( 4
)
when N denotes a truncation order. Then, Eq.( 4) is injected into the initial problem given by Eq. [START_REF] Euler | Institutionum Calculi Integralis[END_REF]. When the operator F is a differential polynomial in u, a recurrence relation of the following form is obtained:

u k+1 = F k+1 (u 0 , . . . , u k ), k ∈ N (5)
where F k+1 is another (differential) polynomial operator in u. Therefore, an analytical approximation of the solution is then continuous and differentiable on a temporal validity domain to be determined. The most interesting property of this time perturbation approach is that the computed solutions obtained are continuous in time. However, one must be very careful when manipulating time series, since there exist some series in physics that are divergent [START_REF] Costin | On Borel summation and Stokes phenomena for rank one nonlinear systems of ODE's[END_REF]- [START_REF] Michalik | Summability of formal solutions to the n-dimensional inhomogeneous heat equation[END_REF].

In such a case, the temporal validity domain of ANM becomes extremely small which leads to difficulties to advance the solution and reach the final time. For instance, the time series solutions of the heat equation, the Burgers equation, and the Navier-Stokes equations are generally divergent [START_REF] Costin | On Borel summation and Stokes phenomena for rank one nonlinear systems of ODE's[END_REF]- [START_REF] Michalik | Summability of formal solutions to the n-dimensional inhomogeneous heat equation[END_REF]. This does not imply that we should neglect the divergent series since the latter can lead to very satisfactory results if they are properly handled. For this reason, several methods of summing divergent series exist in the literature. One may cite: Borel summation, Ramanujan's irregular summation, Weierstrass summation, Abel summation, Lindelöf summation, etc... [START_REF]Divergent Series Redux[END_REF]. The most popular one is the Borel-Laplace summation. Towards the end of the 19th century, E. Borel [START_REF] Borel | Mémoire sur les séries divergentes[END_REF] extended the interesting role of divergent series in analysis by calculating the "sum" of a divergent series. It converts the coefficients of a divergent model into a table of approximate values that approach the Borel sum of a divergent series. This method extended from the Cezaro sum is particularly used by mathematicians to compute the sum of the divergent series to study the solutions of differential equations within the theoretical framework. In 1991, it has been proved that the Borel summation (or Borel-Laplace summation) is the most suitable numerically [START_REF] Thomann | Procédés formels et numériques de sommation de séries solutions d'équations différentielles[END_REF]. By using the well-known Padé approximant [START_REF] Brezinski | Rationnal approximation to formal power serie[END_REF], [START_REF] Graves-Morris | Problems and progress in vector Padé approximation[END_REF] combined with the latter, Borel-Padé-Laplace has become the most popular approach for summing divergent series, owing to its algorithmic simplicity.

For this reason, recently Hamdouni and Razafindralandy [START_REF] Razafindralandy | Time integration algorithm based on divergent series resummation, for ordinary and partial differential equations[END_REF] have extended this algorithm numerically on applications to real time-dependent problems in mechanics, limited nevertheless to problems of small size in an exploratory framework. The proposed numerical algorithm based on the time series perturbation ANM, followed by the Borel summation technique of divergent series, is considered as a time integration scheme. It is named Borel-Padé-Laplace and is denoted by "BPL".

Introduction

It has been successfully applied for the resolution of differential equations as well as some partial differential equations [START_REF] Razafindralandy | Contribution à l'étude mathématique et numérique de la simulation des grandes échelles[END_REF], [START_REF] Razafindralandy | Time integration algorithm based on divergent series resummation, for ordinary and partial differential equations[END_REF]- [START_REF] Razafindralandy | Some robust integrators for large time dynamics[END_REF]. In addition, numerical results have shown that the BPL algorithm reduces the number of iterations compared to the classical schemes, and is able to respect the properties of the Hamiltonian systems [START_REF] Razafindralandy | Some robust integrators for large time dynamics[END_REF]. However, many disadvantages of BPL need to be alleviated. One may cite the appearance of some numerical poles that can arise with the use of the Padé approximant, which sometimes block the continuation technique. To avoid these problems, an orientation towards another technique of resummation of divergent series and considered as a time integration scheme has been proposed in [START_REF] Razafindralandy | Considering inverse factorial series as time integration method[END_REF]. It is the Inverse factorial series denoted by "IFS" [START_REF] Razafindralandy | Considering inverse factorial series as time integration method[END_REF]. A first study on the effectiveness of the Inverse factorial series scheme, used as temporal integrator has been initiated just on some academic examples in [START_REF] Deeb | Intégrateurs temporels basés sur la resommation des séries divergentes[END_REF], [START_REF] Razafindralandy | Considering inverse factorial series as time integration method[END_REF].

Finally, Mera et al. [START_REF] Mera | Fast summation of divergent series and resurgent transseries from Meijer-G approximants[END_REF] proposed in 2018 a new method for summing divergent series, called "Meijer-G approximant" and denoted by "MG". The resulting Meijer-G approximations are easily parameterized using an hypergeometric ansatz. They applied the latter to different examples in quantum mechanics and quantum field theory [START_REF] Mera | Fast summation of divergent series and resurgent transseries from Meijer-G approximants[END_REF]. For these examples, results have shown that Meijer-G is more powerful than Borel summation. These approximants are very well suited for problems where the Borel plane has a branch cut as mentioned in [START_REF] Mera | Fast summation of divergent series and resurgent transseries from Meijer-G approximants[END_REF]. The exploitation of the possibility to use this method for timedependent problems will be initiated in this manuscript.

Goals

Since the development of numerical methods for the resolution of nonlinear evolution problems is currently a growing field of research, the main objective of this thesis is to address several needs regarding the development of the recent time perturbation and numerical resummation techniques used as time integration schemes.

• The first goal is to understand more the behavior, develop the feasibility of such integration as well as consolidate the results obtained previously. The properties and qualities of the Inverse factorial series used as a time integration scheme should be studied in detail, since the latter has just been applied as a first prospective study on some simple examples. The potential of the temporal perturbation approaches compared to the classical schemes should be highlighted more on ordinary differential equations, by providing examples where classical schemes suffer from many weaknesses. The stiff or highly non-linear problems should be studied by time perturbation methods.

• The second goal is to advance their use on more complex models represented by PDEs and finally test their relevance and their robustness. The aim is to see how these alternative numerical methods behave for models involving a high number of degrees of freedom. Therefore, we must determine whether they are simultaneously efficient over a large interval time and with an acceptable computation cost. Hence, the efficiency of the temporal methods is not the only objective. The comparison of CPU times between classical and time perturbation methods is an important point to be discussed.

• The third goal of this thesis is to explore the potential of using the recent proposed Meijer-G for summing divergent series, in the context of time perturbation method. The numerical difficulties of the method, as well as its capacity, for transient problems should be presented. The development of certain aspects necessary for the adaptation of MG to the time-dependent problem should be introduced.

equations (ODEs).

• Part III: Application of time perturbation and resummation methods to partial differential equations (PDEs).

In Chapter 1, first, we introduce some necessary theoretical notions on time perturbation and resummation techniques of divergent series. Then, we describe the numerical algorithm of these methods used as time integration schemes, as well as some developments that we have established. Finally, we introduce the recent "Meijer-G approximant" divergent series summation method as an algorithm for solving time-dependent problems, by developing many points necessary for this exploitation. These methods are applied to examples of nonlinear ODEs in Chapter 2 to highlight their advantages of use over classical schemes in some cases. The following examples are handled: Euler equation, simple harmonic oscillator, free and forced Van der Pol oscillators, Elastic hardening spring problem, combustion equation, Lorenz system, nonlinear single pendulum problem,... The interesting properties and difficulties of each time perturbation approach will also be presented.

To be able to apply, for future works, time perturbation methods to simulate Fluid-Structure Interaction (FSI) problems where equation governing the fluid motion (Navier-Stokes equation) are coupled with the equations governing the evolution of the free surfaces (for example level set equation), the new temporal approaches are applied to solve the Navier-Stokes equations in Chapter 3, and the popular level set evolution equations in Chapter 4. The progression difficulties will be done in stages. First of all in Chapter 3, we start to study the efficiency of the methods in terms of feasibility and computation time on the heat equation, the Burgers equation and finally on the incompressible Navier-Stokes equations. We compare the results obtained with those of the literature. For instance, an application to the Lid-Driven cavity problem was chosen in order to discuss the behavior and the efficiency of these solvers, as it is a common benchmark problem in fluid dynamics. In this manuscript, we present the spatial discretization of the numerical resolution of PDEs by the finite element method and some finite difference schemes. Then, in Chapter 4, benchmarks level set problems, are studied by time perturbation and resummation methods. A discussion of a popular technique for level set called the reinitialization procedure, in the context of continuous transient solutions, is presented. Finally, the conclusion will briefly recall the most important results obtained and propose some orientations for future developments.

Chapter 1

TIME SERIES DEVELOPMENT AND RESUMMATION OF DIVERGENT SERIES

This chapter consists of the description of the time perturbation methods and resummation of divergent series that we will use, develop and apply in this manuscript. First, a reminder of some theoretical foundations of the time perturbation and divergent series resummation techniques is presented. Then, the numerical algorithm of these methods used as time integration schemes is also described. ANM with time series development, Borel-Laplace summation, Inverse factorial series, and Meijer-G approximants are presented respectively in subsections 1.1, 1.2, 1.3 and 1.4. It is recalled that as the Meijer-G approximant is used in this manuscript for the first time as a time integration scheme, we will present the necessary adaptation of this transformation in this chapter, before discussing the numerical application. A new version of Meijer-G applied for a series of vectors will be also proposed.

Asymptotic Numerical Method (ANM)

In the framework of this thesis, we are interested to study ANM using the time as its perturbation parameter to solve evolution problems.

Numerical algorithm ANM used as a time integration scheme

Consider the following ordinary differential problem:

M ü + C u + Ku + N l (u) = 0 (1.1)
where u = du dt and ü = d 2 u dt 2 , with the following initials conditions at a given time t 0

   u(t 0 ) = u 0 u(t 0 ) = u 1 (1.2)
where N l (u) denotes a nonlinear term, M , C and K are prescribed scalars. In the case of multiple degrees of freedom, those quantities are matrices and u, u, ü are vectors.

For simplicity, we give the following expression to the nonlinear term N l that is classical in a solid mechanics framework:

N l (u) = u 3 (1.3)
Perturbation: Unknowns of the problem (and time-dependant data) are sought as truncated power series. In this work, time is the perturbation parameter. Global time t is then defined using a given initial time t 0 and a perturbation time denoted by t as

t = t 0 + t (1.4)
Continuous approximation writes as follow:

u( t) = u 0 + N n=1 u n tn (1.5)
with a known u 0 = u(t 0 ). ∀n = {1, . . . , N } , u n are called the series terms and N denotes the truncation order. It is then possible to differentiate this expression as

u( t) = N n=1 nu n tn-1 (1.6) ü( t) = N n=2 n(n -1)u n tn-2 (1.7)
When we replace these expansions into Eq.(1.1), and identify according to the power of t, we notice that it is not easy to find the recurrence formula that allows to calculate the series term. We note that to incorporate polynomial approximations, it is preferable to write all the nonlinear equations to be solved in a quadratic form by introducing new variables [START_REF] Guillot | A Taylor series-based continuation method for solutions of dynamical systems[END_REF]. For this reason, let the following change of variable:

v = u 2 (1.8)
Therefore Eq.(1.3) becomes:

N l (u, v) = uv (1.9)
This yields to the new quadratic equation

M ü + C u + Ku + uv = 0 (1.10)
Writing v( t) as time power series

v( t) = v 0 + N n=1 v n tn (1.11)
The developments of (Eq.(1.5), Eq.(1.6), Eq.(1.7) and Eq. (1.11)) are injected into Eq.(1.10) and by identifying the terms according to the powers of t, a cascade of linear problems is then obtained: t0 :

   2M u 2 + Cu 1 + Ku 0 + u 0 v 0 = 0 v 0 = u 2 0
(1.12) t1 :

   6M u 3 + 2Cu 2 + Ku 1 + u 1 v 0 + u 0 v 1 = 0 v 1 = 2u 1 u 0 (1.13) t2 :    12M u 4 + 3Cu 3 + Ku 2 + u 2 v 0 + u 0 v 2 + u 1 v 1 = 0 v 2 = 2u 2 u 0 + u 2 1 (1.14)
. . . tN-2 :

             M N (N -1)u N + C(N -1)u N -1 + Ku N -2 + N -2 i=0 (N -1 -i)u N -1-i v i = 0 v N -2 = N -2 j=0 u j u N -2-j (1.15)
Series terms u n , ∀n = 1, • • • , N are obtained by the numerical resolutions of those systems. Once these series terms are calculated, the approximate solution Eq.(1.5) is evaluated. The latter is valid, not at discrete time steps, but in a whole interval, corresponding to the numerical radius of convergence of the series.

Radius of convergence R: It represents the limit between convergence and divergence of the series and can be estimated with different techniques. The quality of the solutions improves with the truncation order N for t < R but degrades for t > R. After some practice tests for different examples with t values close to R, it has been shown that it is more interesting to introduce the notion of the domain of validity (or called also a range of validity) for a truncated series.

Range or domain of validity: Solutions of truncated series are valid in the sense of a desired precision, only on a specific range or domain of validity [0, t max ]. By a numerical observation that two successive polynomial approximations move away briskly beyond a certain domain, Cochelin [START_REF] Cochelin | A path-following technique via an asymptotic-numerical method[END_REF] proposed an explicit expression of t max based on the idea that the difference between two solutions evaluated at two consecutive truncation orders is less than a small parameter δ called the tolerance or the precision/accuracy parameter (see Fig. ( u n tn (series truncated at order N ), and u t,N -1 = N -1 n=0 u n tn (series truncated at order N -1).

Approximating the denominator by t||u 1 || 2 , we get the explicit expression of t max proposed in [START_REF] Cochelin | A path-following technique via an asymptotic-numerical method[END_REF]:

t max = δ ||u 1 || 2 ||u N || 2 1 N -1 (1.17)
where δ denotes a given tolerance and . 2 is the Euclidean norm. Eq.(1.17) depends on two parameters, the truncation order N and the accuracy parameter or tolerance δ. The order N is usually chosen a priori to minimize the computation time. This criterion requires only an insignificant computation time for the algorithm resolution.

Note that another way to evaluate the domain of validity exists in [START_REF] Cochelin | A path-following technique via an asymptotic-numerical method[END_REF], [START_REF] Baguet | Stability of thin-shell structures and imperfection sensitivity analysis with the Asymptotic Numerical Method[END_REF]. As solutions are continuous in time, they can be evaluated at any time t ∈ [0, t max ].

As can be seen in Fig. (1.2), we note that for this example the more terms added in the representation, the wider is the range of validity. u( t) is not able to lead the exact solution for different truncation orders N , it may not be possible to reach u( t) for a large value of t until a final time t. Therefore, a continuation procedure is proposed to reach a specific given time t.

Continuation:

To reach high values of t, a continuation procedure is performed, that is, the problem is solved again with the last valid solution. The continuation procedure is easy, it is sufficient to determine a new starting point u 0 after using the perturbation technique, and then re-apply the latter technique from this point. The new starting point, denoted by u j+1 0 on the piece j of the solution, is calculated using its asymptotic series for the value of t max at each calculation step ( see Fig. (1.3)). It can be summarized as follows:

Algorithm 1 Continuation technique for ANM with time as perturbation parameter

Inputs:

1. The initial conditions u 0 and u 1 .

2. The truncation order N . 3. The tolerance or the accuracy parameter δ. 4. The final time t end or the domain in which we want to obtain the curve, for example: [0, t end ].

Step 1: Compute series terms until a given order N .

Step 2: Evaluate the upper bound of the validity domain t max using Eq.(1.17).

Step 3: Evaluate a new starting point u 0 at a new time t 0 = t 0 + t max . The new point u 0 = u(t 0 + t max ) is denoted by u j+1 0 on the piece j of the solution at time t j+1 .

Step 4: Go to Step 1 and repeat the algorithm until reaching the final time t end .

The number of times to a continuation iteration will be called the number of continuation steps. This time series development allows obtaining a numerical approximation of the solutions continuous in time, which is different from a discrete approach. The performance of the algorithm is therefore based on the domain of validity. If the latter is too small, a high number of continuous steps becomes necessary.

Time series and Padé approximants

It is possible to enlarge the domain of validity of the series: [0, t max ] by replacing the polynomial with a rational fraction (Padé approximant in this work) which is asymptotically equivalent to it. This procedure has proven to be efficient in many cases [START_REF] Najah | A critical review of asymptotic numerical methods[END_REF]- [START_REF] Mallil | An iterative process based on homotopy and perturbation techniques[END_REF]. The polynomial approximation can be replaced by a rational one, using Padé approximants.

The latter are very frequently used in computer sciences and are well known for their efficiency even to extend a series out of its convergence disk. This new representation is written in the following form:

u P adé,N ( t) = u 0 + N -1 i=1 R N -1-i ( t) Q N -1 ( t) ti u i (1.18)
where R i and Q i are polynomials of degree i. For more details about the popular Padé, readers are refered to [START_REF] Brezinski | Rationnal approximation to formal power serie[END_REF], [START_REF] Graves-Morris | Problems and progress in vector Padé approximation[END_REF].

According to several test calculations carried out, rational approximations for different truncation orders give very close results when the parameter t is within a certain domain, but beyond this domain, they diverge strongly. Thus, the simple way to define the range of validity of this rational representation is to extend the criterion used in Eq.(1.17). In this case, it is sufficient to require that the divergence between two rational solutions at two consecutive orders remains small at the end of the step. Therefore, the range of validity of the rational representation is sought by dichotomy via a tolerance criterion δ P , as in [START_REF] Elhage-Hussein | A numerical continuation method based on Padé approximants[END_REF]: With Eq. (1.19), t max of Padé approximants is computed and a new starting point for the asymptotic expansions Eq.(1.5) can be defined by introducing the value of t max in the expression Eq. (1.18). Thus, a continuation method is defined. Remark 1.1.2 If one needs to approximate numerically the derivatives of the Padé solution (i.e the derivative of Eq. (1.18)), one can proceed as follows: First, consider for example a truncation order N = 5, therefore we have:

δ P = u P adé,N ( t) -u P adé,N -1 ( t)
u(t) = N n=0 u n t n = u 0 + u 1 t + u 2 t 2 + u 3 t 3 + u 4 t 4 + u 5 t 5 (1.20)
Therefore the derivative of u(t) is

u(t) = N n=1 nu n t n-1 = u 1 + 2u 2 t + 3u 3 t 2 + 4u 4 t 3 + 5u 5 t 4 (1.21)
On the other hand, u(t) can be written in the following form:

u(t) = u0 + u1 t + u2 t 2 + u3 t 3 + u4 t 4 (1.22) where u0 = u 1 , u1 = 2u 2 , u2 = 3u 3 , u3 = 4u 4 u4 = 5u 5 (1.23)
To calculate the derivatives of the Padé solution (i.e the derivative of Eq. (1.18)), one can evaluate the Padé solution of Eq. (1.22) with the new series terms u0 , u1 , u2 , u3 , u4 . In the same way, one can evaluate the second derivative of the Padé solution of u(t) by writing ü(t) in the same way as Eq. (1.22).

Remark 1.1. [START_REF] Courant | Über die partiellen Differenzengleichungen der mathematischen Physik[END_REF] We will not present the ANM with the time series followed by the Padé approximants in all of the examples in this document because we have found numerical instability in the majority of ODEs cases when using the ANM+Padé approximant for some temporal large problem. Note that the concept of Padé approximant has been defined here for its later use in Borel-Padé Laplace subsection 1.2.2.

Divergent series

The drawback of the ANM approach using the time as its perturbation parameter is that the time series might be divergent. In such a case, the range of validity of this representation in the temporal domain becomes extremely small. Therefore, classical continuation is no longer efficient. An example of divergent time series will be presented [START_REF] Razafindralandy | Time integration algorithm based on divergent series resummation, for ordinary and partial differential equations[END_REF], [START_REF] Deeb | Intégrateurs temporels basés sur la resommation des séries divergentes[END_REF]. The behavior of ANM for different truncation orders N , without any continuation procedure is presented in Fig. (1.4). It is observed that when the truncation order N increases, the computed solutions tend to diverge earlier, contrary to what has been observed in Fig. (1.2).

Example 1.1.4 Consider the Euler equation

t 2 u + u = t (1.24)
where u = du dt , with the following initial condition

u(t 0 = 0) = 0 (1.25)
The formal power series solution of this equation is

u(t) = k≥0 (-1) k k!t k+1 (1.26)
This series is divergent for any t > 0. ANM would not be able to provide an analytical solution. In the absence of theoretical information on the type of time series, we cannot be certain whether the infinite numerical series is convergent or divergent because the solution u(t) cannot be obtained numerically when N tends towards infinity for any computer memory.

How can we study divergent series ?

As there are physical equations that lead to divergent asymptotic series, these cases should not be ignored. In the study of equations, whether they are ordinary differential equations (ODEs) or partial differential equations (PDEs), the aim is to recover information as much as possible regarding the solution, by summing up the divergent asymptotic series.

To obtain an accurate solution, one must be interested in the summation methods of divergent series that give the highest information about the associated function, even if the latter is truncated. Therefore, ANM is unusable when we have a very small or zero radius of convergence, hence the need to apply the summation techniques of divergent series. As mentioned in the introduction, one of the most popular summation methods is the Borel-Laplace summation that will be exposed in the next section.

Borel-Laplace (BL)

First, a theoretical necessary reminder of the summation of Borel-Laplace (or simply Borel summation) denoted by BL is presented. The reader is referred to [START_REF] Michalik | Summability of formal solutions to the n-dimensional inhomogeneous heat equation[END_REF], [START_REF] Borel | Mémoire sur les séries divergentes[END_REF], [START_REF] Malgrange | Fonctions multisommables[END_REF], [START_REF] Ramis | Séries divergentes et théories asymptotiques[END_REF] for more details concerning the definition of Borel summation method. Then the numerical algorithm of Borel-Padé-Laplace (BPL) [START_REF] Razafindralandy | Time integration algorithm based on divergent series resummation, for ordinary and partial differential equations[END_REF] is exposed. We start with some useful necessary definitions and theorems.

Theoretical reminder

First, we note that Watson [START_REF] Watson | The transformation of an asymptotic series into a convergent series of inverse factorials [Memoir crowned by the Danish Royal Academy of Science[END_REF] classified all divergent series according to a so-called Gevrey index, after noticing that the general term of the most frequently identified divergent series always behaved as CA n (n!) 1/k for certain constants C, A, and k > 0 when n is large enough, which yields to the following definition:

Definition 1.2.1 A formal power series u(t) = ∞ n=0
u n t n is a Gevrey series of index or order k, if there exists two constants C > 0 and A > 0 such that

|u n | ≤ CA n (n!) 1/k , ∀n ∈ N (1.27)
In particular, the Gevrey series of order 0 is the convergent series. This definition will be an interesting condition for the Borel transform, that we will introduce in the following, in order to have a non-zero radius of convergence. Note that, it is well-known that most of the series arising in engineering problems are k-Gevrey series for some positive rational number k. In what follows, we will present the Borel-Laplace resummation only for the case of k = 1.

Note that Maillet [START_REF] Maillet | Sur les séries divergentes et les équations différentielles[END_REF] has proven that formal series that are solutions of an algebraic ODE are Gevrey series. Malgrange [START_REF] Malgrange | Sur le théorème de Maillet[END_REF] has extended this observation to the analytic nonlinear ODEs. The theory for PDEs is not quite complete. Some equations in fluid mechanics have been investigated, such as the heat equation and its generalizations and the Burgers equation [START_REF] Lysik | Borel summable solutions of the Burgers equation[END_REF], [START_REF] Michalik | Summability of formal solutions to the n-dimensional inhomogeneous heat equation[END_REF], [START_REF] Łysik | Formal solutions of semilinear heat equations[END_REF]- [START_REF] Balser | Divergent solutions of the heat equation: On an article of Lutz, Miyake and Schäfke[END_REF]. For the Navier-Stokes equations, summability conditions to the formal power series are provided recently by Costin et al. [START_REF] Costin | Divergent expansion, Borel summability and three-dimensional Navier-Stokes equation[END_REF].

Remark 1.2.2 Note that in the case of a formal series solution of a linear equation with polynomial coefficients, the Gevrey index k can be determined from the Newton polygon of the equation. We will not present the details for this remark since most of the cases that will be considered in the following chapters are nonlinear. The interested reader is invited to consult the following references [START_REF] Miyake | Newton polygons and formal Gevrey indices in the Cauchy-Goursat-Fuchs type equations[END_REF], [START_REF] Hashimoto | Newton polygons and Gevrey indices for linear partial differential operators[END_REF].

We recall that we will present the definitions required for Borel-Laplace resummation only for the case of the Gevrey series of order k = 1.

The approximate obtained solution will be an analytical function in a sector and has the series u(t) = ∞ n=0 u n t n as the Gevrey asymptotic, as defined below: Definition 1.2.3 An analytic function u(t) in an open sector S, is Gevrey-asymptotic of order I to a series u(t) = n≥0 u n t n at the origin in S, if for any compact subsector T ⊂ S, there exists two constants C and A such that

∀n ∈ N, ∀t ∈ T, u(t) - N -1 k=0 u k t k ≤ CA n |t| n (n!) I (1.28) Definition 1.2.4 The Borel transform of the time power series u(t) = ∞ n=0 u n t n is Bu(ξ) = ∞ n=0 u n+1 n! ξ n (1.29)
The Borel transform Bu(ξ) has a non-zero convergence radius if u is a Gevrey series analytical in the neighborhood of 0.

To incorporate the notion of Laplace transformation, one must first establish the concept of growth at most exponential to infinity. Definition 1.2.5 Let P (ξ) be an analytic function in a domain containing a semi-line starting from the origin denoted by d. We say that P (ξ) has an exponential decay at the infinity in the direction d if there exists two constants A and C such that

|P (ξ)| ≤ Ae C|ξ| (1.30)
as |ξ| tends to infinity along d.

We can now establish the definition of the Laplace transform.

S d u(t) = u 0 + LP (t) = u 0 + d P (ξ)e -ξ/t dξ. (1.32) Remark 1.2.8
The Laplace transform can be considered as the inverse of the Borel transform when it is applied formally to a power series.

After the required definitions and theorems, we can establish the following three steps of Borel-Laplace to calculate the Borel sum S d u(t) in a direction d:

(i) Compute the Borel transform Bu(ξ).

(ii) Prolongate naturally Bu(ξ) into an analytical function P (ξ).

(iii) Compute the Laplace transform of P along a direction d which permits to go back to the originally space.

These three steps constitute the Borel-Laplace summation method. It is simple to understand: the coefficients of a series with a zero radius of convergence typically grow factorially to high orders. This factorial growth is eliminated in the first step, and the transformed Borel series is easy to handle because it has a non-zero finite radius of convergence. By summing the Borel-transformed series one finds the function P (ξ), and then the Laplace transform can be computed to find finally the Borel sum. Remark 1.2.9 When d = R + , we will note the Borel sum S d u(t) by Su(t).

Remark 1.2.10

The generalization of all the above definitions for Gevrey series of index k is possible. For more information about the resummation of Borel-Laplace for k-Gevrey series, the reader is invited to consult the following references [START_REF] Malgrange | Fonctions multisommables[END_REF], [START_REF] Ramis | Séries divergentes et théories asymptotiques[END_REF].

Numerical algorithm BPL used as a time integration scheme

The previous three steps are translated into the following numerical algorithm BPL. It is organized around the following three stages: Borel transform (B) + Approximation of Padé (P ) + Laplace transform (L) (BPL). We present this numerical algorithm for Gevrey series of index k = 1 in a direction d = R + . Since the singular directions of the Borel transform are not known a priori, therefore this direction is set as the real positive axis in this work as in [START_REF] Razafindralandy | Time integration algorithm based on divergent series resummation, for ordinary and partial differential equations[END_REF].

• Fist step: Borel transform

Numerically, for the first stage, series terms are computed up to an order N

u( t) = N n=0 u n tn (1. 33
)
where t is the perturbation time. This yields that numerically Borel transform B(ξ) is of order N -1. It is given by:

Bu(ξ) = N -1 n=0 u n+1 n! ξ n (1.34)
• Second step: Prolongation via Padé approximant For the second stage, an extension of the Borel transform is obtained by establishing an approximation of Padé:

P (ξ) = A 0 + A 1 ξ + A 2 ξ 2 + • • • + A L ξ L 1 + B 1 ξ + • • • + B M ξ M (1.35)
where L and M are two integers such that L + M = N -1 (because the Borel transform is of order N -1), A 0 , A 1 , . . . , A L and B 1 , . . . , B M are the coefficients of Padé approximants. In the cases of large size problems (for example when solving a PDE), it is recommended to use a vector version of Padé [START_REF] Xu | Matrix Padé approximation: definitions and properties[END_REF]. Note that prolongation of Borel transform can be done with different methods other than Padé.

The interested reader is referred to [START_REF] Thomann | Procédés formels et numériques de sommation de séries solutions d'équations différentielles[END_REF] for furthers details, but these methods are designed for linear equations.

• Third step: Laplace transform For the last stage, the Laplace transform is applied to the Padé approximants to return to the temporal space, with a Gauss-Laguerre quadrature [START_REF] Kzaz | Convergence acceleration of the Gauss-Laguerre quadrature formula[END_REF]. Note that the Gauss-Laguerre method has been chosen for the evaluation of the Laplace transform because it is one of the most well-known quadrature methods applied to semi-infinite integrals. But definitely, other methods can be considered [START_REF] Deun | On computing rational Gauss-Chebyshev quadrature formulas[END_REF]- [START_REF] Kumar | Quadrature Formulas for Semi-Infinite Integrals[END_REF]. The Borel sum is

Su( t) = u 0 + ∞ 0 P (ξ)e -ξ/ tdξ (1.36)
then applying Gauss-Laguerre quadrature to approximate the integral contained in the Borel sum (the Laplace transform). Before elaborating how to process the Borel sum numerically, a brief review of Gauss-Laguerre quadrature is presented.

Recall 1.2.11

In numerical analysis, Gauss-Laguerre quadrature is an extension of the Gaussian quadrature method for approximating the value of integrals of the following kind:

∞ 0 f (x)e -x dx (1.37)
by the following sum

N G i=1 w i f (x i ) (1.38)
where x i points (or nodes) are real, distinct, unique and are the roots of N G Laguerre polynomials denoted by L n (x) (N G number of Gauss points), and the weight w i is given by

w i = x i (n + 1) 2 [L n+1 (x 2 i )] (1.39)
Using the following relation:

P (ξ)e -ξ/ t = P (ξ)e -ξ/ te ξ e -ξ (1.40)
then applying the Gauss-Laguerre quadrature to Eq.(1.40) yield to

∞ 0 P (ξ)e -ξ/ tdξ ≈ N G i=1 P (ξ i )e -ξ(
Su( t) = u 0 + t N G i=1 P (ξ i t)w i (1.43)
Note that changing the order of this scheme is very simple. We just have to replace the value of N with another one.

This method can be considered as a method to extend the range of validity of the series when the latter is convergent and plays its task as a resummation method when the series diverges [START_REF] Razafindralandy | Time integration algorithm based on divergent series resummation, for ordinary and partial differential equations[END_REF].

Since the initial series Eq.(1.33) is truncated up to an order N , and because of the numerical effects, numerical resummation methods are not able to lead the exact solution, i.e by applying the resummation only once, it may not be possible to reach Su N ( t) for a large value of t until a final time t. Therefore, Borel sum Eq.(1.43) allows approximating analytically the exact solution only in a certain time interval. Moreover, the definition of the range of validity [0, t max ] of BPL is still not evaluated explicitly like that of ANM (see Eq.(1.17)). Therefore, a continuation procedure is proposed to reach specific given time t. The continuation of BPL algorithm might be carried out using a criterion based on the residual since it is the most useful criterion in engineering science. The residual Res(t) of the differential problem (Eq.(1.1)) is defined by

Res(t) = M S(t) + C Ṡ(t) + KS(t) + N l (S(t)) (1.44)
where -S(t) = Su(t) is the Borel sum presented in Eq.(1.43).

-Ṡ(t) = dSu(t) dt is the first derivative of Borel sum.

-S(t) = d 2 Su(t) dt 2 is the second derivative of Borel sum. replaced in the initial equation of the problem (Eq.(1.1)), M, C and K are defined before. Those latter should be determined in order to evaluate the residual. Here it should be mentioned that to calculate the derivatives there are (as for the Padé approximant) at least two ways of doing it:

1. Calculate the derivatives of the series and establish an identification between the terms. (see remark 1.1.2). 2. Establish an analytical expression of the derivatives of the "Borel expression" which is an original contribution of this thesis. Those latters are established in Appendix A.

We will evaluate numerically some examples of the two ways of calculating the derivatives for high simulation times in the next chapter.

The analytical approximations are examined as long as

Res(t) 2 < (1.45)
where is a small positive parameter called the residual tolerance or the accuracy parameter and . Step 2: Apply the three steps of BPL (Borel transform, prolongation via Padé approximant, Laplace transform) to calculate the Borel sum Su(t).

Step 3: Apply the first and second derivatives approximations of BPL to calculate Ṡ(t) and S(t). (See Appendix A).

Step 4: Evaluate the residual Res(t) given in Eq.(1.45).

Step 5: Evaluate the validity domain t max such that t max is the last value of t for which Res(t) 2 < .

Step 6: Evaluate a new starting point u 0 at a new time t 0 = t 0 + t max .

Step 7: Go to Step 1 and repeat the algorithm until reaching the final time t end .

This approach with the residual criterion is quite a time expensive, but we will see in the numerical applications if it is fast enough to outperform other numerical schemes. An evaluation of the range of validity of BPL (t max of BPL) without using the residual computation is preferable but is not yet available. This type of solver is promising to make the numerical simulations more efficient in terms of capture of the fast dynamic phenomenon and numerical stability over a large number of cycles and computation time. This algorithm has been successfully applied on academic examples and some examples of Hamiltonian systems in [START_REF] Razafindralandy | Contribution à l'étude mathématique et numérique de la simulation des grandes échelles[END_REF], [START_REF] Razafindralandy | Time integration algorithm based on divergent series resummation, for ordinary and partial differential equations[END_REF]- [START_REF] Razafindralandy | Some robust integrators for large time dynamics[END_REF].

Strengths and weaknesses of BPL

Numerical results in [START_REF] Razafindralandy | Time integration algorithm based on divergent series resummation, for ordinary and partial differential equations[END_REF], [START_REF] Deeb | Intégrateurs temporels basés sur la resommation des séries divergentes[END_REF], [START_REF] Deeb | Comparison between Borel-Padé summation and factorial series, as time integration methods[END_REF] have shown that the BPL algorithm presents some interesting properties. To illustrate the BPL algorithm, we use one example that have already been discussed in [START_REF] Razafindralandy | Time integration algorithm based on divergent series resummation, for ordinary and partial differential equations[END_REF], [START_REF] Deeb | Intégrateurs temporels basés sur la resommation des séries divergentes[END_REF], [START_REF] Deeb | Comparison between Borel-Padé summation and factorial series, as time integration methods[END_REF].

Example 1.2.12 Consider again the Euler equation

t 2 u + u = t (1.46)
where u = du dt , with the following initial condition

u(t 0 = 0) = 0 (1.47)
The formal power series solution of this equation is

u(t) = n≥0 (-1) n n!t n+1 (1.48)
which is a divergent series. The Borel transform of Eq. (1.48) is

Bu(ξ) = ∞ n=1 (-1) n ξ n (1.49)
This series can be naturally extended into the analytic function

P (ξ) = 1 1 + ξ (1.50)
Therefore, the Laplace transform in the direction d = R + is

LP (t) = ∞ 0 P (ξ)e -ξ/t dξ = ∞ 0 e -ξ/t 1 + ξ dξ (1.51)
this yield to the Borel sum Many advantages can be obtained by using BPL [START_REF] Razafindralandy | Contribution à l'étude mathématique et numérique de la simulation des grandes échelles[END_REF], [START_REF] Razafindralandy | Time integration algorithm based on divergent series resummation, for ordinary and partial differential equations[END_REF]- [START_REF] Razafindralandy | Some robust integrators for large time dynamics[END_REF]. One may cite: (i) Ability to solve the problem over a large time interval.

Su(t) = u 0 + ∞ 0 e -ξ/t 1 + ξ dξ = ∞ 0 e -ξ/t 1 + ξ dξ. ( 1
(ii) BPL holds the properties of Hamiltonian systems (symplecticity, iso spectrality).

(iii) BPL needs less steps than the classical scheme to reach the same final time.

Padé approximants constitute simultaneously the strength as well as the weakness of the BPL algorithm presented. The latter makes it possible to identify the possible poles of the function to be approximated, which thus enables to choose the direction of the integration of the Laplace transformation. But when the Padé table is not normal, the linear problem resulting from the Padé algorithm becomes non-invertible. And if it is not the case, the linear system is ill-conditioned in general for many examples, which can lead to the appearance of spurious null pole pairs [START_REF] Stahl | Spurious poles in Padé approximation[END_REF]. To avoid the problems induced by the Padé approximants, several techniques exist in the literature. For example, one can cite the use of the "Froissart doublet" [START_REF] Stahl | Spurious poles in Padé approximation[END_REF]. One can also mention the use of the Padé approximant where the denominator (and the poles) are prescribed. Another technique that can be also listed is the use of conformal mapping to extend the Borel transform series, but the latter requires the knowledge of the nearest pole of Bu.

In [START_REF] Razafindralandy | Considering inverse factorial series as time integration method[END_REF], [START_REF] Deeb | Comparison between Borel-Padé summation and factorial series, as time integration methods[END_REF], the weaknesses of the BPL algorithm for the numerical computation of the Borel sum have been discussed through different examples introducing some cases where:

)).

(iii) There is some discretization errors due to the numerical integration of Gauss-Laguerre. To avoid these problems by Padé approximants, another numerical resummation method namely Inverse factorial series and denoted by "IFS" is proposed as a time integration scheme in [START_REF] Razafindralandy | Considering inverse factorial series as time integration method[END_REF] and is presented in the following subsection.

Inverse factorial series (IFS)

The Borel sum can be represented by a generalized factorial series which is computed directly from the terms of the series, without passing through the Borel space. The key point of this approach is that it avoids the three-step calculation and there is no need to evaluate P (ξ). The introduction of the inverse factorial series as an asymptotic of a power series was initiated at the beginning of the last century: [START_REF] Watson | The transformation of an asymptotic series into a convergent series of inverse factorials [Memoir crowned by the Danish Royal Academy of Science[END_REF], [START_REF] Nielsen | Sur les séries factorielles et la fonction gamma (extrait d'une lettre adressée à M. N. de Sonin, à Saint-Pétersbourg)[END_REF], [START_REF] Nielsen | Les séries de factorielles et les opérations fondamentales[END_REF]. More recently, the effective summation by the inverse factorial series has been numerically investigated in [START_REF] Thomann | Procédés formels et numériques de sommation de séries solutions d'équations différentielles[END_REF], [START_REF] Delabaere | Sommation effective d'une somme de Borel par séries de factorielles[END_REF], [START_REF] Weniger | Summation of divergent power series by means of factorial series[END_REF], and very recently used as a time integration scheme in [START_REF] Razafindralandy | Considering inverse factorial series as time integration method[END_REF], [START_REF] Deeb | Comparison between Borel-Padé summation and factorial series, as time integration methods[END_REF] for which the details will be presented in the following.

First, we introduce some useful definitions and theoretical results. Then, we present the numerical algorithm of IFS used as a time integration scheme. Note that the link between Borel sum and Inverse factorial series is shown in [START_REF] Razafindralandy | Considering inverse factorial series as time integration method[END_REF] and will be recalled here since it is important to illustrate it. First, we consider again equation Eq. (1.32):

Su(t) -u 0 = ∞ 0 P (ξ)e -ξ/t dξ (1.53) Let ξ = -ln µ, therefore ∞ 0 P (ξ)e -ξ/t dξ = 1 0 P (-ln µ)µ 1 t -1 dµ = ∞ n≥0 b n 1 0 (1 -µ) n µ 1 t -1 dµ (1.54)
where b n are the coefficients of the Taylor expansion of P (-ln µ) at µ = 1 given by:

b n = 1 n! n+1 k=1 |S(n, k -1)|u k (1.55)
where S(n, k) denotes the Stirling numbers of the first kind. Therefore, Eq.(1.53) becomes:

Su(t) -u 0 = ∞ n≥0 b n 1 0 (1 -µ) n µ 1 t -1 dµ (1.56)
After integration by part, we get

Su(t) -u 0 = ∞ n≥0 b n 1 0 (1 -µ) n µ 1 t -1 dµ (1.57) = ∞ n≥0 b n n! 1 t ( 1 t + 1) . . . ( 1 t + n) (1.58) = ∞ n≥0 b n n!t n+1 (1 + t) . . . (1 + nt) (1.59)
Finally, we have

Su(t) = u 0 + ∞ n≥0 b n n!t n+1 (1 + t) . . . (1 + nt) (1.
I(t) = u 0 + ∞ n=0 b n n!t n+1 (1 + t) . . . (1 + nt) (1.61) with b n = 1 n! n+1 k=1 |S s (n, k -1)|u k (1.62)
where S s (n, k) denotes the Stirling numbers of the first kind.

I(t) is unique and absolutely convergent for sufficiently high values of the real part of (1/t) [START_REF] Nielsen | Les séries de factorielles et les opérations fondamentales[END_REF], [START_REF] Thomann | Resommation des series formelles. Solutions d'équations différentielles linéaires ordinaires du second ordre dans le champ complexe au voisinage de singularités irrégulières[END_REF], [START_REF] Nielsen | Recherches sur les séries de factorielles[END_REF] and converges towards S û(t). Note that Eq.(1.61) is valid for the 1-summable series in a direction d = R + .

Numerical algorithm IFS used as a time integration scheme

Numerically, I(t) is truncated up to an order N and is denoted by I N ( t) with t the perturbation time: Therefore, the Inverse factorial series summation

I N ( t) = u 0 + N n=0 b n n! tn+1 (1 + t) . . . (1 + n t) (1.
I( t) of S û( t) satisfies |Su N ( t) -I N ( t)| ≤ C N tN+1 N ! (1 + t) • • • (1 + N t)(Re 1 t -1 T ) (1.64)
where Su N ( t) is the Borel sum truncated at an order N , and

C N = CT 1 T (N + 1 + 1 T ) N +1+ 1 T (N + 1) N (1.65)
Remark 1.3.3 When the truncation order N is large, Eq. (1.62) can not be used because of the numerical explosion of the Stirling numbers S(n, k) [START_REF] Deeb | Intégrateurs temporels basés sur la resommation des séries divergentes[END_REF].

Hence, an equivalent method to evaluate I N -1 ( t), proposed in [START_REF] Thomann | Resommation des series formelles. Solutions d'équations différentielles linéaires ordinaires du second ordre dans le champ complexe au voisinage de singularités irrégulières[END_REF], is described:

Suppose that u( t) is 1-summable in a direction d, with an angle θ with the positive half-axis.

Let τ l = s l e iθ with l = (1, 2, . . . , N ), where {s 1 , s 2 , . . . , s N } is any sequence of complex numbers. Let z = 1 t , y = ze iθ (1.66)

a 1 = u 0 z , a 2 = u 1 z 2 , . . . , a m = u m-1 z m , m = 1, . . . , N (1.67)
Note that if the series are summable in the R + direction i.e d = R + , we have y = z and we will choose s l = l for l ≥ 1.

We call the Inverse factorial series the following

I N -1 ( t) = 1 t N -1 n=0 v n+1 = z N -1 n=0 v n+1 = z(v 1 + v 2 + • • • + v N ) (1.68)
where v n is the nth term of IFS which are calculated using the following recursive algorithm based on

v (j) n+1 = τ n-1 v (j) n + yv (j+1) n y + τ n n ≥ 1, j ≥ 1, (1.69)
with

v (1) 1 = a 1 v (2) 1 = a 2 . . . v (N ) 1 = a N (1.70)
The terms v

(1)

1 , v (1)
2 , . . . , v (1) n , v

(1)
n+1 are the terms of the IFS denoted by: v 1 , v 2 , . . . , v n , v n+1 .

Remark 1.3.4 Note that the evaluation of the IFS sum with this approach should start from t 0 = 0 (because of Eq. (1.66)). For this reason, for all the examples that we will present in this manuscript, we will fix numerically the initial time t 0 to a very small value close to 0.

The IFS algorithm can overcome the problems encountered by BPL, therefore it is a good alternative to BPL. Note that only one step of this algorithm has been considered in [START_REF] Deeb | Comparison between Borel-Padé summation and factorial series, as time integration methods[END_REF] i.e. no continuation has been applied, and only one example with continuation procedure in [START_REF] Deeb | Intégrateurs temporels basés sur la resommation des séries divergentes[END_REF]. Below is an example extracted from [START_REF] Deeb | Intégrateurs temporels basés sur la resommation des séries divergentes[END_REF] where there are poles in Padé approximant using BPL, showing that IFS (denoted by Factorial in [START_REF] Deeb | Intégrateurs temporels basés sur la resommation des séries divergentes[END_REF] where there is a pole in Padé approximant.

We propose a continuation of IFS algorithm using the criterion of the residual. First, note that relation Eq.(1.64) can be used to evaluate the error useful to apply a continuation as mentioned in [START_REF] Deeb | Comparison between Borel-Padé summation and factorial series, as time integration methods[END_REF]. The residual Res(t) is defined by:

Res(t) = M Ï(t) + C İ(t) + KI(t) + N l (I(t)) (1.71)
where İ(t) the first derivative of IFS and Ï(t) is the second derivative of IFS. Therefore the approximations of derivatives of the IFS algorithm are required for the residual evaluation replaced in the initial equation of the problem (Eq.(1.1)). Here it should be mentioned that to calculate the derivatives there are (as for the Padé approximant) at least two ways of doing it: 1. Calculate the derivatives of the series and establish an identification between the terms. (see remark 1.1.2). 2. Establish an analytical expression of the derivatives of the IFS sum which is an original contribution of this thesis. Those latter are established in Appendix B. We will evaluate numerically some examples of the two ways of calculating the derivatives of IFS for high simulation times in the next chapter. Therefore, the analytical approximations are examined as long as

Res(t) 2 < (1.72)
where is a small positive parameter called the tolerance and .

Meijer-G approximant (MG)

According to the recent research [START_REF] Mera | Fast summation of divergent series and resurgent transseries from Meijer-G approximants[END_REF], the more powerful method of summing divergent series than Borel-Laplace is the new Meijer-G summation with Hypergeometric approximants denoted by "MG". It is a simple algorithm proposed in 2018 by Mera et al. [START_REF] Mera | Fast summation of divergent series and resurgent transseries from Meijer-G approximants[END_REF] that allows a fast and accurate low-order Borel summation. It replaces the Padé approximants in the Borel plane with more general and flexible hypergeometric functions that are, in principle, of an arbitrarily high order. MG is a generalized Borel-Hypergeometric approach where the hypergeometric ansatz is used to transform the coefficients of a divergent series into an array of hypergeometric approximations in the Borel plane [START_REF] Mera | Fast summation of divergent series and resurgent transseries from Meijer-G approximants[END_REF]. This technique has been successfully applied to various examples in quantum mechanics and quantum field theory (QFT) (for example: zero-dimensional φ 4 theory, the quartic anharmonic oscillator, and the three-dimensional self-avoiding walk case in the N vector model,...) used as reference models for summing divergent series in these fields [START_REF] Mera | Fast summation of divergent series and resurgent transseries from Meijer-G approximants[END_REF], [START_REF] Antipin | Resummation in QFT with Meijer G-functions[END_REF]. It has also been applied for test cases where the perturbed expansion is not Borel summable along the positive real axis due to some obstacles (for example: summation of transseries/multi-instanton expansions). In all of these cases, the Meijer-G summation approach works effectively, provides a fast way to calculate the sum of divergent series, and outperforms Borel-Padé-Laplace technique [START_REF] Mera | Fast summation of divergent series and resurgent transseries from Meijer-G approximants[END_REF].

In this manuscript, we will transform MG resummation to a time perturbation approach. According to the literature [START_REF] Mera | Fast summation of divergent series and resurgent transseries from Meijer-G approximants[END_REF], this resummation method is particularly suitable for series with a branch cut. However, the encouraging results presented in [START_REF] Mera | Fast summation of divergent series and resurgent transseries from Meijer-G approximants[END_REF] motivated us to adapt this technique to sum divergent time series and investigate it further. Since MG is a recent summation method, we will not limit the analysis in this manuscript only for BPL and IFS summations, MG will be also considered as a new time integration scheme. First of all, we review the Meijer-G approximant algorithm as it was proposed initially in [START_REF] Mera | Fast summation of divergent series and resurgent transseries from Meijer-G approximants[END_REF]. Then, we will transform it into a time integration scheme by adding the necessary explanations on the subtle points of this transformation. Finally, we propose a new method for applying Meijer-G to a vector series, which has not been addressed in any of the previously listed studies of MG resummation.

Numerical algorithm of Meijer-G (MG)

The Meijer-G algorithm is described here [START_REF] Mera | Fast summation of divergent series and resurgent transseries from Meijer-G approximants[END_REF]. First, we assume in the description of this algorithm that the truncation order N is an odd number (because we will use (N -1)/2 in the following). The generalization for even orders will be presented at the end of this subsection. Let

u(t) = N n=0 u n t n (1.73)
with normalized coefficients, u 0 = 1. This algorithm consists of the following four steps:

• Fist step: Borel transform The first stage is to compute Borel-transformed coefficients denoted by b n and given by

b n = u n n! (1.74)
and the ratios of the consecutive Borel-transformed coefficients denoted by r N (n) and given by:

r N (n) = b n+1 b n (1.75)
• Second step: Hypergeometric ansatz The second stage is to make the ansatz that r N (n) is a rational function of n defined as:

r N (n) = l m=0 p m n m 1 + l m=1 q m n m (1.76)
where n is from 0 to N -1, l = (N -1)/2 and for this reason, the truncation order N is assumed to be an odd number. The N unknowns p m and q m are determined from N equations in Eq.(1.76). A system of N linear equations with N unknowns is obtained and solved numerically.

• Third step: Hypergeometric approximants in the Borel plane After calculating the N unknowns p m and q m in the previous step, the third step is to find the two vectors (x 1 , . . . , x l ) and (y 1 , . . . , y l ) solutions of the following two equations:

l m=0 p m x m = 0 (1.77) 1 + l m=1 q m y m = 0 (1.78)
We call those latter hypergeometric vectors. This is derived from the definition of hypergeometric functions. The hypergeometric vectors thus found identify in a unique way the hypergeometric function l+1 F l [START_REF]NIST Digital Library of Mathematical Functions[END_REF]:

B N (ξ) ≡ l+1 F l x, y, p l q l ξ (1.79)
where x = (1, -x 1 , . . . , -x l ), y = (-y 1 , . . . , -y l ) and B N (ξ) is the hypergeometric approximant in the Borel plane.

• Fourth step: Meijer-G approximants Finally, one must return to the original series in the variable t and eliminate the 1/n! by the Laplace transform of Eq. (1.79). This gives the sum of Meijer-G denoted by S M G (t) represented as:

S M G (t) = ∞ 0 e -ξ B N (tξ)dξ (1.80)
The Laplace transform of any hypergeometric function can be calculated analytically using Meijer-G functions. This yields to the following sum of Meijer-G approximant:

S M G (t) = l i=1 Γ(-y i ) l i=1 Γ(-x i ) G l+2,1 l+1,l+2   1, -y 1 , . . . , -y l 1, 1, -x 1 , . . . , -x l - q l p l t   (1.81)
where Γ(x) is the Euler Gamma function and G p,q m,n   c 1 , . . . , c p d 1 , . . . , d q z   is the Meijer-G function [START_REF] Beals | Meijer G-functions: a gentle introduction[END_REF] with c 1 , . . . , c p and d 1 , . . . , d q are the parameters of the function.

The hypergeometric ansatz and the resulting hypergeometric vectors are all that are required to parametrize these extremely complex functions. After determining the latter, the Meijer-G approximants can be parameterized for arbitrarily large truncation orders.

Recall 1.4.1

The Meijer-G function has the following integral representation [START_REF] Bateman | Higher transcendental functions[END_REF]:

G p,q m,n   c 1 , . . . , c p d 1 , . . . , d q z   = 1 2πi L n k=1 Γ(s -c k + 1) m k=1 Γ(d k -s) p k=n+1 Γ(-s + c k ) q k=m+1 Γ(s -d k + 1) z s ds (1.82)
where z = 0 and m, n, p, q are integer numbers such that 0 ≤ m ≤ q and 0 ≤ n ≤ p. L is the path of integration that can be from -i∞ to +i∞ [START_REF] Bateman | Higher transcendental functions[END_REF]. Note that Meijer-G function is an analytic function of z. This integral is called the Mellin-Barnes integral.

Remark 1.4.2

The idea of generalizing the MG algorithm to an even truncation order N can be found in [START_REF] Mera | Fast summation of divergent series and resurgent transseries from Meijer-G approximants[END_REF], [START_REF] Antipin | Resummation in QFT with Meijer G-functions[END_REF]. It is straightforward to generalize to an even number of N. Simply subtract the constant term from the original series, factor the linear component, and then apply the technique for odd N as described previously. Finally, the result must be re-multiplied with the linear term and the re-added constant. For example, suppose that the truncation order N is an even number, for example consider N = 4, therefore:

u(t) = u 0 + u 1 t + u 2 t 2 + u 3 t 3 + u 4 t 4 (1.83) = u 0 + u 1 t 1 + u 2 u 1 t + u 3 u 1 t 2 + u 4 u 1 t 3 (1.84)
We apply the algorithm previously described for odd-N on 1 +

u 2 u 1 t + u 3 u 1 t 2 + u 4 u 1 t 3 , then
we multiply the result by u 1 t and finally we add u 0 to the final result.

Finally, it should be noted that a recent work in 2020 proposed a resummation algorithm based also on Meijer-G but without using the Borel or Padé techniques. The suggested algorithm in [START_REF] Shalaby | Precise critical exponents of the O(N )-symmetric quantum field model using hypergeometric-Meijer resummation[END_REF], [START_REF] Shalaby | Weak-coupling, strong-coupling and large-order parametrization of the hypergeometric-Meijer approximants[END_REF] uses hypergeometric functions that have a zero radius of convergence by summarizing them using the equivalent integral representation of the Meijer-G function as an approximation of the given perturbation series.

Advantages/Disadvantages of MG

The advantages and disadvantages of MG that will be listed here are deduced from the examples that have been discussed in [START_REF] Mera | Fast summation of divergent series and resurgent transseries from Meijer-G approximants[END_REF] in the field of quantum mechanics. Thus, we recall that these advantages led us to treat the Meijer-G scheme for temporal cases. We will check whether we would benefit from the same advantages for temporal cases.

Many advantages were obtained by using the Meijer-G approximants in QFT fields. A brief review of the advantages is recalled here [START_REF] Mera | Fast summation of divergent series and resurgent transseries from Meijer-G approximants[END_REF]:

(i) The parameterization of Meijer-G approximants is simple. The hypergeometric ansatz provides a fast parameterization of Meijer-G approximants. The ratios between consecutive Borel transform coefficients are approximated by a rational function. This then allows defining of a hypergeometric function in the Borel plane. The Laplace transform of any hypergeometric function can be calculated analytically using Meijer-G functions. (ii) The Meijer-G approximants are more precise than the Borel-Padé. The examples presented in [START_REF] Mera | Fast summation of divergent series and resurgent transseries from Meijer-G approximants[END_REF] have shown that the Meijer-G approximant converges to the exact Borel sum in various QFT examples. (iii) Meijer-G functions also provide an approximation to the generalized Borel sums.

Generalizations of the Borel sum can be easily constructed using Meijer-G functions.

Although the use of the Meijer-G approximant has obvious advantages, it also has some disadvantages that should be taken into consideration. One may cite the following:

(i) Since Meijer-G approximants are new, their convergence properties are not yet developed and generalized for all truncation orders while the properties of the Padé approximants are well known. (ii) For large truncation orders, Meijer-G approximants might be difficult to evaluate numerically. (iii) There are difficulties with the numerical access to the Meijer-G function which does not exist in all the programming languages.

We now address the new part, compared to what has been done in the previously listed references, for the Meijer-G approximant:

1. Meijer-G for temporal test cases. 2. Meijer-G for vector series. 3. Meijer-G associated with the finite difference and finite element methods.

We start with the first point. The third point will be explored in the following chapters in the numerical applications.

Numerical algorithm MG used as a time integration scheme

In this section, we transform the numerical algorithm Meijer-G described previously into a time integration scheme. Therefore, all steps remain the same, except the last step where the perturbation parameter becomes t. Therefore, the sum of Meijer-G approximant becomes:

S M G ( t) = l i=1 Γ(-y i ) l i=1 Γ(-x i ) G l+2,1 l+1,l+2   1, -y 1 , . . . , -y l 1, 1, -x 1 , . . . , -x l - q l p l t  (1.85)
As already mentioned for the other methods, by applying the Meijer-G sum Eq.(1.85) just once, it may not be possible to reach S M G ( t) for a large value of t up to a final time t. Hence the definition of the notion of the validity domain [0, t max ] for the Meijer-G sum. Note that t max of Meijer-G is not yet evaluated explicitly. Therefore, a continuation procedure is proposed to reach a given specific time t. We propose the continuation of the MG algorithm by using a criterion based on the residual.

The residual Res(t) by using MG, is defined by:

Res(t) = M SMG (t) + C ṠMG (t) + KS M G (t) + N (S M G (t)) (1.86)
where ṠMG (t) the first derivative of Meijer-G sum and SMG (t) is the second derivative of Meijer-G sum. Therefore the approximations of derivatives of the MG algorithm are required for the residual evaluation replaced in the initial equation of the problem (Eq.(1.1)).

Here it should be mentioned that to calculate the derivatives, there are (as for the Padé approximant) at least two ways of doing it: 1. Calculate the derivatives of the series and establish an identification between the terms. (see remark 1.1.2). 2. Establish an analytical expression of the derivatives of the Meijer-G sum. The latter are established in Appendix C. We will see in the numerical example which one is the the most appropriate.

The continuation procedure of the MG algorithm used as a time integration scheme is summarized by:

Algorithm 4 Continuation technique for MG

Inputs:

1. The initial conditions u 0 and u 1 .

2. The truncation order N . 3. The residual tolerance or the accuracy parameter . 4. The final time t end or the domain in which we want to obtain the curve, for example: [0, t end ].

Step 1: Compute series terms until a given order N .

Step 2: Apply the four steps of MG to calculate the Meijer-G sum which is S M G (t).

Step 3: Apply the first and second derivatives approximations of MG to calculate ṠMG (t) and SMG (t). (See Appendix C).

Step 4: Evaluate the residual Res(t) given in Eq.(1.86).

Step 5: Evaluate the validity domain t max such that t max is the last value of t for which Res(t) 2 < .

Step 6: Evaluate a new starting point u 0 at a new time t 0 = t 0 + t max .

Step 7: Go to Step 1 and repeat the algorithm until reaching the final time t end .

For the Step 5 of the Algorithm 4 of the continuation technique of Meijer-G, it can be noted that the new starting point u 0 is u(t max ). To be able to apply Meijer-G algorithm with a normalized coefficient u 0 as described in subsection 1.4.1, one can proceed as follows: Firstly, we have:

u(t) = u 0 + u 1 t + u 2 t 2 + • • • + u N t n (1.87)
Secondly, we multiply and divide the right hand side of Eq.(1.87) by u 0 to obtain the first normalized coefficient. Therefore,

u(t) = u 0 u 0 + u 1 u 0 t + u 2 u 0 t 2 + • • • + u N u 0 t n u 0 (1.88) = 1 + u 1 u 0 t + u 2 u 0 t 2 + • • • + u N u 0 t n u 0 (1.89)
Thirdly, we apply the Meijer-G algorithm previously described on 1 +

u 1 u 0 t + u 2 u 0 t 2 + • • • + u N u 0 t n ,
and finally we multiply the obtained result by u 0 . This is also valid for the case where u 0 is not 1.

New version of MG approximant for a series of vector

This subsection is dedicated to the new part that we propose, regards introducing the Meijer-G approximant for a series of vectors. Suppose now that

u(t) = N n=0 u n t n (1.90)
where u n is either a vector of R n or a continuous field. We will take inspiration from what has already been done for the Padé approximants for a vector series [START_REF] Elhage-Hussein | A numerical continuation method based on Padé approximants[END_REF] and we will make the necessary adaptations to apply successfully Meijer-G to a vector series.

Example: First, we restrict the details of the development to an expansion of order N = 5. The idea is to orthogonalize the vectors u 1 , u 2 , u 3 , u 4 , u 5 by constructing an orthogonal basis using the classical Gram-Schmidt procedure in the following form:

u 1 = α 11 u * 1 , (1.91
)

u 2 = α 21 u * 1 + α 22 u * 2 , ( 1
u M G (t) -u 0 = S 1 M G (t)tu * 1 + S 2 M G (t)t 2 u * 2 + S 3 M G (t)t 3 u * 3 + S 4 M G (t)

Generalization:

We will generalize what we have just introduced for an expansion of order N = 5 to a given order N by providing the following two steps:

• The Gram-Schmidt method is used to define a set of orthonormal vectors:

u 1 = α 11 u * 1 u i = i j=1 α ij u * j , ∀i ≥ 2 (1.99)
The Gram-Schmidt coefficients α ij , (1 ≤ j ≤ i) are qualified by the orthonormality condition:

u * i , u * j = δ ij (1.100)
where δ ij is the symbol of Kronecker.

• By replacing Eq.(1.99) in Eq.(1.90), one can obtain:

u M G (t) -u 0 = N i=1 t i f i (t)u * i (1.101)
where the f i (t) are polynomials of degree N -i calculated from the Gram-Schmidt coefficients:

f i (t) = N -i k=0 t k α i+k,i (1.102)
We replace each of these polynomials by their corresponding Meijer-G sum denoted by S i M G (t). Therefore, the truncated series Eq.(1.101) has been replaced by a new Meijer-G representation which is of the following form

u M G (t) -u 0 = N i=1 t i S i M G (t)u * i (1.103)
Remark 1.4. [START_REF] Courant | Über die partiellen Differenzengleichungen der mathematischen Physik[END_REF] We would like to make it clear that in the rest of the manuscript, we will not apply Meijer-G resummation to all the examples we deal with in this manuscript, because we have already implemented some numerical tests in languages where the numerical access of the Meijer-G function G p,q m,n is not available (e.g. the C language). For example, the Meijer-G resummation will not be applied to the level set function in Chapter 4 because the latter was written in the C programming language.

Having outlined the methods, developments and our proposals, we now proceed to their numerical applications to different examples of ODEs and PDEs in the following chapters.

Chapter 2

APPLICATION ON NONLINEAR ODES

After exposing the time perturbation and summation approaches as time integration schemes in Chapter 1, we seek to validate the latter on various nonlinear ODE examples. This chapter aims to understand the behavior and the feasibility of these integrations, as well as to validate the first and second derivatives of BPL and IFS that we have established, and finally to consolidate the previous results. This leads us to highlight the potential of time perturbation approaches compared to the classical schemes by providing examples where the latter suffers from many issues. The properties and qualities of IFS will be studied in detail since the latter has already been applied as a first prospective study on some simple examples. The Meijer-G (MG) algorithm used in temporal cases will be studied on two ODE examples in this chapter. We recall that MG will not be applied to all the examples because of an access difficulty of the Meijer-G function in some programming languages.

In this chapter, different examples of nonlinear ODEs will be presented. First, we will consider again the popular Euler equation in order to validate and study the behavior of MG on this example. Then, the simple harmonic oscillator will also be studied by MG to validate the proposal of studying MG in temporal cases. BPL, IFS, and ANM will be used for further examples. The details will be performed only for a chosen solid mechanics problem which is the free and forced Van der Pol oscillators. Then, the proposed techniques are applied to the nonlinear combustion equation, the hardening spring problem, and the nonlinear simple pendulum problem. These three examples aim to investigate the behavior of these methods for problems presenting large nonlinearity where classical explicit schemes suffer from different constraints. Finally, the efficiency of these methods is also established for the most famous chaotic system which is the Lorenz system (three degrees of freedom DOF example). The goal of this example is to check whether the temporal approaches can respect the different phase portraits obtained by the Lorenz system and to validate their solving ability for a 3 DOF example before dealing with problems involving several DOFs in the next chapter.

Since we have a priori no information on the convergence type of the obtained numerical series, the BPL and IFS summation techniques are also applied. Therefore, even if the obtained series are convergent, we will use resummation techniques to observe their interest on a convergent series and study their accuracy and robustness.

For a fair comparison: the accurate first, third and fourth-order Runge-Kutta (RK) methods using fixed or adaptive time steps schemes are considered in this chapter.

Application to the Euler equation

Consider again the Euler equation

t 2 u + u = t (2.1)
with the following initial condition

u(t 0 = 0) = 0 (2.2)
The formal power series solution of this equation is

u(t) = ∞ n=0 (-1) n n!t n+1 (2.3)
which is a divergent series.

As the resulting series is divergent, this example is important to be tested by the summation algorithms. Since BPL has been applied to this example in Chapter 1 (subsection 1.1.3) and in [START_REF] Razafindralandy | Time integration algorithm based on divergent series resummation, for ordinary and partial differential equations[END_REF], [START_REF] Deeb | Intégrateurs temporels basés sur la resommation des séries divergentes[END_REF], the goal of this example is to present a first validation of the Meijer-G algorithm used in the temporal case. To be able to apply the Meijer-G approximants, the initial conditions u 0 cannot be zero. The series are calculated up to a truncation order N = 5. We compare the solutions obtained by the MG, ANM, BPL, and IFS without any continuation in Fig. (2.1). The solution obtained by the Meijer-G approximant is in good agreement with the exact solution and presents a better approximation of the solution than ANM, BPL, and IFS. Now, we maintain the truncation order N = 5 for Meijer-G (MG) and increase it to N = 10 for ANM, BPL, and IFS. Even with a higher truncation order than the Meijer-G approximant, MG provides a better solution than the other methods as can be seen in Fig.

(2.2) using a small truncation order. We recall that the general definition of the Meijer-G function is given by a line integral in the complex plane, hence the possibility to obtain a complex result. (See the definition of the Meijer-G function in Chapter 1). We note that we are interested in the real part of the solution. The obtained imaginary part of the Meijer-G solution of the Euler equation is zero.

Application to the simple harmonic oscillator

Consider the following form of the simple harmonic oscillator given by ü + u = 0 (2.4)

where u = du dt and ü = d 2 u dt 2 , with the following initials conditions

   u(0) = 1 u(0) = 2 (2.5)
which has an exact solution for t 0 = 0 on the form:

u(t) = u(0) cos(t) + u(0) sin(t) (2.6)
Series terms u k , for k ≥ 2 are given by the following recurrence formula

k(k -1)u k = -u k-2
The Meijer-G approximant is applied to this example to test the effectiveness of the proposal way to use it in a temporal case. The series are calculated up to a truncation order N = 5. As can be seen in Fig. (2.4), the imaginary part of the solution obtained by Meijer-G solution is zero until t = 1s. We note that in some cases, it is quite simple to give a physical meaning to the imaginary part. In this case, it is not so obvious and the physical interpretation is not our main objective. Our goal is to obtain a precise approximation of the solution. In the following, we just plot the real part of the Meijer-G solution that interests us.

Application to the free Van der Pol oscillator

In this section, we are interested to study the effectiveness of time perturbation and resummation methods on the free Van der Pol oscillator problem. The details of the study will be presented only for this example. The Van der Pol oscillator has attracted considerable interest from many researchers for a long time since it models many problems in different areas [START_REF] Marios | Theoretical and Numerical Study of the Van der Pol equation[END_REF]. It is a dynamic system described by a variable x(t) satisfying the following nonlinear differential equation of the form

ẍ -µ(1 -x 2 ) ẋ + x = 0 (2.7)
where ẋ = dx dt and ẍ = d 2 x dt 2 , with the following initials conditions

   x(0) = a ẋ(0) = b (2.8)
where a and b are two parameters, and µ ≥ 0 denotes the control parameter or called also the coefficient of nonlinearity. Therefore, the nonlinear behavior of the oscillations depends on µ. Therefore, according to its value, the trajectories have precise characteristics.

For µ = 0, equation Eq.(2.7) will be the simple harmonic oscillator addressed in section 2.2. For µ = 0, equation Eq.(2.7) has no exact solution, therefore it was first integrated by the classical time discretization schemes, to compare the efficiency of time perturbation methods to the classical ones.

Classical discretization schemes

Equation Eq.(2.7) is written as a first order system of differential equation:

   ẋ = y ẏ = µ(1 -x 2 ) ẋ -x (2.9)
A convergence study for different time steps ∆t between Runge-Kutta of 1-st (RK1) and 4th order (RK4) has been carried out. For sake of simplicity, only results obtained with RK4 are presented for this example. The goal of this comparison is to get a good approximation of the solution since it will be the reference solution.

Firstly, constant time step ∆t = 10 -3 s is chosen for RK4. Secondly, a type of adaptive RK4 has been used as presented in [START_REF] Press | Adaptive Stepsize Runge-Kutta Integration[END_REF]. It computes a local error between two solutions obtained after performing two consecutive time steps using step size ∆t and another using only a one-time step but with a larger step size 2∆t. The tolerance parameter on local error denoted ∆ 0 , makes it possible to adjust the step size. Initial conditions at a given time t 0 = 0 are arbitrarily chosen as x(0) = 1 and ẋ(0) = 0 and the coefficient of nonlinearity as µ = 2.

It should be noticed that for every choice of initial conditions, except {x(0) = 0, ẋ(0) = 0}, a unique periodic motion is presented.

The results of numerical integration by Runge-Kutta methods show that the characteristic of the trajectory depends on the value of the coefficient of nonlinearity µ are presented in Fig. (2.7). For small µ, the limit cycle is nearly circular. For large µ, the limit cycle is no longer nearly circular and is representative of bounded periodic oscillatory behavior.

Time perturbation series recurrence formula

In this subsection, time perturbation and resummation methods are applied to the free Van der Pol oscillator. We seek to demonstrate the accuracy of these methods by comparing their results to RK4, check their ability on a problem involving significant temporal variations and limit cycles, and consolidate some results previously achieved. Consider again the Van der Pol equation

ẍ -µ(1 -x 2 ) ẋ + x = 0 (2.10)
Unknowns of this equation x(t), ẋ(t) and ẍ(t) are sought as time power series

x( t) = x 0 + N i=1 x i ti ẋ( t) = N i=1 ix i ti-1 (2.11) ẍ( t) = N i=2 i(i -1)x i ti-2
where N denotes the truncation order of the series. When we replace the expansion Eq.(2.11) into equation Eq.(2.10), and identify according to the power of t, we notice that it is not easy to find the recurrence formula that allows to calculate the series term. For this reason, let the following change of variable

z = x 2 (2.12)
This yields to the new quadratic equation

ẍ -µ(1 -z) ẋ + x = 0 (2.13)
Writing z( t) as time power series

z( t) = z 0 + N i=1 z i ti (2.14)
This yields to the following sequence problems

                         t0 : 2x 2 -µ(x 1 -z 0 x 1 ) + x 0 = 0 t1 : 6x 3 -µ(2x 2 -2x 2 z 0 -x 1 z 1 ) + x 1 = 0 t2 : 12x 4 -µ(3x 3 -3x 3 z 0 -2z 1 x 2 -z 2 x 1 ) + x 2 = 0 . . . tN : (N + 2)(N + 1)x N +2 -µ (N + 1)x N +1 - N i=0 (N + 1 -i)x N +1-i z i + x N = 0
Therefore, series terms x k for k ≥ 2 are given by the following recurrence formula

k(k -1)x k = µ (k -1)x k-1 - k-2 i=0 (k -1 -i)x k-1-i z i -x k-2
where

z i = i j=0 x j x i-j (2.15)
Once series terms x k are numerically evaluated, range of validity t max of the perturbed representation Eq.(2.11) is sought. To do so, two techniques are used in this study. The first uses Eq.(1.17), when using this criterion the method will be noted ANM(C). For the sake of comparison with other perturbation methods, the residual criterion will be used to evaluate this t max (see remark 1.3.6). This method is noted by ANM(R) or simply ANM.

Numerical study of the two possible ways to evaluate the derivatives of the IFS sum

We recall that in Chapter 1, we proposed two possible ways to evaluate the derivatives of the Borel sum, the IFS sum and the Meijer-G sum which allow to calculate the residual and thus to apply the continuation based on the latter.

First possibility:

It is based on calculating the initial series terms, the first and second derivative series terms, and finally applying the resummation algorithm for each of these terms. Therefore, this possibility is related to the three following steps:

• Calculate the series terms denoted by x k , ∀k = 1, • • • , N corresponding to each problem, then apply the resummation algorithm to these terms for a truncation order equal to N .

• Calculate the first derivative series terms denoted by ẋk , ∀k = 1, • • • , N -1 that are linked to the initial series terms x k by the following relation:

ẋk = (k + 1)x k+1 ∀k = 1, • • • , N -1 (2.16)
then apply the resummation algorithm to these new terms ẋk , for a truncation order equal to N -1.

For a more detailed explanation of how Eq.(2.16) is obtained, see the remark 1.1.2 in Chapter 1.

• Calculate the second derivative series terms denoted by ẍk , ∀k = 1, • • • , N -2 that are linked to the initial series terms x k by the following relation:

ẍk = (k + 1)(k + 2)x k+2 ∀k = 1, • • • , N -2 (2.17)
then apply the resummation algorithm to these new terms ẍk , for a truncation order equal to N -2.

For a more detailed explanation of how Eq.(2.17) is obtained, see the remark 1.1.2 in Chapter 1.

As a result, this way requires three serial term calculations (x k , ẋk , and ẍk ) and three resummation algorithms applied for each of these. Since this possibility is based on a series derivation, it will be called "Series.D" in what follows in this section.

Second possibility: It consists of calculating the initial series terms once, then applying the resummation algorithm to these terms, and finally establishing the expression of the first and second derivatives to the sum formulas obtained through the resummation algorithms. For more details, this possibility is based on the following three steps:

• Calculate the series terms corresponding to each problem x k , ∀k = 1, • • • , N , then apply the resummation algorithm to these terms.

• Establish an analytical expression for the first derivatives of the BPL, IFS, and MG sum determined respectively in Appendices A, B, and C.

• Establish an analytical expression for the second derivatives of the BPL, IFS, and MG sum determined respectively in Appendices A, B, and C.

As a result, this way requires one calculation of series terms, one resummation algorithm, and two approaches of first and second derivation applied to the summation formulas of BPL, IFS, and MG. Since this possibility is based on a resummation derivation, it will be called "Resummation.D" in what follows in this section.

Therefore, it is quite natural to wonder about the difference between "Series.D" and "Resummation.D". The questions that can arise are:

-Which one is more appropriate to be applied? -Which one is more suitable for long-time simulation?

-Which one is faster in terms of CPU times and requires fewer time steps to reach the same final time? Hence, we will indicate the answers through some numerical tests applied to the example of Van der Pol by the Inverse factorial series (IFS) with continuation. Therefore, for "Resummation.D", Appendix B is used (first and second derivatives approximations of the exact sum of IFS).

Test 1: IFS with continuation to a final time t = 10s using a small value of the coefficient of nonlinearity µ = 2

The first test is a continuation of IFS for a final time t = 10s using a small value of the coefficient of nonlinearity which is fixed to µ = 2. The two ways previously mentioned to calculate the derivatives and then the residual, are used for the continuation technique. The truncation order is fixed to N = 15 and the residual tolerance to = 10 -4 for the both ways. The reference solution is RK4 with a discrete-time step ∆t = 10 -3 s. The solutions of displacement (x(t)), velocity ( ẋ(t) or the first derivative of the displacement), and acceleration (ẍ(t) or the second derivative of the displacement) obtained by IFS are compared to the reference solutions in Fig. (2.8). A good agreement of the IFS solutions (using both ways) with the reference solutions can be seen in the Fig. (2.8), which shows the effectiveness of "Series.D" and "Resummation.D" for this test.

After validating the accuracy of "Series.D" and "Resummation.D", we are now interested in comparing the number of time steps required by both ways to reach the same final time t = 10s. As can be observed in Tab.(2.1), "Series.D" requires less steps than "Resummation.D" to reach t = 10s. Another test is presented afterwards to check whether these results hold true. The second test is a continuation of IFS for a final time t = 1000s using a small value of the coefficient of nonlinearity µ = 2. The truncation order and the residual tolerance are maintained respectively to N = 15 and = 10 -4 .

It can be seen in Fig. (2.9) that the solutions of displacement, velocity, and acceleration obtained by IFS using "Series.D" show a drift with the reference solution, whereas the solutions obtained using "Resummation.D" are in good agreement with the reference solution. We recall that for the sake of comparison, the same numerical parameters are used for the both methods. This test indicates that the application of "Series.D" over a long simulation period time may present a risk of phase shift with respect to the reference or the exact solution. This may create other significant phase differences at longer intervals.

63

Étude numérique de méthodes d'intégration temporelle basées sur des techniques de perturbation. : application à des problèmes de mécanique Cynthia Tayeh 2022 The residual tolerance is now decreased from = 10 -4 to = 10 -5 for the "Series.D" to check if it allows to obtain a faithful solution as "Resummation.D" and the reference solution, for a long final time. After various studies on the range of tolerances, it was found that "Series.D" with a residual tolerance = 10 -6 provides an efficient solution compared to the reference solution as can be appreciated in Fig. (2.10). Regarding the number of steps, IFS using "Series.D" and a residual tolerance = 10 -6 is compared to the IFS using "Resummation.D" with a residual tolerance = 5 × 10 -4 (since these two optimal choices of the residual tolerance provide an efficient solutions for each way, until a final time t = 1000s). As can be seen in Tab.(2.2), the number of steps of "Resummation.D" is in the same order of magnitude than that of "Series.D" for this test.

Application on nonlinear ODEs

IFS with "Series.D" and = 10 -6 IFS with "Resummation.D" and = 5 × 10 -4 10617 10223

Table 2.2 -Number of steps required to reach t = 1000s for the free Van der Pol oscillator for µ = 2, using time perturbation-resummation method IFS with N = 15 with different values of the residual tolerance .

Test 3: IFS with continuation to a final time t = 1000s using a large value of the coefficient of non linearity µ = 1000

We have already addressed the case of a small nonlinearity coefficient for a final time t = 10s in Test 1, and t = 1000s in Test 2. In this third test case, we increase the value of the nonlinearity coefficient to µ = 1000. The goal of this test is to study the behavior of using a continuation of IFS via the derivatives of "Series.D" and "Resummation.D" for a problem presenting large nonlinearity. The simulation is carried out until a final time t = 1000s. We will only display in the figures the part where there is a jump in the solution, which is equivalent to the time interval [807.38, 807.385]s. As can be seen in Fig. (2.11), there is a shift between the reference solution that is either the RK4 with ∆t = 10 -5 s or the ANM solution (this part will be described in detail in subsection 2.3.5) and the IFS with "Series.D" with a residual tolerance = 10 -4 . With the same choice of , "Resummation.D" is able to provide an efficient solution.

-1 -0.5 0 0.5 Indication of this study: After these different test cases, one can in principle state that the application of the IFS continuation using the exact approximations of the deriva-tives in Appendix A (named here by "Resummation.D") guarantees a good quality of the solution for different nonlinearity parameters of the problem. Moreover, the solution obtained by the latter is also efficient for long-time intervals. Calculate the derivatives of the series and establish an identification between the terms (named here by "Series.D") may also provide an efficient solution by using fewer time steps than "Resummation.D" to reach the same final time in some cases. However, one must be careful because "Series.D" may cause some shifts in the solution obtained for simulations over long time intervals. As a result, the residual tolerance must be reduced, which increases the number of steps required to reach the final time and, as a result, the CPU time too. This may also lead to the possibility of a large shift in the solution over time. It is expected that the same results occur with the other resummation methods BPL and MG using "Series.D" and "Resummation.D" based on Appendix A for Borel sum and Appendix C for Meijer-G sum. In addition, it should be noted that similar numerical tests were carried out by using BPL. The results revealed that numerical difficulties could appear in particular in the numerical choice of the truncation orders of the 3 series (the initial series, the first derived series, and the second derived series) to avoid the appearance of poles of the Padé approximations. Therefore, after this simple study, we note that for the rest of the manuscript we will only apply the second possibility ("Resummation.D") to calculate the residual of BPL, IFS, and MG and finally to perform an efficient continuation. Therefore, Appendices A, B and C will be used.

Numerical results for time perturbation methods

Initial conditions are set to x(0) = 1 and ẋ(0) = 0, and the coefficient of nonlinearity is fixed in this subsection to µ = 2.

Numerical results of ANM

Firstly, ANM behavior is presented without any continuation procedure. Effects of the truncation orders N are depicted in Fig. (2.12). Continuous ANM solutions are compared to the discrete solutions obtained using RK4. It is observed that ANM can reproduce correct solutions for a large range of time when compared to RK4. Then the continuation technique is used. The truncation order is fixed to N = 15 and the tolerance of ANM(C) is set to δ = 10 -6 .

Remark 2.3.1

We note that we have not represented RK4 with the cross symbols as it will produce many points due to the solution for each discrete-time step which can make the figure less clear. It should be noted that in between each dot on the plot, solutions are continuous. This is one of the main features of such an approach.

Continuation is depicted in Fig. (2.13). It also highlights the accuracy of the continuation of the time series, since it presents the same numerical result as RK4.

Figure 2.12 -Approximated solution using time perturbation method ANM (series) using different truncation order without continuation, for the free Van der Pol oscillator for µ = 2. In black the reference RK4 solution using ∆t = 10 -3 s. 2.14b) that is outside the limit cycle is not a deviation of the solution, but is just the starting part of the solution before reaching the limit cycle.

Numerical results of IFS

It is proposed to use IFS in the case of N = 8, µ = 2. Poles appear for BPL method using those parameters (see Fig.

Simulation over a large time

Finally, the simulations are executed by ANM, BPL, and IFS up to a final time t = 10000s. Numerical results show the ability of time perturbation and resummation methods to reproduce the dynamics of the system over a large time. The limit cycle is respected by these three methods even for a large final time. The effiency of the system dynamic obtained by BPL is presented in Fig. (2.18) for the interval time [9900; 10000]s. It should be noted that the long time behavior of BPL is investigated also through numerical experiments on the Hamiltonian and non-Hamiltonian systems in [START_REF] Razafindralandy | Some robust integrators for large time dynamics[END_REF]. A second fully detailed article from the same research team is in progress on the stability of the BPL scheme (see https://arxiv.org/pdf/2007.00390.pdf). As IFS used as a time integration scheme has never been applied for long time. This example can be used to illustrate the capability of IFS for long time simulations.

ANM, BPL and IFS: t max comparisons

We now propose a comparison of ranges of validity t max . For the sake of comparison, ranges of validity are evaluated using the residual criterion.

Evolution of t max obtained using ANM, BPL and IFS is proposed in Fig. (2.19). It is observed an automatic adaptation of the range of validity [0, t max ] accordingly to the cyclic solution. Ranges of validity of ANM and BPL oscillate around 0.4s and stay bounded in [0.1, 0.78]s and is larger than that of IFS.

ANM, BPL, and IFS: Number of continuation steps comparisons

We now propose a comparison of the number of continuation steps needed to reach a specific final time. We recall here that for each continuation step, continuous approximated solutions are obtained using time perturbation methods. For the sake of comparison, ranges of validity are evaluated using the residual criterion. As a reminder, RK4 needs 10000 steps to reach a final time of t = 10s using a discrete-constant time step of ∆t = 10 -3 s.

A small study on number of step is presented in Tab.(2.3) for adaptive RK4 [START_REF] Press | Adaptive Stepsize Runge-Kutta Integration[END_REF]. We recall that adaptive RK4 computes a local error between two solutions obtained after performing two consecutive time steps using step size ∆t and another using only a onetime step but with a larger step size 2∆t. The tolerance parameter on local error denoted ∆ 0 , adjusts the step size. The number of steps for ANM, BPL, and IFS with each other is given in Tab. (2.4). It is noted that ANM and BPL reach the final time using approximately the same number of steps. BPL is slightly better than ANM in this example. IFS performs two to four times more steps than ANM and BPL. Interestingly, these 4 methods need fewer steps than the classic and the adaptive RK4 in Tab. (

Computational efforts

To the best of our knowledge and as already mentioned in Chapter 1, no explicit formula of temporal range of validity (t max ) of BPL and IFS has been established yet. Therefore, it is proposed to evaluate it via a residual criterion as in [START_REF] Razafindralandy | Time integration algorithm based on divergent series resummation, for ordinary and partial differential equations[END_REF]. The latter guarantees a very good quality of the solution which is continuous but sometimes timeconsuming because of a large number of numerical residue evaluations. Some global information on the computational effort is given in this subsection. For this test case, the comparison is only an indication since the Van der Pol oscillator is considered a problem of small size.

We compare the performance of the numerical methods in terms of their computational efficiency (as measured both by computational load i.e the CPU time using here the GNU Octave) and in terms of numerical accuracy (as quantified by evaluating the residual of the computed numerical solution) for a final time of t = 100s.

As a reference, RK4 using a constant time step ∆t = 10 -3 s needs 7.8709s of CPU time to reach a final time of t = 100s. When the RK4 adaptive time step is used with ∆ 0 = 10 -6 , the CPU time needed is 10.258s.

Remark 2.3.3

We just compare the CPU time given for ∆ 0 = 10 -6 since this choice of ∆ 0 allows the residual of RK4 to stay in the same residual range as those of the perturbation methods.

CPU times for temporal perturbation approaches are given in Tab. (2.5). ANM(C) with δ = 10 -5 , is in the same residual ranges as other methods. Complete computation is performed in 0.29s of CPU time. It should be noted that ANM(R) (ANM using residual norm to evaluate t max ) and BPL require more CPU time than ANM(C). It is noted that IFS is relatively less efficient in terms of computational time as can be seen in Tab.(2.5). Interestingly, ANM(C), ANM(R), BPL, and IFS are faster than RK4 adaptive and RK4 with constant step size. Therefore, the results show that ANM(C) with δ = 10 -5 is the fastest in CPU time while the residual stays in the same range as others methods. ANM(C) is 4 time faster than ANM(R), however, ANM(R), BPL, and IFS ensure a prescribed residual tolerance. Moreover, most of the CPU time is spent in the evaluation of the range of validity.

Method CPU time (s) ANM(C) 0.29951 ANM(R) 1.3056 BPL(R) 3.9281 IFS(R) 7.0321 Adaptive RK4
10.258 Table 2.5 -CPU time needed to reach t = 100s using the RK4 adaptive time step for a local error tolerance ∆ 0 = 10 -6 , time perturbation method ANM, time perturbationresummation methods BPL and IFS for a truncation order N = 8. (R) stands for ranges of validity evaluated using a residual tolerance criterion = 10 -3 , while (C) stands for the criterion proposed by Cochelin [START_REF] Cochelin | A path-following technique via an asymptotic-numerical method[END_REF] using δ = 10 -5 .

Finally as conclusion regarding computational efforts and the number of steps for adaptive RK4 and time perturbation methods, it is possible to say that wider step sizes are performed, in a slightly better amount of time and for a better-prescribed accuracy of the solution.

Remark 2.3. [START_REF] Van Zuijlen | Implicit and explicit higher order time integration schemes for structural dynamics and fluid-structure interaction computations[END_REF] We note that only one of the methods will be displayed on the figures for the remaining test cases of this chapter (for clarity). However, we emphasize that we have investigated the efficiency of all methods: ANM, BPL, and IFS.

Behavior for a higher value of µ

In this subsection, we consider the same equation of the free Van der Pol oscillator by increasing the non-linearity coefficient µ to 1000. In this particular case, the solution changes radically and presents oscillations on a much longer time scale.

Goal of this example:

We are interested to see how the polynomial approaches can deal with large values of the nonlinearity coefficient where the equation becomes stiff and some very fast dynamic in the system is expected. It is observed that the use of the resummation methods BPL and IFS tend to increase the range of validity of the ANM series in this case. It can also be seen that the validity domain of the four methods is larger than the discrete-time step of RK4 which is 10 -5 . Regarding the adaptive time step of RK4, and as mentioned in [START_REF] Cartwright | The dynamics of Runge-Kutta methods[END_REF], if we integrate a stiff problem using adaptive Runge-Kutta, the initial step size chosen generally leads to a large local error estimate. This then causes the routine to reduce the step size until the main local error is reduced to its specified limit. The routine sometimes works well to integrate the problem, but uses a much larger number of steps than seems logical. As a result, rounding error and computation time are constraints when using classical techniques to integrate such problems. These constraints are surmounted with the time perturbation and resummation methods. This example is a very good illustration to realize that these methods seems to be appropriate for stiff problems.

Application to the forced Van der Pol oscillator

An external force can be simulated as acting on the system by adding a new periodic force to the second member of the Van der Pol equation [START_REF] Parlitz | Period-doubling cascades and devil's staircases of the driven Van der Pol oscillator[END_REF]. This forced equation is studied in some papers [START_REF] Hinvi | A Study of the Forced Van der Pol Generalized Oscillator with the Renormalization Group Method[END_REF]- [START_REF] Cooper | Controlling Chaos-Forced Van der Pol Equation[END_REF], and is in the following form:

ẍ -µ(1 -x 2 ) ẋ + x = F (t) (2.18)
This equation is called forced Van der Pol.

The external action will be assumed sinusoidal periodic force, therefore, the second member F (t) is chosen as A cos(ωt), where A denotes the forcing amplitude and ω the forcing frequency.

Goal of this example:

While we use time perturbation for the unknowns, the righthand side (RHS) F (t) also needs to be represented by a power series. To obtain the best numerical solution of the complete forced equation Eq.(2.18), it is necessary to have the best approximation of this RHS. Therefore, a time perturbation study of this RHS is realized by comparing the Taylor method and time perturbation methods using specific change of variables [START_REF] Guillot | A Taylor series-based continuation method for solutions of dynamical systems[END_REF] based on quadratic recast. It is interesting to note that there could be differences between these two methods in the sense of the range of validity due to the truncation order. One has to find out which of these two methods has a larger range of validity to reach the final time.

Time perturbation of the forcing term

The chosen time-dependent forcing term is F (t) = A cos(ωt). It is necessary to look for a polynomial approximation of it, as:

F ( t) = F 0 + N i=1 F i ti (2.19)
with F 0 = F (t 0 ). Two specific procedures are compared to evaluate the F i series terms: Variable changes and Taylor decomposition.

Variable change:

The following change of variables U = cos(ωt) and V = sin(ωt) proposed in [START_REF] Guillot | A Taylor series-based continuation method for solutions of dynamical systems[END_REF] leads to

   U = cos(ωt) V = sin(ωt) =⇒    dU = -ωV dt dV = ωU dt (2.20)
Writing U and V as series:

U = U 0 + N i=1 U i ti and V = V 0 + N i=1 V i ti , where U 0 = cos (ωt 0 )
and V 0 = sin (ωt 0 ) are the initial terms with t 0 a given initial time.

Injecting the latter into Eq.(2.20), and then equating terms with identical powers of t, we get for i ≥ 1:

   U i = -ω i V i-1 V i = ω i U i-1 (2.21)
Thus, the recurrence formula for the RHS, ∀i ≥ 0, is:

F i = AU i (2.22)
Taylor decomposition:

The general formula of Taylor series for cos(ωt) also makes it possible to determine the F i series terms:

cos (ωt) = N i=0 (cos(ωt 0 )) (i) i! (t -t 0 ) i (2.23)
where (cos(ωt 0 )) (i) denotes the ith derivative of the cosine function evaluated at the point (ωt 0 ). This formula has been directly used for numerical validation. The recurrence formula for the F i terms is not provided here as it was not necessary.

Numerical validation of the RHS series representation

The exact solution is compared to the solution given by the Taylor method and the change of variable method, for A = 17, ω = 4. It is done without a continuation procedure to fully understand the behavior of the obtained series. Firstly, influence of the initial time t 0 is depicted in Fig. (2.24) using truncation order N = 20. No difference can be seen between both methods, nor influence of the initial time. Secondly, truncation order influence is presented in Fig. (2.25). It can be seen that the domain of validity t max of the RHS, ≈ 1.5s, is wider than the previously computed branch of solutions.

Numerical results without continuation procedure show that the two methods: the Taylor formula and the method based on change of variable, have the same behavior, the same range of validity, and the same error with the exact solution. Since the proposed method leads to a much simpler recurrence formula Eq.(2.22) for evaluating the F i terms, we decide to use this method for the RHS.

Numerical study of the forced Van der Pol equation

Classical discretization schemes

The forced Van der Pol equation Eq.(2.18) is solved by converting it into a first order system of differential equation:

   ẋ = y ẏ = µ(1 -x 2 ) ẋ -x + F (t) (2.24)
A convergence study for different time steps ∆t between RK4 and RK1 schemes has been also carried out in this case. The corresponding time step for x(0) = 1, ẋ(0) = 0 and µ = 2 is ∆t = 10 -3 s.

Time perturbation representation

To solve the forced Van der Pol equation using time perturbation, we proceed with the same procedure as in the free case. Power series representation of the unknowns Eq.(2.11), using time as the perturbation parameter, is injected in the equation to solve Eq.(2.18). Moreover, series terms F i of the RHS are evaluated using Eq. (2.22). Therefore, by identifying the terms according to the powers of t, the following recurrence formula is obtained:

k(k -1)x k = µ (k -1)x k-1 - k-2 i=0 (k -1 -i) x k-1-i z i -x k-2 + F k-2 (2.25)
where

z i = i j=0 x j x i-j (2.26)

Numerical results

For all the computations, initials conditions are set to x(0) = 1, ẋ(0) = 0, the amplitude to A = 17 and the frequency to ω = 4.

To understand the behavior of the ANM, solutions are plotted without continuation procedure in Fig. (2.26). It can be seen as a standard evolution in the perturbation method. After some time, continuous solutions diverge. Moreover, the more terms added in the representation, the wider is the range of validity. Remark 2.4.1 The proposed time perturbation methods have also been successfully applied to the forced Van der Pol oscillator with two degrees of freedom. The example considered concerns a triode generator with an additional circuit which is the best-known example of a self-oscillatory system with two degrees of freedom. As it does not provide any new results, we made the choice not to add extra redundant information in the manuscript. Another model with three degrees of freedom has been studied in section 2.7 and the more next chapters of the manuscript.

Application to the nonlinear combustion equation

This example is studied by Larry Shampine, one of the authors of MATLAB's suite of ordinary differential equations [START_REF] Shampine | The Matlab ODE suite[END_REF]. This problem is characterized by its transitory aspect in the middle of the integration interval. The solution switches from non-stiff to stiff, and then it becomes non-stiff again.

Goal of this example:

As numerical drift is known to occur in such cases, we are interested in studying the performance of time perturbation and resummation methods in this example. We consider a model of flame propagation. If a match is lit, the fireball spreads rapidly until it reaches a critical size. It then preserves this size because the amount of oxygen absorbed by the combustion inside the ball balances the amount accessible through the surface of the fireball. To illustrate this further, let us consider the following nonlinear combustion equation [START_REF] Abelman | Comparison of some recent numerical methods for initial-value problems for stiff ordinary differential equations[END_REF]:

ẏ = y 2 (1 -y) (2.27)
where ẏ = dy dt , with the initial condition

y(0) = R (2.28)
The scalar variable y(t) represents the radius of the ball. The terms y 2 and y 3 refer to surface area and volume. The critical parameter R corresponds to the initial radius which is small. The solution is sought over a length of time such that t ∈ [0, 2/R].

The exact solution of this differential equation is [START_REF]The Math Works[END_REF]:

y(t) = 1 W (a exp(a -t)) + 1 (2.29)
where a = 1/R -1 and W (t) is the Lambert W function (which is the inverse of the function: t exp(t)).

Time perturbation series recurrence formula

Series terms of time perturbation method for equation Eq.(2.27) are given by the following recurrence formulas, ∀k = 1, • • • , N we have:

ky k = - k-1 i=0 z i y k-1-i
(2.30)

z i = i j=0 y j y i-j (2.31)

Numerical results

For the numerical test, we take R = 0.0001. We provide the solution of the combustion equation by time perturbation methods and MATLAB ode solvers. We start by ode45 (based on an explicit Runge-Kutta (4, 5) formula and the Dormand-Prince pair [START_REF] Dormand | A family of embedded Runge-Kutta formulae[END_REF]) since in general, it is the best solver to apply as a first try for the majority of problems. Note that good results can be obtained with MATLAB solver ode23s (based on a modified Rosenbrock formula of order 2 [START_REF] Shampine | The Matlab ODE suite[END_REF]). However, it only solves some kinds of stiff problems and none of all ode solvers is exempt from the step size restriction. Moreover, as temporal perturbation methods are explicit methods and explicit schemes are known to present some difficulties when solving stiff problems, we have shown that time perturbation-resummation methods solves well this kind of problem. Thus, the qualitative characteristics of the solution obtained by the time perturbation method are shown as well as the good performance of the solution up to a large final time (t = 20000s). This was proven by evaluating the error between the exact solution and the one obtained by the time perturbation method over time which is of order 10 -5 .

Application to the Elastic hardening spring problem

Consider the nonlinear elastic hardening spring problem [START_REF] Liu | An assessment of the differential quadrature time integration scheme for nonlinear dynamic equations[END_REF], [START_REF] Kim | Higher-order explicit time integration methods for numerical analyses of structural dynamics[END_REF] in the following form:

ü

+ s 1 (1 + s 2 u 2 )u = 0 (2.32)
where u = du dt and ü = d 2 u dt 2 , with the following initials conditions

   u(0) = u 0 u(0) = v 0 (2.33)
where s 1 > 0 and s 2 > 0.

Goal of this example:

This problem is well-known for being a great example for evaluating the performance of new temporal integration techniques as mentioned in [START_REF] Kim | Higher-order explicit time integration methods for numerical analyses of structural dynamics[END_REF]. It is an example of conservative systems which maintain constant total energy in their exact solutions. The solution to this problem presents periodic oscillations.

The exact period and solution can be found analytically in [START_REF] Wood | Stability properties of some algorithms for the solution of nonlinear dynamic vibration equations[END_REF].

Time perturbation series recurrence formula

Series terms of time perturbation method for Eq.(2.32) are given by the following recurrence formulas, ∀k = 1, • • • , N we have:

k(k -1)u k = -s 1 k-2 i=0 (1 + s 2 y i )u k-2-i
(2.34)

y i = i j=0 u j u i-j (2.35)

Numerical results

For numerical results, we take u 0 = 1.5, v 0 = 0.0, s 1 = 100, s 2 = 10 as presented in [START_REF] Kim | Higher-order explicit time integration methods for numerical analyses of structural dynamics[END_REF]. The nonlinear period of the problem is T = 0.15153. Classic time integration schemes use a discrete-time step ∆T = T /32 as in [START_REF] Kim | Higher-order explicit time integration methods for numerical analyses of structural dynamics[END_REF]. Solutions and phase portraits of the reference solution and the one obtained by the time perturbation method (ANM) are given in Fig. (2.30) and Fig. (2.31) (the same result is obtained by BPL and IFS). It is seen that good agreement is obtained. However, results obtained with classical discrete-time step ∆t show the inaccuracy of the numerical solution. In conservative systems, the exact displacement-velocity portrait should be a closed cycle which is not satisfied by the classical integration schemes Runge-Kutta of third (RK3) and fourth-order (RK4) as can be deduced from Fig. (2.31). Note that similar results are obtained by using the following classical schemes: Newmark, the Houbolt method, α-method, the average acceleration method, the central difference method, and the Wilson θ-method. This has been confirmed by the deviation of the total energy obtained by these schemes from its initial value in [START_REF] Liu | An assessment of the differential quadrature time integration scheme for nonlinear dynamic equations[END_REF].

Application to the three degrees of freedom (3 DOF) example: Lorenz system

Goal of this example: A three-degree-of-freedom nonlinear ODE problem (Lorenz system) is used to verify the performance of the proposed methods for multi-degree-offreedom problems. The ability of time perturbation techniques to follow the chaotic and non chaotic behavior of the Lorenz system is examined One of the most famous chaotic systems was developed by Edward Lorenz [START_REF] Lorenz | Deterministic Nonperiodic Flow[END_REF] who was interested in the fluid flow patterns of the Earth's atmosphere and noticed an unexpected chaotic behavior. Depending on the choice of the parameters that exist in the equations, the Lorenz system can exhibit chaotic and nonchaotic behavior. The Lorenz equations govern, in a lower order, the dynamics of convection in a heated fluid layer and present particular challenges due to its high sensitivity to small variations in initial conditions, from continuous convection to low turbulence (chaos). The Lorenz equations are:

ẋ = σ(y -x) (2.36) ẏ = x(ρ -z) -y (2.37) ż = xy -βz (2.38)
where x, y, and z are dynamical variables, with the following initials conditions

         x(0) = x 0 y(0) = y 0 z(0) = z 0 (2.39)
The first component of the solution x is related to the convection in the atmospheric flow, while the other two components y and z are related to horizontal and vertical temperature variation. The parameters σ, ρ, and β are the related constants. Note that this system of differential equations is not linear and it is very complicated to obtain exact solutions with the presence of two nonlinear terms.

Bifurcation studies show that with the following parameters σ = 10 and β = 8/3, the critical chaos parameter is ρ = ρ cr = 27.74 [START_REF] Belokolos | Chaos in a generalized Lorenz system[END_REF]. The behavior of the solution is non-chaotic when ρ < ρ cr and chaotic otherwise. Solutions are completely dependent on initial parameters and conditions and it is very difficult to predict their behavior. For some parameter values, the orbit of the solution in the three-dimensional space is a strange attractor. For other values of the parameters, the solution tends either to converge towards a fixed point or diverge to infinity or follow a periodic oscillation motion.

When resolving this system numerically, the numerical methods provide the solutions only at a discrete point in time. To ensure a good convergence of the solution, these methods require a very small discrete-time step. However, even by reducing the integration step ∆t, the result can not be improved. This is because the integration error presents an extremum as a function of ∆t as mentioned in [START_REF] Lozi | A new accurate numerical method of approximation of chaotic solutions of dynamical model equations with quadratic nonlinearities[END_REF]. The importance of handling the Lorenz system with a time perturbation study is outlined also in [START_REF] Lozi | A new accurate numerical method of approximation of chaotic solutions of dynamical model equations with quadratic nonlinearities[END_REF].

Time perturbation series recurrence formula

Series terms of time perturbation method for this system are given by the following recurrence formulas, ∀k = 1, • • • , N we have:

x k = 1 k σ (y k-1 -x k-1 )
(2.40)

y k = 1 k ρx k-1 - k-1 i=0 x i z k-1-i -y k-1 (2.41) z k = 1 k k-1 i=0 x i y k-1-i -βz k-1 (2.42)

Numerical results

For the numerical results, the following initials conditions are set first to x 0 = 5, y 0 = 5, z 0 = 400 and the constant ρ is set to ρ = 350. .

Application to the nonlinear single pendulum problem

It is well-known that time integration methods can perform very different results in nonlinear cases. For this reason, the temporal perturbation and numerical resummation methods considered as time integration schemes are applied on the nonlinear oscillating pendulum.

Goal of this example:

Since different information about the motions of the pendulum can be precisely evaluated using given initial conditions, it is well recognized that this problem is an appropriate example for testing the performance of temporal integration schemes.

The equation of motion modeling the free, undamped simple pendulum is:

θ + ω 0 2 sin(θ) = 0 (2.43)
where θ is the angle between the rigid rod and the vertical line at time t, and ω 0 is defined as

ω 0 = g l (2.44)
where l denotes the length of the pendulum and g the acceleration due to the gravity. Note that equation (2.43) is nonlinear, due to the presence of sin(θ). The oscillations of the pendulum are subjected to the initial conditions:

   θ(0) = θ 0 θ(0) = θ0 (2.45)
where θ 0 denotes the amplitude of oscillations.

The energy in the system is conservative and can be expressed as:

1 2 θ2 -ω 2 cos(θ) = constant (2.

Time perturbation series recurrence formula

When we write θ(t) as time power series, replace the obtained expansions in the initial equation Eq.(2.43) and identify according to the power of t, we notice that it is hard to find the recurrence formula that allows to calculate the series term θ k , ∀k = 1, • • • , N . For this reason, let the following change of variables

U = cos(θ) (2.47) V = sin(θ) (2.48) F = ω 0 2 sin(θ) (2.49)
that lead to the following recurrence formulas of series terms. Therefore, ∀k = 1, • • • , N we have:

U k = 1 k k-1 j=0 V j F k-1-j (2.50) V k = 1 k i-1 j=0 U j F k-1-j (2.51) F k = - w 2 0 k U k-1 (2.52) θ k = 1 k F k-1 (2.53) (2.54)

Numerical results

For the numerical test, we consider g = 1 and l = 1. We will present two cases of initials conditions to illustrate two different motions of the single pendulum as presented in [START_REF] Kim | An improved explicit time integration method for linear and nonlinear structural dynamics[END_REF] (see Fig. (2.34)).

First choice of initials conditions

Firstly, we consider θ(0) = 0.0 and θ(0) = 1.999999238456499 as given in [START_REF] Kim | An improved explicit time integration method for linear and nonlinear structural dynamics[END_REF] to synthesize a highly nonlinear case, where the pendulum performs an oscillating motion. The period T is 33.72102056485366 [START_REF] Kim | An improved explicit time integration method for linear and nonlinear structural dynamics[END_REF]. For further information about the oscillating pendulum, the nonlinear period, and the maximum angle, the reader is invited to consult the following references [START_REF] Kim | An improved explicit time integration method for linear and nonlinear structural dynamics[END_REF], [START_REF] Fung | Solving initial value problems by differential quadrature method-part 2: second-and higher-order equations[END_REF].

If θ 0 = 0 t = 2 |Ω 0 | Ei (sin(θ/2), κ) (2.55)
where κ = 2ω/Ω 0 and Ei(z, κ) is the elliptical integral of the first kind and Ω 0 is the initial angular velocity.

As can be appreciated in Fig. (2.35), BPL solution is superposed to the reference one, while classic schemes present very large error for different step size ∆t. Even with small time step, the classical schemes RK are not able to provide a valid solution, whereas with a residue tolerance = 10 -4 , BPL reproduces the exact solution very easily. The same result is obtained with ANM and IFS. The reference solution is in black [START_REF] Kim | An improved explicit time integration method for linear and nonlinear structural dynamics[END_REF], [START_REF] Fung | Solving initial value problems by differential quadrature method-part 2: second-and higher-order equations[END_REF].

Second choice of initials conditions

Secondly, we consider θ(0) = 0.0 and θ(0) = 2.000000761543501 as given in [START_REF] Kim | Higher-order explicit time integration methods for numerical analyses of structural dynamics[END_REF] to synthesize a highly nonlinear case where the pendulum performs a complete rotation instead of oscillation. The period is T = 33.7210. Also for the second set of initial conditions, BPL is superposed to the exact solution as shown in Fig. (2.36). However, the third-order RK (RK3) method shows a huge period error and provides a completely inaccurate solution. As can be seen in Fig. (2.36), the pendulum with RK3 oscillates between the two highest points instead of performing full rotations. RK4 method with ∆t = T /400 shows a noticeable period error while BPL is in very good agreement with the exact solution.

Conclusion

Nonlinear differential equations have been successfully solved using the time perturbation and summation methods as time integrators. The first and second-order derivatives that we have developed explicitly for BPL and IFS are validated. The properties of these schemes have been well illustrated, in particular for the IFS scheme as it has been applied as a first prospective study on simple examples.

The recent resummation algorithm Meijer-G (MG) has been applied on the popular Euler equation and the simple harmonic oscillator. The proposal of studying MG in temporal cases has been validated. The conclusions that can be drawn about Meijer-G so far from the two examples are: In the case of divergent series, we noticed that Meijer-G provides an approximation of the solution using a truncation order less than BPL, ANM, and IFS. It also overcomes the pole of Padé in BPL. The continuation of Meijer-G that allows to obtain a complete continuous solution in time and the vector aspect of Meijer-G will be handled in the next chapter. Numerical examples have shown also that the time perturbation and resummation methods (ANM, BPL, and IFS) are well suited to simulation over long times (free Van der Pol oscillator study for a final time t = 10000s, nonlinear combustion equation for a final time t = 20000s). As we have seen, classical schemes often introduce a small time lag at each period, whereas time perturbation and summation methods reproduce the limit cycle more faithfully (Elastic hardening spring and nonlinear simple pendulum problems). Moreover, these methods are also appropriate for stiff problems and can deal with cases that exhibit very fast dynamics in the system (Van der Pol with µ = 1000, nonlinear combustion equation). The efficiency of these methods has also been validated for a 3 DOF example (Lorenz system). The application to models with a large number of degrees of freedom, more accurate computational efforts, and comparisons of the proposed time perturbation integrators will be studied in the next chapters. These encouraging results lead us to apply these schemes to partial differential equations (PDEs) in Chapters 3 and 4, in particular those solving fluid problems.

98

Étude numérique de méthodes d'intégration temporelle basées sur des techniques de perturbation. : application à des problèmes de mécanique Cynthia Tayeh 2022 Chapter 3

TOWARDS APPLICATIONS TO THE NAVIER-STOKES EQUATIONS

After the validation of the time perturbation and resummation algorithms on ODEs in Chapter 2, we will discuss interesting examples of time-dependent PDEs using these methods towards applying these approaches to the Navier-Stokes equation in this chapter. We consider first the heat equation and then the Burgers equation that are more suitable forms of the Navier-Stokes equations having the same convective or diffusion terms as the incompressible Navier-Stokes equations. We recall that validating the time perturbation and resummation methods on the Navier-Stokes equation, is a study that we are interested in.

We note that Hamdouni and Razafindralandy have already performed a first study on the heat equation and the Burgers equation by BPL in one dimension in [START_REF] Razafindralandy | Time integration algorithm based on divergent series resummation, for ordinary and partial differential equations[END_REF]. The extension to the two-dimensional case was treated by Deeb in [START_REF] Deeb | Intégrateurs temporels basés sur la resommation des séries divergentes[END_REF]. Concerning the application of the time perturbation resummation methods to the Navier-Stokes equations, a reduced model of the Navier-Stokes equation by the Proper Orthogonal Decomposition (POD) method is presented in [START_REF] Razafindralandy | Numerical divergent series resummation in fluid flow simulation[END_REF] using BPL. In addition, a simple example of the complete Navier-Stokes equation (without applying the reduction methods), for which the exact solution is known, has been treated in [START_REF] Deeb | Intégrateurs temporels basés sur la resommation des séries divergentes[END_REF] by BPL. Another non-reduced Navier-stokes problem is presented in [START_REF] Deeb | Borel-Laplace summation method used as time integration scheme[END_REF] by BPL. In this chapter, we will first review some existing applications of time perturbation methods on some examples in the literature by developing some points, and then discuss new ones. We further note that the IFS used as a time integration scheme has not been applied to PDEs in the previously cited work. For this reason, its properties on PDEs (efficiency in terms of feasibility and computation time) will also be discussed in this chapter. The first goal of this chapter is to advance the use of temporal perturbation methods on models represented by PDEs. These methods have not been widely applied to an extensive number of PDEs. The second goal is to see how these different numerical methods perform for models evolving a large number of degrees of freedom. Therefore, we need to determine if they are both efficient over a large time interval and with an acceptable computational cost. Therefore, the efficiency of time perturbation-resummation methods for time-dependent PDEs is not the only objective. A comparison of CPU times between classical schemes and time perturbation-resummation methods is essential. Therefore, we will present time studies for different parameters of the probem. The third goal is to show the possibility of associating these temporal perturbation methods to different spatial discretization methods and then determine whether there are any issues with this association. In this chapter, we will use the two well-known spatial discretization methods which are the finite element methods and the MAC method [START_REF] Harlow | Numerical Calculation of Time-Dependent Viscous Incompressible Flow of Fluid with Free Surface[END_REF] (a popular finite difference method). The last goal is to explore the weaknesses and strengths of each of these methods applied on PDEs. The Meijer-G approximant will also be tested on examples programmed in Python language, where access to the Meijer-G function is available. The advantages and difficulties of such a use will be presented.

For all the test problems treated in this chapter, we note that if the spatial discretization method used is the finite element method, we will use a popular open-source platform called FEniCS. The latter was founded in 2003 for solving PDEs by the finite element method and is currently being developed in collaboration with researchers from various universities and research institutes. For more information about FEniCS, as well as the most recent updates to the FEniCS software and tutorials, the interested reader is invited to visit the FEniCS web page [START_REF]Fenics project[END_REF]. Note that FEniCS can be programmed in C++ and Python. However, in this chapter, only Python is used, to be able to test Meijer-G for some examples (since the Meijer-G function is available on Python). Before proceeding to the numerical applications, it should be noted that to the best of our knowledge, there is currently no unique result concerning the summability of the formal series obtained for all the PDEs (conditions to apply BPL and IFS algorithms). That is, it is not yet known whether any set of formal solutions of a PDE is a Gevrey set and, if so, the Gevrey index is also unknown [START_REF] By | Formal Power Series Solutions of Nonlinear First Order Partial Differential Equations[END_REF]- [START_REF] Miyake | Newton polygons and gevrey indices for linear partial differential operators[END_REF]. On the other hand, the formal solution of the equations that we have already mentioned (heat equation, Burgers equation, and Navier-Stokes equation) is known to be summable in most cases. The interested reader is invited to consult the following references [START_REF] Lutz | On the Borel Summability of Divergent Solutions of the Heat Equation[END_REF]- [START_REF] Michalik | Summability of formal solutions to the n-dimensional inhomogeneous heat equation[END_REF], [START_REF] Balser | Divergent solutions of the heat equation: On an article of Lutz, Miyake and Schäfke[END_REF], [START_REF] Costin | Divergent expansion, Borel summability and three-dimensional Navier-Stokes equation[END_REF]. We recall that the Meijer-G approximant can be applied to the series that are not Borel summable [START_REF] Mera | Fast summation of divergent series and resurgent transseries from Meijer-G approximants[END_REF].

Application to the heat equation 1D

Consider the following time-dependent heat equation in a one dimensional domain Ω:

∂u ∂t = ν ∂ 2 u ∂x 2 (3.1)
with the following initial condition

u(x, t 0 ) = u 0 (x) (3.2)
associated to the Dirichlet boundary conditions. The parameter ν denotes the thermal diffusivity coefficient. Note that this coefficient is often not included in mathematical studies of the heat equation, although its value can be very important in engineering sciences.

Time perturbation series recurrence formula

We note that the temporal perturbation method is performed first and then followed by the spatial discretization for each of the series terms. The unknown of the problem u is written as time power series as

u(x, t) = u 0 + N k=1 u k (x) tk (3.3)
where u 0 is the initial known solution. Series terms of time perturbation method for Eq.(3.1) are given by the following recurrence formula, ∀k = 1, • • • , N we have:

(k + 1)u k+1 = ν d 2 u k dx 2 (3.4)
that is

u k = ν k u (2k) 0 k! (3.5) such that u (2k) 0
is the (2k)th derivative of u 0 with respect to x.

Remark 3.1.1 In general cases, the formal series of Eq.(3.1) is divergent but it is convergent if u 0 is an integer function growing at most exponential of order 2 [START_REF] Lutz | On the Borel Summability of Divergent Solutions of the Heat Equation[END_REF], i.e if:

|u 0 (x)| ≤ Ce A|x| 2 , x ∈ C (3.6)

Spatial discretization

We will use the central finite difference for the spatial discretization. Therefore, Ω is discretized into an uniform grid of size denoted by n x , and the recurrence formula for obtaining series terms (Eq.(3.4)) becomes after discretization, as follows:

(k + 1)u i k+1 = ν u i-1 k -2u i k + u i+1 k ∆x 2 , i = 1, • • • , n x and k = 1, • • • , N (3.7)
where u i are the values at the nodes.

As x belongs to Ω and the latter is discretized into a uniform grid, a vector noted by x is obtained. Therefore u(x, t) is also a vector.

Numerical example

For the numerical example, we take ν = 1, n x = 16, and Ω = [0, π]. The exact solution is [START_REF] Razafindralandy | Time integration algorithm based on divergent series resummation, for ordinary and partial differential equations[END_REF]:

u(x, t) = e -t sin x (3.8)
Therefore, the initial condition is u 0 (x) = sin x (3.9)

Goals of this example:

This example has different purposes. Firstly, it is a first application of Meijer-G to a time-dependent PDE. Secondly, it is also a first validation of Meijer-G approximant to a vector series, the notion that we have proposed in Chapter 1, based on the Gram-Schmidt procedure (see subsection 1.4.4 of Chapter 1). Thirdly, it consists in associating Meijer-G approximant to a spatial discretization technique (for example here is the finite difference method) also for the first time in the literature to the best of our knowledge. Finally, it allows validating the notion of the continuation technique for Meijer-G to obtain the complete solution, instead of only a part of the solution (the classical aspect in previous articles employing Meijer-G [START_REF] Mera | Fast summation of divergent series and resurgent transseries from Meijer-G approximants[END_REF]). To satisfy these objectives and since this test case has been realized by BPL in [START_REF] Razafindralandy | Time integration algorithm based on divergent series resummation, for ordinary and partial differential equations[END_REF], we will just present here the result by Meijer-G. Therefore, we have to check first if the new version of the Meijer-G approximant for a vector series is efficient in terms of feasibility, then we will present some comparisons.

Numerically:

The series is truncated up to a truncation order N = 5. The Gram-Schmidt procedure is used to evaluate Meijer-G for a series of vectors. First, the solution is evaluated without applying any continuation. As can be appreciated in Fig. One can deduce from this part of study that the first use of Meijer-G for timedependent PDE associated with the finite central difference method, is well validated. The use of Meijer-G for a series of vectors using the Gram-Schmidt procedure, as introduced in the subsection 1.4.4 in Chapter 1, is also well validated.

After the approval of Meijer-G on this example, we proceed now to a simple comparison between BPL and Meijer-G for N = 5, with continuation based on the residual criterion to reach a final time of t = 3s. With a residual tolerance fixed to = 3 • 10 -2 , BPL needs 82 steps and Meijer-G needs 51 steps. No CPU time studies are presented for this example since it is not the main objective and since we are still studying the efficiency of Meijer-G summation for temporal problems.

The conclusions that can be drawn from this example are the following. To the best of our knowledge, for the first time in the literature, we have a:

• Validated application of Meijer-G for time-dependent problems.

• Robust application of Meijer-G for a PDE.

• Association of Meijer-G to a spatial discretization (here is the finite difference method).

• Validated application of Meijer-G for a vector series using the Gram-Schmidt procedure.

• Efficient continuation of Meijer-G which allows obtaining a complete solution instead of only a part of the solution.

• Validated residual and first derivative calculation of Meijer-G (see Appendix C).

Application to the heat equation 2D

Now, we consider the heat equation in a two dimensional domain Ω with a source term f in the following form:

∂u ∂t = ν∆u + f (3.10)
with Dirichlet boundary conditions

u |∂Ω = g (3.11)
and with the following initial condition u(x, y, t 0 ) = u 0 (x, y) (3.12)

where ∆ is the Laplace operator. The unknown u is now a function of two spatial variables

x and y defined over a two-dimensional domain Ω.

Classical numerical resolution

For the classical resolution of this problem, the RK1 scheme is used for the temporal discretization and the finite element method for the spatial discretization using FEniCS. Therefore, first and foremost, the details of the variational formulation will be given only for this example.

Variational formulation:

The basic point of the finite element method is to express a PDE in a variational form [START_REF] Mejdi | Eléments finis pour les fluides incompressibles[END_REF]. To convert a PDE to a variational problem, three steps are required. First, we multiply the equation by a test function v, then we integrate both members of the resulting equation on the domain Ω, and finally, we perform integration by parts for the terms with second-order derivatives. The function v that multiplies the PDE is called the test function. The unknown function u that should be approximated is called the trial function. The trial and test functions belong to a so-called function space denoted by V , that defines the properties of the functions. In what follows, V is the well-known Sobolev function H 1 0 (Ω) and the definition of the integral related to the scalar product, is used. The variational problem of Eq.(3.10) is as follows: find u ∈ V such that

( ∂u ∂t , v) = (ν∆u, v) + (f, v), ∀v ∈ V (3.13)
We note that the first derivative ∂u ∂t is approximated by the classical time integration scheme RK1. After some mathematical operations by using the Gauss theorem and by considering that the test function v vanishes on the Dirichlet boundary, Eq.(3.13) becomes

( ∂u ∂t , v) + ν(∇u, ∇v) = (f, v), ∀v ∈ V (3.14)
This equation is the variational form of the original problem. It is a continuous problem that defines the solution u in the infinite dimensional function space V . The finite element method finds an approximate solution of the variational problem Eq.(3.14) by replacing the infinite space of functions with discrete or called finite-dimensional spaces V h ⊂ V where h denotes the characteristic size of an element of the mesh. Therefore, the discrete variational problem is as follows: find

u h ∈ V h ⊂ V such that ( ∂u h ∂t , v h ) + ν(∇u h , ∇v h ) = (f h , v h ), ∀v h ∈ V h . (3.15)
Following various steps of assemblage of elementary matrices, these equations are transformed into a matrix equation. Therefore, one can obtain the following semi-discrete equation in the following form:

M du h dt = Ku h + f h in Ω h (3.16)
where M is the assembled mass matrix, K is the rigidity matrix and is associated with the operator ∆, whereas u h is the vector containing the values of the nodes, and f h is the discretized source term.

Details about the mesh: For this example, we use a triangular mesh of the 2D rectangle (-2, 2) × (2, 2) with n x cells in each direction. Therefore the total number of triangles used is 2 × n x × n x and the total number of vertices is (n x + 1) × (n x + 1), (and n x will be fixed in the following).

Concerning the type of elements, the standard P1 linear Lagrange element, which is a triangle with nodes at the three vertices of degree 1, is used.

Time perturbation series recurrence formula

After presenting the classical resolution of the initial problem in the previous section using RK1 for the temporal discretization and the finite element method for the spatial discretization, we now write the initial problem in the form of a temporal perturbation and then discretize the recurrence formula that allows obtaining the series terms using the finite element method. The unknown of the problem u is written as time power series as u(x, y, t) = u 0 (x, y)

+ N k=1 u k (x, y) tn (3.17)
where u 0 is the known initial data. The time-dependent source term f and g (if it is time-dependent) are also written as time power series as

f (x, y, t) = f 0 (x, y) + N n=1 f n (x, y) tn (3.18) g(x, y, t) = g 0 (x, y) + N n=1 g n (x, y) tn (3.19)
where f 0 and g 0 are respectively the initial source term and the initial boundary data.

Continuous formulation in time:

The following equations Eq.(3.17), Eq.(3.18), and Eq. (3.19) are injected in the initial problem Eq.(3.10). By identifying the terms according to the power of t, series terms of time perturbation method for Eq.(3.10) are given by the following recurrence formula, ∀k = 1, • • • , N we have:

ku k + ν∆u k-1 = f k-1 (3.20)
Semi-discrete formulation in Ω: For the spatial discretization by the finite element method, we multiply Eq.(3.20) by the test function v which yields to the following equation:

k(u k , v) + ν(∇u k-1 , ∇v) = (f k-1 , v), ∀v ∈ V (3.21)
The approximation of the series terms at the nodes of the mesh denoted by u h k is given by the following semi-discrete formulation:

u h k = 1 k M -1 (Ku h k-1 + f h k-1 ) in Ω h (3.22)
where M is the mass matrix, and K is the stiffness matrix. We note that only one M-matrix inversion is required for the whole simulation.

Numerical example

For the numerical example, an interesting problem called the "Diffusion of a Gaussian function" is chosen. The initial data is as follows u 0 (x, y) = e -ax 2 -ay 2 (3.23)

Homogeneous Dirichlet boundary conditions are used (g = 0). We take a = 5 and the domain Ω = [-2, 2] × [2, 2] on the following grid n x × n x with n x = 30. The number of degrees of freedom denoted by ndofs, corresponding to this grid is ndofs = 961. The series terms are first evaluated and the Euclidean norm of each term of the series at the zeroth step (without continuation procedure) is presented in Fig. (3.3). The latter reports a steady increase as truncation order increases. Note that high values of the Euclidean norm of series terms are not a classical behavior in the ANM with perturbation parameters others than time.

Towards applications to the Navier-Stokes equations We also tested the Meijer-G resummation for this test case for N = 5. However, we were unable to compute the Meijer-G sum on some points on the domain. This could be due to the small value of t max in the Meijer-G relation sum, or to the fact that the implementation of the Meijer-G function does not always give the right answers in some cases of specific parameters [104]; or the Euler Gamma function is infty for some parameters. Therefore, it would be interesting to study the cause.

Comparisons:

After testing the effectiveness of the methods in this example, we want to compare them to the classical scheme and then to each other. The final time is fixed to t = 1s. For sake of comparisons, the CPU number is fixed to 1 (as FEniCS would use as many CPUs as possible by default), the Sparse Linear Equation Solver (SuperLU) is used for the resolution [105]. The residual criterion is used for continuation, and t max is sought using dichotomy. The residual norm of the classical scheme is calculated and will be used as the prescribed residual norm tolerance of ANM, BPL, and IFS. All computations have been performed on an Intel Xeon Gold 5120 @ 3.2 GHz, using only one thread.

For the rest of the manuscript, we note by: • Total (s): The total CPU time spent in all the continuation procedures to reach the final time.

• Series (s): The CPU time spent to calculate the series term during all the continuation procedures, to reach the final time.

• Residual Series (s): The CPU time spent for evaluating the polynomial residual (dichotomy algorithm to perform the continuation technique and evaluate the range of validity via the residual criterion).

• Residual BPL (s): The CPU time spent for evaluating the BPL residual.

• Residual IFS (s): The CPU time spent for evaluating the IFS residual.

• #Residual: The number of times the residual is calculated during the continuation.

• < t max >: The average of the upper bound of the validity domain [0, t max ] to reach the final time.

• Steps: The number of the continuation steps needed to reach the final time.

The CPU time is all measured in seconds. The software used is FEniCS where the mass matrix M assembles once for all.

Remark 3.2.1 We note that since the obtained system-matrix to solve is always the same for each time step, it is not necessary to factorize the mass matrix at each increment.

It should be possible on FEniCS to store the mass matrix one time, and then use it every time, reducing "solver" to a linear algebra multiplication. This avoids assembly and LU factorization in every solve, which can save CPU time. This is done by defining a "LUSolver" object and asking for reusing the matrix factorization and then saving the factorization between time steps. The part of the code for direct solvers on FEniCS is considered here: solver = LUSolver() solver.set_operator(M) solver.parameters["reuse_factorization"] = True However, as far as we look into the documentation of FEniCS, the option has been removed. "Remove reuse_factorization. This parameter is now deprecated, there are some demo which uses it, it must be removed" as mentioned in [106]. Therefore, the CPU time for series computation in FEniCS is not optimal, since factorizing the mass matrix is performed for every truncation order and every step. But, it is possible to estimate an optimal CPU time by considering some assumptions for the computation of the time perturbation approaches and the classical scheme. Therefore, for the rest of the manuscript, we consider that the factorization time is 90% of the solving time [START_REF] Zadeh | Performance evaluation of SuperLU and PAR-ADISO in power system load flow calculations[END_REF] The "CPU Series(s)" are computed in the same way.

Classical scheme:

The classical scheme RK1 with a fixed discrete-time ∆t = 10 -3 s needs 1000 steps to reach the final time t = 1s with a total CPU time equal to 110.05s. The residual norm of RK1 is ||Res|| 2 ≈ 2.2 • 10 -3 at the final time t = 1s. It was not possible to get a better quality of the solutions with RK1. This time step is then the only one where it is possible to perform comparisons. We will use this value of the residual norm, as the residual tolerance , to perform comparisons of perturbation methods. Therefore, for the following comparisons, is fixed to = 2.2 • 10 -3 .

ANM:

For time perturbation method ANM, we present in Tab.(3.1) some statistics concerning the CPU times spent to reach the final time (Total (s)), the time dedicated to calculate the series terms (Series (s)) during all the continuation techniques and the residual computation time (Residual Series(s)) (dichotomy loop-algorithm). The number of steps required to reach the final time and the average of t max are also presented in this table. As already mentioned, the residual tolerance is fixed to = 2.2 • 10 -3 and the study is carried out with different truncation orders N . It can be observed that the CPU time for evaluating the series terms is between 13s and 23s for different truncation orders N , varying from N = 5 to N = 25. The total CPU time taken by the ANM is approximately 25s for different truncation orders, while the one taken by the classical scheme is 110.05s for the same residual norm (that is fixed as the residual tolerance for ANM, BPL, and IFS). Therefore, the ANM method is about 5 times faster than the classical scheme for this example. It can be observed also that the CPU time seems to decrease and the t max increase by increasing the truncation orders.

N Steps

BPL:

The same statistic is presented for BPL in Tab. (3.2). For BPL, one can notice that it is fast in terms of CPU time for both the evaluation of the series terms and the calculation of t max using the residual criterion. The t max values are higher than those of the previous table (Tab.(3.1)). We also notice that for N = 20, BPL seems to be faster and more optimal than the other orders, requiring only a total CPU time of 3.68s to reach the final time. and a viscosity ν = 1. The grid used is 30 × 30.

N Steps

Comparison between RK1, ANM, BPL, and IFS:

Finally, a brief summary of the methods comparisons is derived from the tables Tab. ), one can conclude that for the same fixed residual tolerance and same truncation order, BPL is the fastest in CPU time, requires the fewest steps, and thus has the largest range of validity to reach the final fixed time. We have seen that increasing the truncation orders increases the CPU time of IFS significantly, but decreases the CPU time of ANM and BPL for some orders. Finally, ANM and BPL are much faster than the classical scheme for this fixed residual tolerance = 2.2 • 10 -3 , but the latter is faster than IFS in terms of computing time. We would like to point out that the three methods are efficient in terms of feasibility. Total (s) RK1 Total (s) ANM Total (s) BPL Total (s) IFS 110.05 14.28 3.68 729.27

Table 3.4 -CPU times required to reach a final time t = 1s for the heat equation 2D, using RK1 and ANM, BPL, IFS with N = 20 for a fixed residual tolerance = 2.2 • 10 -3 and a viscosity ν = 1. The grid used is 30 × 30.

Study for different mesh size

We are now interested in investigating the sensitivity of these methods to the mesh size. We perform the study for a fixed truncation order of N = 20 and three different mesh size n x × n x . The residuals of the classical schemes are calculated at the final time. These residual norms were used as residual tolerance for the temporal perturbation resummation methods computation. The following grids are considered:

• For n x = 30: The mesh size is 30 × 30. The number of degrees of freedom denoted by ndofs, correspoding to this grid is ndofs = 961. The residual tolerance set here is = 2.2 • 10 -3 .

• For n x = 60: The mesh size is 60 × 60. The number of degrees of freedom is ndofs= 3721. The residual tolerance set here is = 1.9 • 10 -3 .

• For n x = 120: The mesh size is 120 × 120. The number of degrees of freedom is ndofs = 14641. The residual tolerance set here is = 3.7 • 10 -2 .

The evolution of the CPU time and the t max with the grid size, are presented respectively in Tab.(3.5) and Fig. (3.8). By increasing the grid size, one can see that the total CPU time of time perturbation-resummation has increased significantly. For this example and this residual tolerance, one can see that RK1 has a lower sensitivity to the grid size than the others approaches. For n x = 60, we notice that BPL is always the fastest among RK1, ANM, and IFS. However for n x = 120, RK1 becomes the fastest. We will see if this behavior remains valid for all the following test cases of this chapter. Finally, increasing the grid size from 30×30 (n x = 30) to 120×120 (n x = 120) decreases the range of validity t max of the three methods, as can be seen in Fig. (3.8).

n

Application to the Burgers equation

The two-dimensional coupled viscous Burgers equation is a more suitable form of the Navier-Stokes equations having the same convective and diffusion form as the incompressible Navier-Stokes equations. It is a simple model for understanding various flows and physical problems such as hydrodynamic turbulence, shock wave theory, wave processes in thermoelastic media, and vorticity trajectories. Consider the two dimensional Burgers equations of the following form:

∂u ∂t + u • ∇u = ν∆u (3.25)
with Dirichlet boundary conditions

u |∂Ω = g (3.26)
and with the following initials conditions u(x, t 0 ) = u 0 (x) = u 0 (x, y) (3.27)

Note that the velocity u = (u, v) where u and v are the components of u in the x-and y-directions respectively and x = (x, y). The function g is a given Dirichlet boundary data that can be time-dependent.

Classical numerical resolution

We also use the finite element method for the spatial discretization using FEniCS and RK1 for the temporal discretization. After different operations, one can finally obtain the following semi-discrete equation of Eq.(3.25) in the following form:

M du h dt + N(u h , u h ) = νKu h , in Ω h (3.28)
where M is the assembled mass matrix, K is the rigidity matrix associated with the operator ∆, N is the nonlinear discretized operator corresponding to the convective term, whereas u h is the vector containing the values of the nodes.

Details about the mesh:

We use a triangular mesh of the 2D rectangle (0, 0.5)×(0, 0.5) with n x cells in each direction (n x will be chosen later). Therefore the total number of triangles used is 2 × n x × n x and the total number of vertices is (n x + 1) × (n x + 1) with a crossed diagonal directions. Concerning the type of elements, the standard P1 linear Lagrange element is used.

Time perturbation series recurrence formula

Having presented the classical solution of the initial problem in the previous section using RK1 for the temporal discretization and the finite element method for the spatial discretization, we now write the initial problem in the form of a temporal perturbation and then discretize the recurrence formula that permits to obtain the series terms using the finite element method. The unknown of the problem u and g (if it is time-dependent) are written as time power series u(x, y, t) = u 0 (x, y)

+ N k=1 u k (x, y) tk (3.29) g(x, y, t) = g 0 (x, y) + N k=1 g k (x, y) tk (3.30)
where u 0 and g 0 are respectively the initial solution and the initial boundary data.

Remark 3.3.1

In the one-dimensional case, Lysik showed that the formal solution of the Burgers equation is generally divergent, and is Gevrey of index 1 [START_REF] Lysik | Borel summable solutions of the Burgers equation[END_REF]. In this section, the two-dimensional case is studied. One can assume that the intervening series diverge no faster than a Gevrey series of index 1 as in [START_REF] Deeb | Intégrateurs temporels basés sur la resommation des séries divergentes[END_REF] and that we can apply the resummation to this example.

Continuous formulation in time:

Injecting Eq.(3.29) and Eq.(3.30) in Eq.(3.25) and identifying the terms according to the power of t, one can obtain the following recurrence formula that permits to evaluate the series terms, ∀k = 1, • • • , N :

(k + 1)u k+1 + k n=0 u k-n ∇u n = ν∆u k (3.31)
Semi-discrete formulation in Ω: For the spatial discretization by the finite element method, we multiply Eq.(3.31) by the test function v which yields to the following equation:

(k + 1)(u k+1 , v) + ( k n=0 u k-n ∇u n , v) = ν(∆u k , v), ∀v ∈ V (3.32)
After some operations, the approximation of the series terms at the nodes of the mesh denoted by u h k+1 is given by the following semi-discrete formulation:

(k + 1)Mu h k+1 + k n=0 N(u h k-n , u h n ) = νKu h k , in Ω h (3.33)
where M is the mass matrix, K is the stiffness matrix, and N is the nonlinear discretized operator corresponding to the convective term. We note that only one M-matrix inversion is required for the whole simulation (as time perturbation-resummation methods are explicit schemes).

Numerical example

The computational two-dimensional domain is Ω = [0 × 0.5] × [0 × 0.5] with the initial conditions

   u(x, y, 0) = sin(πx) + cos(πy) v(x, y, 0) = x + y (3.34)
and the following boundary conditions:

∀ t ≥ 0 and 0 ≤ y ≤ 0.5, we have:

              
u(0, y, t) = cos(πy) u(0.5, y, t) = 1 + cos(πy) v(0, y, t) = y v(0.5, y, t) = 0.5 + y

(3.35)
∀ t ≥ 0 and 0 ≤ x ≤ 0.5, we have:

               u(x, 0, t) = 1 + sin(πx) u(x, 0.5, t) = sin(πx) v(x, 0, t) = x v(x, 0.5, t) = x + 0.5 (3.36)
The Euclidean norm of each term of the series at zeroth step (without continuation procedure) is plotted in Fig. (3.9). The same behavior of the norms of the series terms as the previous case (heat equation) is depicted here as well. However, the norm of the series terms are much larger than that of the heat equation.

Efficiency and feasibility study:

This example has been discussed in several papers by many authors, one may cite Jain and Holla [START_REF] Jain | Numerical solution of coupled Burgers' equations[END_REF], Bahadir [START_REF] Bahadır | A fully implicit finite-difference scheme for two-dimensional Burgers' equations[END_REF], Srivastava et al. [START_REF] Srivastava | An implicit logarithmic finitedifference technique for two dimensional coupled viscous Burgers' equation[END_REF], [START_REF] Srivastava | Numerical solutions of coupled Burgers' equations by an implicit finite-difference scheme[END_REF], and Tamsir et al. [START_REF] Tamsir | Numerical Solution of two dimensional coupled viscous Burgers Equation using the Modified Cubic B Spline Differential Quadrature Method[END_REF]. We will compare the effectiveness of the feasibility of the current work with the result of the previously cited work. For this reason, we use the same parameters as those presented in these papers. The grid is fixed to 20 × 20 (n x = 20), the viscosity to ν = 0.02, and the final time to t = 0.625s as in [START_REF] Jain | Numerical solution of coupled Burgers' equations[END_REF]- [START_REF] Tamsir | Numerical Solution of two dimensional coupled viscous Burgers Equation using the Modified Cubic B Spline Differential Quadrature Method[END_REF]. Tab. (3.6) and Tab. (3.7) show the comparisons of the velocity numerical results at some typical mesh points obtained using time perturbation and resummations schemes at t = 0.625s with the methods of Jain and Holla [START_REF] Jain | Numerical solution of coupled Burgers' equations[END_REF] (based on cubic spline function technique), Bahadir [START_REF] Bahadır | A fully implicit finite-difference scheme for two-dimensional Burgers' equations[END_REF], Srivastava et al. [START_REF] Srivastava | An implicit logarithmic finitedifference technique for two dimensional coupled viscous Burgers' equation[END_REF], [START_REF] Srivastava | Numerical solutions of coupled Burgers' equations by an implicit finite-difference scheme[END_REF] (denoted by I-LFDM and Expo-FDM in their works) and Tamsir et al. [START_REF] Tamsir | Numerical Solution of two dimensional coupled viscous Burgers Equation using the Modified Cubic B Spline Differential Quadrature Method[END_REF] (denoted by MCB-DQM, based on a the modified cubic B-spline differential quadrature method). Since the evaluated error between ANM, BPL and IFS is very small (approximately ≥ 10 -6 depending on the tolerance), we note by the "present work" the results of the three methods in Tab. (3.6) and Tab.(3.7). From both tables, it can be noticed that our present work achieves similar results as those of the literature. Efficient numerical results indicate that the time perturbationresummation schemes are reliable and robust techniques for the numerical solutions of the Burgers equation. The results obtained by these methods are consistent with the literature and references that have treated this example.

(x, y)

Ref. [START_REF] Srivastava | An implicit logarithmic finitedifference technique for two dimensional coupled viscous Burgers' equation[END_REF] Ref. [START_REF] Srivastava | Numerical solutions of coupled Burgers' equations by an implicit finite-difference scheme[END_REF] Ref. [START_REF] Bahadır | A fully implicit finite-difference scheme for two-dimensional Burgers' equations[END_REF] Ref. [START_REF] Jain | Numerical solution of coupled Burgers' equations[END_REF] Ref. [START_REF] Tamsir | Numerical Solution of two dimensional coupled viscous Burgers Equation using the Modified Cubic B Spline Differential Quadrature Method[END_REF]

(x, y)

Ref. [START_REF] Srivastava | An implicit logarithmic finitedifference technique for two dimensional coupled viscous Burgers' equation[END_REF] Ref. [START_REF] Srivastava | Numerical solutions of coupled Burgers' equations by an implicit finite-difference scheme[END_REF] Ref. [START_REF] Bahadır | A fully implicit finite-difference scheme for two-dimensional Burgers' equations[END_REF] Ref. [START_REF] Jain | Numerical solution of coupled Burgers' equations[END_REF] Ref. [START_REF] Tamsir | Numerical Solution of two dimensional coupled viscous Burgers Equation using the Modified Cubic B Spline Differential Quadrature Method[END_REF]

Comparisons:

After checking the efficiency of the methods for the Burgers equation, we are interested in carrying out the comparisons first between the classical methods and the time perturbation-resummation methods, and then between the time perturbation-resummation methods, for an arbitrary final time of t = 0.17s. Like the previous example, for sake of comparisons, the CPU number is fixed to 1 (as FEniCS would use as many CPUs as possible by default), the Sparse Linear Equation Solver (SuperLU) is used for the resolution [105]. The residual criterion is used for continuation, and t max is sought using dichotomy. The residual norm of the classical scheme is calculated and will be used as the prescribed residual norm tolerance of ANM, BPL, and IFS. All computations have been performed on an Intel Xeon Gold 5120 @ 3.2 GHz, using only one thread.

Classical scheme:

The classical scheme RK1 with a fixed time step ∆t = 10 -3 s needs only 170 steps to reach the final time, with a total CPU time equal to 9.53s. The residual norm is ||Res|| 2 ≈ 7.3 • 10 -3 at the final time. To have the same residual range for all methods, we will use this value of residual norm to perform comparisons of perturbation methods. Therefore, for the following comparisons, the residual tolerance is fixed to = 2.2 • 10 -3 .

ANM:

For time perturbation method ANM, we present in Tab.(3.8) some statistics concerning the CPU times spent to reach the final time (Total (s)), the time dedicated to calculate the series terms (Series (s)) during all the continuation technique, and the residual computation time (Residual Series(s)). The number of steps required to reach the final time and the average of t max are also presented in this table. As already mentioned, the residual tolerance is fixed to = 7.3 • 10 -3 and the study is carried out with different truncation orders N . It can be observed that most of the total CPU time, is spent by evaluating the series terms at every step of the continuation. We also notice that the CPU time increases by increasing the truncation order.

N Steps

BPL:

The same statistic is presented for BPL in Tab. (3.9). For BPL, one can notice that it is fast in terms of CPU time for both the evaluation of the series terms and the calculation of t max using the residual criterion.

Remark 3.3.3 We also investigated the CPU time required to compute one-IFS solution for various truncation orders N . It was found that evaluating one-IFS solution alone is

time-consuming as well as its derivative.

Comparison between RK1, ANM, BPL, and IFS:

Finally, a brief summary of the methods comparisons is derived from the tables Tab.(3.8), Tab.(3.9) and Tab.(3.10) and is depicted in Fig. (3.10). One can conclude that, for this example, BPL is the fastest among AMN and IFS. It has the highest range of validity and requires the fewest steps to reach the same final time. IFS is the most time-consuming. The classical scheme was slightly faster than BPL for this example because the calculation of the series terms for the Burgers equation is considered to be time-consuming. However, the t max of BPL is larger than that of RK1. 123

Study for different mesh sizes

We are now interested in investigating the sensitivity of these methods to the mesh size. We perform the study with a fixed truncation order of N = 20 and three different mesh size n x × n x . The residuals of the classical scheme are calculated at the final time. These residual norms are used as residual tolerances in the temporal perturbation methods computation. The following grids are considered:

• For n x = 20: The mesh size is 20 × 20. The number of degrees of freedom denoted by ndofs, corresponding to this size is ndofs = 6562. The residual tolerance set here is = 7.3 • 10 -3 .

• For n x = 40: The mesh size is 40 × 40. The number of degrees of freedom is ndofs = 25922. The residual tolerance set here is = 4.0 • 10 -3 .

• For n x = 80: The mesh size is 80 × 80. The number of degrees of freedom is ndofs = 103042. The residual tolerance set here is = 4.9 • 10 -1 .

The evolution of the CPU time and the t max with the grid size, are presented in Fig. (3.11). By increasing the mesh size, one can see that the t max decreases, and the CPU time increases, for the three methods. We note that the t max of IFS and ANM decreases much more than that of BPL when n x varies from 40 to 80.

Application to Navier-Stokes equations

After studying the time perturbation-resummation methods on the heat equation and the Burgers equation, we now proceed to their application to the Navier-Stokes equations that describe the motion of a viscous fluid flow. The Navier-Stokes equations contain the conservation laws of mass and momentum and describe the spatial and temporal change of the velocity field of the fluid. The numerical solution of the Navier-Stokes equations can be split into two main categories: compressible and incompressible flows. We have chosen to consider the assumptions of the incompressibility of the fluid in this manuscript. Therefore, consider the time-dependent incompressible Navier-Stokes equations (NSE) in a connected domain Ω ∈ R 2 with a piece-wise smooth boundary ∂Ω, in the following form:

ρ ∂u ∂t -µ∆u + ρu • ∇u + ∇p = f (3.37) ∇ • u = 0 (3.38)
with Dirichlet boundary conditions

u(x, t) = g(x, t) in ∂Ω (3.39)
with the following initial conditions

u(x, 0) = u 0 (x) (3.40)
Eq.(3.37) derives from the conservation of momentum, and Eq.(3.38) derives from the conservation of mass and is the incompressibility condition. Therefore, we note that for this problem we have two residues to consider (the momentum residual and the continuity residual linked respectively to Eq.(3.37) and Eq. (3.38)).

The unknowns of the problem are the velocity u = (u, v) and the pressure p. Note that u and v are the components of the fluid velocities in the x-and y-directions respectively. µ denotes the dynamic viscosity, f the given source term and ρ is the density of the fluid. We define the kinematic viscosity as ν = µ ρ .

A description of each term in the Navier-Stokes equations is presented in the following:

• ∂u ∂t : The derivative of the velocity with respect to time. During each time step of the simulation, this term will be calculated at all grid points containing fluid.

• ρ(u • ∇u): The convection term.

This term is derived from the conservation of momentum.

• ∇p: The pressure term. This term refers to the forces produced by pressure differences within a fluid. As the fluid is assumed to be incompressible, we will combine the pressure term with Eq. (3.38) to ensure that the flow remains incompressible.

• µ∆u: The viscosity term.

In thick fluids, friction forces cause the velocity of the fluid to move toward the neighborhood average. This term captures this relationship using the ∆ operator.

• f : The source term. This term can include all external forces, such as gravity or forces of contact with objects.

One can note that the major difficulty in finding a solution to the Navier-Stokes equations is the coupling between velocity and pressure through the incompressibility condition which is Eq. (3.38). Different methods for solving the Navier-Stokes equations are proposed in the literature. One can mention artificial compressibility methods, projection methods, spectral methods, or penalty methods such as the augmented Lagrangian method. Projection methods are the most successful ones because of their practical simplicity of implementation and flexibility of use. They have been widely used and developed. The splitting or projection method proposed by Chorin in 1968 [START_REF] Chorin | Numerical solution of the Navier-Stokes equations[END_REF] and Temam in 1969 [START_REF] Temam | Sur l'approximation de la solution des équations de Navier-Stokes[END_REF], also known as the Chorin method, is one of the oldest methods. It is also known in the literature as the "non-incremental pressure-correction scheme" or the "pressure-Poisson approach". We have chosen to use this original version of the projection algorithm since it is based on a fully explicit formulation like the time perturbation-resummation approaches. The Chorin projection method is described in the following subsection.

Chorin projection method and classical resolution

In the first place, an intermediate velocity field is determined before it is projected to zero divergences by a Hodge decomposition. Given the approximate discrete solution at time t n , we look for the solution at time t n+1 following a time step ∆t such that t n+1 = t n + ∆t. The Chorin projection method is based on the following three steps:

• Prediction: The first step consists in ignoring the pressure gradient term of the momentum equation and then calculating an intermediate velocity u * by

     ρ u * -u n ∆t -µ∆u n + ρu n • ∇u n = f u * = g in ∂Ω (3.41)
where u n is the velocity at n th time step. Note that the first derivative ∂u ∂t was approximated by RK1.

• Pressure solution: The second step consists in correcting the intermediate velocity to obtain the final time step solution.

       ∆p n+1 = ρ ∆t ∇ • u * ∂p n+1 ∂n = 0 in ∂Ω (3.42)
where n is the unit normal vector at ∂Ω.

• Correction of the velocity field: The final step consists in correcting the predicted velocity in order to satisfy the continuity equation.

     u n+1 = u * - ∆t ρ ∇p n+1 ∇ • u n+1 = 0 in ∂Ω (3.43)
After presenting the classical resolution of the problem, we are interested in solving Eq.(3.37) and Eq. (3.38) with the time perturbation-resummation methods in the next subsection.

Time perturbation series recurrence formula

We first introduce the development in temporal series in detail, in order to see how the problem will be decomposed and then the obtained equations will be discretized. The unknowns of the problem defined by Eq.(3.37) and Eq. (3.38), which are u and p are written as time power series:

u(x, t) = u 0 (x, y) + N k=1 u k (x, y) tk (3.44) p(x, t) = N k=0 p k (x, y) tk (3.45)
where u 0 is the initial known solution. We note that p 0 can be unknown in many cases and is calculated from u * 0 . The external forces f and the Dirichlet boundary data g (if it is time-dependent) are also written as time power series:

f (x, t) = f 0 (x, y) + N k=1 f k (x, y) tk (3.46) g(x, t) = g 0 (x, y) + N k=1 g k (x, y) tk (3.47)
where f 0 and g 0 are respectively the initial source term and the initial Dirichlet boundary data. Then, Eq.(3.44), Eq.(3.45), Eq.(3.46) and Eq.(3.47) are injected in the initial problem Eq.(3.37) and Eq. (3.38). Finally, by equating the terms according to the power of t, the following nonlinear systems are obtained: t0 :

         ρu 1 = µ∆u 0 -ρu 0 • ∇u 0 -∇p 0 + f 0 ∇ • u 1 = 0 u 1 = g 1 in ∂Ω (3.48) t1 :          2ρu 2 = µ∆u 1 -ρu 1 • ∇u 0 -ρu 0 • ∇u 1 -∇p 1 + f 1 ∇ • u 2 = 0 u 2 = g 2 in ∂Ω (3.49) . . . tk , k ∈ [2, N ] :              ρ(k + 1)u k+1 = µ∆u k -ρ k r=0 u r • ∇u k-r -∇p k + f k ∇ • u k+1 = 0 u k+1 = g k+1 in ∂Ω (3.50)
To resolve this system, a prediction-correction algorithm is used. Therefore, the concept of the Chorin projection method is applied to every order.

Calculation of u 1 and p 0 :

A detail of the resolution steps of the first system of t0 (Eq.(3.48)) that allows to calculate u 1 and p 0 is presented. First, consider again the first equation of the system presented in Eq. (3.48):

ρu 1 = µ∆u 0 -ρu 0 • ∇u 0 -∇p 0 + f 0 (3.51)
By applying a divergence to the latter on both sides, one can obtain:

∇ • (ρu 1 ) = ∇ • (µ∆u 0 -ρu 0 • ∇u 0 -∇p 0 + f 0 ) (3.52)
Since ∇ • u 1 = 0 (continuity equation, see Eq.(3.38)), the left hand side of Eq.(3.52) vanishes. Therefore, one can get:

∇ • (µ∆u 0 -ρu 0 • ∇u 0 -∇p 0 + f 0 ) = 0 (3.53)
Thus, the development of the gradient in Eq.(3.53) yields to

-∆p 0 + ∇ • (µ∆u 0 -ρu 0 • ∇u 0 + f 0 ) = 0 (3.54)
An auxiliary variable denoted by u * 1 is introduced [START_REF] Deeb | Intégrateurs temporels basés sur la resommation des séries divergentes[END_REF] and defined by

ρu * 1 = µ∆u 0 -ρu 0 • ∇u 0 + f 0 (3.55)
Therefore, Eq.(3.54) becomes

-∆p 0 + ∇ • (ρu * 1 ) = 0 (3.56) which is equivalent to ∆p 0 = ρ ∇ • u * 1 (3.57)
associated to the following Neumann boundary condition

∇p 0 • n = 0 (3.58)
where n denotes the normal to ∂Ω. Then, Eq.(3.51) becomes

ρu 1 = µ∆u 0 -ρu 0 • ∇u 0 -∇p 0 + f 0 (3.59) = ρu * 1 -∇p 0 (3.60)
Therefore,

u 1 = u * 1 - 1 ρ ∇p 0 (3.61)
Finally as summary, the steps that should be resolved to obtain u 1 an p 0 are as follows:

           ρu * 1 = µ∆u 0 -ρu 0 • ∇u 0 + f 0 ∆p 0 = ρ ∇ • u * 1 u 1 = u * 1 - 1 ρ ∇p 0 (3.62)
associated to the following Dirichlet boundary condition

u 1 = g 1 in ∂Ω (3.63)
and to the Neumann boundary condition

∇p 0 • n = 0 (3.64)
where n denotes the normal to ∂Ω.

Calculation of u k and p

k-1 for k ∈ [2, N ]:
The same procedure is applied for u 2 , u 3 , • • • , u N . Therefore, the steps that are resolved to obtain u k and p k-1 for k ∈ [2, N ] are as follows:

t 0 :            ρu * 1 = µ∆u 0 -ρu 0 • ∇u 0 + f 0 ∆p 0 = ρ ∇ • u * 1 u 1 = u * 1 - 1 ρ ∇p 0 (3.65) t 1 :            2ρu * 2 = µ∆u 1 -ρu 0 • ∇u 1 -ρu 1 • ∇u 0 + f 1 ∆p 1 = 2 ρ ∇ • u * 2 u 2 = u * 2 - 1 2 ρ ∇p 1 (3.66) . . . (3.67) t k , k ∈ [2, N ] :                  ρ (k + 1)u * k+1 = µ∆u k -ρ k r=0 u r • ∇u k-r +f k ∆p k = ρ (k + 1) ∇ • u * k+1 u k+1 = u * k+1 - 1 ρ (k + 1) ∇p k (3.

Spatial discretization

The obtained equations Eq.(3.68), Eq.(3.69) and Eq.(3.70) that permit to obtain series terms are discretized. All popular methods, such as the finite difference method, the finite element method, the spectral method, or the spectral element method, can be used for the Navier-Stokes equations. In this section, we have chosen to use the finite difference MAC method for the Lid-driven cavity problem, and the finite element method for the Navier-Stokes in an L-shaped domain, and the Channel flow examples using FEniCS.

Some useful definitions

Before we present the various test cases, we introduce some definitions that will be used for the rest of the document.

• Stream function:

The flow field of some problem is visualized through streamline obtained from stream function denoted by ψ. Stream function is defined from the velocity components u and v. The formulations that relate the stream function ψ to velocity components u and v for a 2D flow are given by:

u = ∂ψ ∂y and v = - ∂ψ ∂x (3.71)
Streamlines are the lines which are tangent to the velocity vector at each point of the domain: • Reynolds number:

u v = dx dy (3.
The Reynolds number characterizes a flow by classifying the latter from laminar (at low Reynolds numbers) to turbulent (at high Reynolds numbers). It is defined as:

Re = ρU L µ = U L ν (3.74)
where U and L are respectively a speed and a reference length characteristic of the application. L denotes the characteristic length of the flow. The viscosities ν and µ are defined previously. This number can range from zero to millions or even tens of millions. This wide range corresponds to very different physical situations.

Numerical examples

First, a problem with a known analytical solution called the Channel or Poiseuille flow is presented by these methods. To validate the terms of the computed series, the error between the analytical solution and the solution obtained by the time perturbation methods is evaluated. As this example is a priori simple and is performed as a first validation of the time perturbation methods on the Navier-Stokes equations, the comparison and computation time studies will not be presented for the Channel flow example. Then, an application of the methods to a flow in an L-shaped domain is presented. Finally, their application to a reference problem in fluid mechanics, which is the LID-driven cavity for different Reynolds numbers, is provided. For the last two examples, the comparison between the time perturbation and resummation methods with each other and with the classical methods is performed.

Test problem 1: Channel flow Presentation of the problem

A steady, laminar flow of a viscous fluid between two horizontal parallel plates separated by a distance, H, is defined as the Channel or Poiseuille flow. The flow is induced by a pressure gradient across the length of the plates, L, and is characterized by a symmetric 2D parabolic velocity profile about the horizontal mid-plane, as shown in Fig. (3.12). Note that there are no internal forces. This problem has a well-known analytical solution in the steady-state. First, the boundary conditions are imposed by the symmetry assumptions of the axisymmetric configuration. The pressure is fixed to p = 8 at the inlet x = 0, and p = 0 at the outlet x = 1. Then, the velocity is fixed to u = (0, 0) on the walls of the channel, that is at y = 0 and y = 1. Finally, we note that without further specifications, this defines the pressure drop and should result in unit maximum velocity at the inlet and outlet, as well as a parabolic velocity profile. Using the boundary conditions, we end up with the following Navier-Stokes equations:

- 1 ρ ∂p ∂x + µ ∂ 2 u ∂y 2 = 0 (3.75) - 1 ρ ∂p ∂y = 0 (3.76)
The Poiseuille equation gives the velocity of the fluid as a function of the radius at which it is placed in the pipe. It reflects the parabolic velocity profile through the channel, symmetric about the central axis, in the following form:

u(y) = -6U 0 y 2 H 2 - y H (3.77)
where U 0 is the characteristic inlet velocity. The pressure gradient along x is also analytical and is in the following form:

∂p ∂x = -12U 0 µ h 2 (3.78)

Discretization

The finite element method is chosen to be used as a spatial discretization for this problem. Some finite element spaces are well known to be unstable for the Navier-Stokes equations. The use of continuous piecewise polynomials of the first degree for both velocity and pressure is an example of an unstable pair of finite element spaces. Using a pair of unstable spaces usually results in a pressure solution with spurious oscillations. Therefore, we use piecewise continuous quadratic elements for viscosity and piecewise linear elements for pressure. These elements combine to form the Taylor-Hood element. We use an uniform finite element triangular mesh over the 2D unit square [0, 1] × [0, 1]. The mesh consists of cells, which in 2D are triangles with straight sides. The square is divided into n x × n x rectangles, each divided into a pair of triangles. The total number of triangles will be 2 × n x × n x and the total number of vertices will be (n x + 1) × (n x + 1). The size grid n x is fixed to 16.

Simulation, implementation and results

The viscosity, the density, and the width of the channel are fixed to 1 (ν = µ = 1 and ρ = 1). To apply time perturbation and resummation methods, we calculate the series terms already evaluated. Then, we apply the algorithms of ANM, BPL, and IFS with continuation. The truncation order is fixed to N = 6 and the tolerance of ANM is fixed to δ = 10 -4 . We calculate the L2-norm error at the nodes to check the effectiveness of the time perturbation method and verify that our implementation is correct. The latter is presented in Fig. (3.14). The error decreases over time to zero and is approximately 10 -3 at time t = 0.5s. The velocity profile at t = 5s via a probe-line is presented in Fig. (3.13). It is consistent with the analytical solution. The parabolic profile is also obtained with BPL and IFS. We do not present CPU time studies nor comparisons between methods for this test case presented just as the first validation of time perturbation resummation methods to the simple form of Navier-Stokes equations.

Test problem 2: Navier-Stokes equations on an L-shaped domain

Consider the Navier-Stokes equations governing the flow of an incompressible viscous fluid on an L-shaped domain. The latter is known as the L-shaped domain since it consists of a subset of the unit square while neglecting the upper right quadrant. This interesting shape of the domain has been used in a variety of articles and theses. One may cite [START_REF] Reyes | Stabilized reduced order models in low speed flows[END_REF]- [START_REF] Hsu | The design of L-shaped runners for gravity casting[END_REF].

Goal of this example

This example aims to study the efficiency of time perturbation and resummation methods for Navier-Stokes on an L-shaped domain geometry that has a corner singularity. The capabilities of these methods in capturing the behavior of 2D flows inside the L-shape geometry, for different values of the viscosity (thus for different Reynolds number), is tested.

Presentation of the problem

The flow is directed by an oscillating pressure p in (t) = sin(3t) at the inlet y = 1 and a constant pressure at the outlet x = 1 denoted by p out = cte [START_REF]Fenics project[END_REF]. The flow boundary condition is set as free for the velocity at the inlet and outlet, and no-slip for the rest of the boundary. We note that the flow topologies also depend on the specific cavity shape, wall motions, and the Reynolds number Re. The geometry of the L-shaped cavity considered in the present example is shown in Fig. (3.15) and the length of the cavity L is fixed to 1.

Spatial discretization

For spatial discretization, the finite element method is employed using FEniCS. We use piecewise continuous quadratic elements for velocity and piecewise continuous linear elements for pressure. A mesh designed for the L-shape geometry, that can be found on the FEniCS website [START_REF]Fenics project[END_REF], is used here.

Simulation, implementation and results

For more simplicity the density of the fluid is fixed to ρ = 1, the kinematic viscosity is fixed to ν = 0.01, the time step for the classic time integration RK1 scheme is ∆t = 10 -3 s. Before proceeding with the continuation of the temporal perturbation schemes, we are interested to see the behavior of the norms of the velocity and pressure series terms in Fig. (3.16) at time t = 0.0002s. It can be observed that the norms have extremely large values, reaching 10 100 for N = 30. The simulation is computed up to a final time t = 5s, the truncation order is fixed to N = 5. The number of degrees of freedom for the used mesh is ndofs = 1147. We apply ANM, BPL, and IFS for this example with continuation. The momentum and the continuity residuals tolerances denoted respectively by M and c , are fixed to M = 10 -4 and c = 9.898 • 10 -2 . The appropriate fluid behavior and velocity vectors are well respected using these methods as shown in Fig.

Comparisons

We are interested in comparing the efficiency of the methods for this test case first between the classical methods and the time perturbation-resummation methods, and then between these latter for an arbitrary final time of t = 2s and a viscosity ν = 0.01. For the sake of comparison, the CPU number is set to 1 in FEniCS and the generalized minimal residual method (GMRES) [START_REF] Saad | GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF] which is an iterative method for the numerical solution of an indefinite nonsymmetric system of linear equations, is used for the resolution. The residual criterion is used for continuation, and dichotomy is used to find t max . The classical schemes residual norms are calculated and will be used as the prescribed residual norm tolerances.

Classical scheme:

The classical scheme RK1 with a fixed discrete-time step ∆t = 10 -3 s needs 2000 steps to reach the final time with a total CPU time equal to 668.82s. The number of degrees of freedom for the used mesh is ndofs = 1147. The residual momentum norm for RK1 during the calculation is ||Res M || 2 ≈ 10 -4 and the continuity residual norm is ||Res C || 2 ≈ 9.898•10 -2 . These values of the residual norm are used as the momentum and continuity residuals tolerances denoted respectively by M and c , to perform comparisons of perturbation methods. Therefore, for the following comparisons, M and c are fixed to M = 10 -4 and c = 9.898 • 10 -2 .

Time perturbation-resummation methods:

Statistics to reach t = 2s for the Navier-Stokes equations on an L-shaped domain, using ANM, BPL and IFS with different truncation orders N , are presented in Tab. (3.11), Tab.(3.12) and Tab. (3.13). The viscosity is ν = 0.01. The number of degrees of freedom is ndofs = 1147. The number of Gauss used for BPL is N G = 6. For this example, one can see from Tab.(3.11), Tab.(3.12) and Tab.(3.13) that the t max is almost constant when the truncation order varies for each of these methods. BPL has the largest t max , the fewest number of steps, and the smallest CPU time. One can as well notice that the computation of series terms for the Navier-Stokes equation is considered time-consuming. Concerning IFS, we only made the statistic for N = 5 because the code diverges for large truncation orders. 3.13 -Statistics to reach t = 2s for the Navier-Stokes equations on an L-shaped domain, using IFS with N = 5, for a fixed momentum residual tolerance M = 10 -4 and a continuity residual tolerance c = 9.898 • 10 -2 . The viscosity is fixed to ν = 0.02.

N Steps

Comparisons between RK1, ANM, BPL, and IFS:

If we compare the results of ANM, BPL, and IFS in Tab. (3.11), Tab.(3.12) with the classical scheme RK1 that requires 2000 steps and a total CPU time of 668.82s to reach the final time, one can deduce:

• For the different truncation orders tested, BPL requires fewer steps than the classical scheme RK1, whereas ANM and IFS require significantly more steps than RK1 (2000 steps).

• The CPU time of BPL for different truncation orders is less than that of the classical scheme (which is 668.82s). For example, for N = 5, BPL needs 162.57s to reach the final time, whereas the classical scheme RK1 takes 668.82s for the same set of parameters.

• For N < 10, ANM requires less time than the classical scheme, while for N > 10, the classical scheme becomes faster.

• IFS is slower than the classical scheme RK1 and requires many more steps for even a small truncation order.

Study for different mesh size

We are now interested in investigating the sensitivity of these methods to the mesh size. We perform the study for a fixed truncation order of N = 5 and three different refining mesh sizes, evolving three different degrees of freedom ndofs. The residuals of the classical scheme are calculated and used here also as residual tolerance for the temporal perturbation resummation methods computation.

Tab. (3.14) shows the typical effect of mesh refinement on the total CPU time of the methods. This table shows that IFS is the most sensitive to the mesh change. BPL appears to be slightly more sensitive to the mesh change than ANM and RK1, but it is still the fastest in terms of CPU time among the other methods. Table 3.14 -The total CPU time required to reach t = 2s for the Navier-Stokes equations on an L-shaped domain, using RK1 with discrete-time step ∆t = 10 -3 s, ANM, BPL and IFS for N = 5. Different number of degrees of freedom ndofs are considered. The viscosity used is ν = 0.01.

Test problem 3: Lid-driven cavity by MAC method

The Lid-driven cavity example is a primordial problem in fluid dynamics that has attracted the attention of the scientific community for a long time. It was first suggested by Ghia et al. [START_REF] Ghia | High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method[END_REF]. Over the years, the problem has given rise to a large number of papers, the majority of the latter have been devoted to the development of computational algorithms for solving the problem, focusing on the two-dimensional rectangular or threedimensional cubic form of the problem. For the sake of brevity, we will limit our analysis in this example to calculations in two-dimensional rectangular. This problem serves as a classical reference model for evaluating the effectiveness of numerical methods as well as for studying fundamental aspects of incompressible flows and the pure physics of two-dimensional flows. It is a suitable problem to be tested for several reasons. First, there is a wealth of literature to compare with. Second, the laminar solution is steady. Third, the boundary conditions are very easy. In addition, the fully developed flow exhibits almost all fluid mechanics phenomena, with increasingly complex aspects emerging as the Reynolds number increases. One can cite AbdelMigid et al. [START_REF] Abdelmigid | Revisiting the lid-driven cavity flow problem: Review and new steady state benchmarking results using GPU accelerated code[END_REF], who provides a recent and comprehensive review of the literature on this subject, presenting and discussing the work of several authors. Note that the good data set for comparison is that of Ghia et al. [START_REF] Ghia | High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method[END_REF], as it includes tabular results for different Reynolds numbers. These simulation results are obtained by the non-primitive variable approach.

Goal of this example

The Lid-driven cavity is chosen to test the effectiveness of time perturbation and resummation methods, as it is a common benchmark problem in fluid dynamics and has been long used as a validation test for new numerical methods. The goal of this example is to test the capabilities of time perturbation and resummation methods in capturing the behavior of 2D flows inside confined volumes for different Reynolds numbers. We note that this example does not aim to present an optimized numerical model neither for CPU time nor for computational memory consumption. An indication of the comparison between the perturbation methods will be presented.

Presentation of the problem

The Lid-driven cavity problem is modeled by the Navier-Stokes equation presented in Eq.(3.37) and Eq. (3.38). It has no-slip boundary conditions on all walls. On the bottom wall and the side walls, u = (0, 0), while on the top wall (the lid), u = (U = 1, 0), where the driving velocity U in x-direction is one of the input parameters of the simulation that is chosen to be 1m/s. It consists of a cubic cavity, whose lid (the top layer) has a uniform velocity, and whose interior is filled with a fluid undergoing incompressible flow. Here, the cavity of size [0, 1] is considered. Note that, we simplify the problem by not considering the gravitational field force. The schematic of Lid-driven cavity flow is presented in Fig. (3.20). Figure 3.20 -Schematic of the Lid-driven cavity flow [START_REF] Huang | Simulation of Lid-Driven Cavity Flow with Internal Circular Obstacles[END_REF].

Spatial discretization

Regarding spatial discretization, a wide variety of different numerical methods have already been applied to the Lid-driven cavity problem. Here, we use the most popular finite difference method centered on a shifted grid (also known as the MAC mesh or MAC method) for the approximation of the Navier-Stokes equations. A brief background on the MAC method is presented in the following.

Background-MAC method:

The MAC method is well-known for its simplicity, its efficiency, and its remarkable mathematical properties. It is a finite difference technique based on a staggered grid. It was developed by Harlow and Welch [START_REF] Harlow | Numerical Calculation of Time-Dependent Viscous Incompressible Flow of Fluid with Free Surface[END_REF] in 1965 specifically for free surface flows. They used a particle marker to mark the cell containing the fluids and used interpolation to track the movement of the surface. One of the most notable features is the use of virtual Lagrangian particles, whose coordinates are saved and which move from one cell to the next based on the latest evaluated speed field. Some of the notable characteristics of the MAC method include the superposition of mesh rectangular cells on the flow domain and the use of virtual particles to represent the existing fluid. The particles move from one cell to another as a function of flow velocities, keeping track of the fluid location. If a cell contains a certain number of particles, it is considered to contain fluid. Thus, in the MAC method, the flow domain is divided into a certain number of cells that are labeled as empty, full, surface, or border cells, and the labels must be updated as time evolves. The idea of MAC, Marker, and Cell, is to place the unknowns (u, p) = (u, v, p) in different locations. Fig. (3.21) shows an illustration of the MAC mesh near the border [START_REF] Samelson | Surface Pressure Poisson Equation Formulation of the Primitive Equations: Numerical Schemes[END_REF]. Specifically the pressure p is located in the center of each cell and evaluated at square points (i ∓ 1/2, j ∓ 1/2), the x-component velocity u on the middle points of vertical edges at the triangle points (i, j ∓ 1/2) and the y-component velocity v on the middle points of horizontal edges at the circle points (i ∓ 1/2, j).

Simulation, implementation and results

Here, the time integration scheme Euler (RK1) with a discrete time-step ∆t = 10 -4 s is used for this problem. The method has been numerically implemented with a solver coded in C++ language. We note that Meijer-G has not been tested for this problem (Meijer-G function is not available on C). As already mentioned, the driving speed and the fluid density ρ are fixed to 1 in this problem. The following grid size 64 × 64 is used. Concerning the viscosity of the fluid, various parameters leading to different Reynolds numbers have been tested. The obtained results will be analyzed in the following.

The truncation order is fixed to N = 10. The simulation is running until a final time t = 760s by time perturbation method ANM, BPL, and IFS with continuation. We will first present the velocity stream function for different Reynolds numbers (Re) to show the ability of time perturbation and resummation methods to track the correct behavior of the fluid. Then, we will perform a detailed study of these methods, for this example for a Reynolds number fixed to Re = 12500. The velocity stream function is illustrated in Fig. (3.22) for various Reynolds numbers. These velocity curves match the results obtained by Wahba [START_REF] Wahba | Multiplicity of states for two-sided and four-sided lid driven cavity flows[END_REF] and Arumuga Perumal and Anoop Dass [START_REF] Perumal | Multiplicity of steady solutions in twodimensional lid-driven cavity flows by Lattice Boltzmann Method[END_REF] for very similar Reynolds numbers. The flows correspond to those found in the literature [START_REF] Ghia | High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method[END_REF]. Therefore, as can be seen, ANM, BPL, and IFS allow to determine the solution up to very high Reynolds numbers and show the same trend and same flow behavior. A sequence of vortices appears as the Reynolds number increases. As can be observed in Fig. (3.22), the primary vortex forms near the cavity geometric center, while secondary vortices form near the cavity upper left and lower right corners. Higher Reynolds numbers lead to a much greater dispersion of the result, in particular for larger node arrangements. The flow is steady up to a Reynolds number close to 8000 and then becomes periodic in time. Detailed numerical study for Re = 12500: A detailed numerical study for ν = 0.00008, which lead to Re = 12500, is presented.

Classical scheme:

Evolution of the v component velocity and the phase portrait for a probe point P1, located in the lower left corner, are presented in Fig. (3.23) obtained by the classic RK1 scheme. As can be seen in this figure, the limit cycle is reached. Therefore, one might wonder whether time perturbation-resummation methods are able to reproduce the solution with such a large number of Reynolds, for a large final time. We have verified that the fluid behavior is well respected using these methods, by the velocity stream function plots. Two tests are carried out to respond to the question. The first one is for an initial time t 0 = 100s, i.e. a starting point that is not located yet on the limit cycle. The second one is for an initial time t 0 = 400s i.e a starting point that is located on the limit cycle.

Time perturbation and resummation methods for an initial time t 0 = 100s

and Re = 12500: The truncation order is fixed to N = 10. The simulation is running from an initial time t 0 = 100s until a final time t = 700s by time perturbation method ANM, BPL, and IFS with continuation. The velocity solution v obtained by ANM based on the continuation proposed by Cochelin [START_REF] Cochelin | A path-following technique via an asymptotic-numerical method[END_REF], with a tolerance δ < 10 -32 (the choice of δ will be presented in details in the next subsection), is displayed in Fig. 3. Time perturbation and resummation methods for an initial time t 0 = 400s and Re = 12500: For an initial point t 0 = 400s that is already located on the limit cycle, we are interested in determining if the methods maintain the solution on the limit cycle for a large time, for Re = 12500. With continuation: ANM with continuation is now applied from t 0 = 400s using a truncation order of N = 10 and a tolerance δ = 10 -32 . As can be seen in Fig. (3.25), the ANM solution is shifted from the RK1 solution. The limit cycle is no longer respected. Note that the evolution of the validity domain t max of ANM, corresponding to this choice of δ = 10 -32 is presented in Fig. (3.26). It oscillates between 0.00009s and 0.00011s. Therefore, to get a smaller t max than 0.00011s, one need to decrease the tolerance δ more. We are interested to see if time perturbation resummation methods improve the ANM solution.

Indication of the comparison between RK1, ANM, BPL, and IFS

We have noticed that the residual of the temporal perturbation methods has very large values overtime for this test case. Therefore, we forced the upper right bound of the validity domain t max . The goal is to provide an indication of comparison between the methods and to advance the solution with the continuation technique. Due to the inappropriate residual value, t max was not calculated from the residual criterion, as we have done for all the test cases in this manuscript. The forced t max was chosen to be the most appropriate so that the solution does not deviate from the limit cycle. Therefore, it should be noted that the following comparisons are only for a forced t max .

The selected time interval is [400 ; 405]s, this is the interval to obtain 2 periods. From Tab. (3.15), one can indicate that for this test case, ANM and IFS are much more time consuming than BPL. RK1 seems faster than BPL but we note that this study is just an indication of comparison. We recall that the three methods are efficient in terms of feasibility to reach the final time t = 750s, with a forced corresponding t max for different Reynolds number.

Conclusion of the Lid-driven cavity example

This study was made to test the efficiency and accuracy of the time perturbation and resummation methods. The latter is applied to the Lid-driven cavity problem with a Reynolds number that varies in the low and high range, from 100 to 12500. The calculated results are in good agreement with the literature. A numerical study of two-dimensional (2D) incompressible steady flow in the Lid-driven cavity for Re = 12500 was carried out in detail in the present study. Time perturbation and resummation methods have shown great efficiency in terms of stability, robustness, and accuracy for steady solutions at high Reynolds numbers. The various computations at high Reynolds numbers show the robustness and efficiency of the approximation and the computational method. For the comparative studies, we only indicated the CPU time between the methods, for a forced t max rather than a more real robust study: We forced the upper bound of the validity domain t max to advance the continuation of the solution instead of using the residual criterion. This is due to the large residual value of the time perturbation methods over time. Since the starting residual is not optimal, we need to develop a method to have a corrector, and then have an appropriate value for the residual for every step.

Conclusion

In this chapter, a study was conducted on the performance of temporal perturbation and summation methods (ANM, BPL, IFS, and MG) to numerically solve the heat equation, Burges equation, and Navier-Stokes equations. Various tests have been carried out in two-dimensional geometry. We used the classical central finite difference method, the MAC method, and the finite element method for spatial discretization. For the Navier-Stokes equations, we started with simple configurations where the analytical solution is known so that we could compare the solutions and evaluate errors (e.g. The Poiseuille problem). Then, we performed the simulations on more common configurations such as the flow around a domain of shape L and the flow of the cavity driven for different Reynolds numbers. In all of the numerical tests conducted in this chapter, it can be proven that ANM, BPL, and IFS have demonstrated their ability to provide appropriate solutions for PDEs even for a large final time, for high numbers of degrees of freedom and when changing the parameters of the problem (the grid size, the viscosity, the Reynolds number). The obtained results were in good agreement with the literature. We have shown also that these methods can be associated with various spatial discretization techniques. In our time computation studies, we found that when we compare the time perturbation and resummation schemes for the same comparison criterion and parameters, we revealed that: IFS is the most expensive in terms of time, BPL was faster than ANM for all of the examples presented in this chapter with a larger t max requiring fewer steps to reach the final time. When compared to the classic RK1 scheme, RK1 was significantly faster than IFS in most cases. For the processed heat equation and the Navier-Stokes case on an L-shaped domain, RK1 was much slower than ANM and BPL. In the Burgers equation case, RK1 was slightly faster than BPL but much faster than ANM. Because the residual had incorrect values in the LID-driven cavity case for Re = 12500, only a comparison indication was presented rather than a true comparison based on the residual criterion. According to the studies, calculating the series terms is time-consuming for some cases of the Navier-Stokes problems. As a result, we must devise a method to accelerate the calculation of series terms for time perturbation and summation methods of Navier-Stokes equations. In terms of changing the grid and the mesh size, we observed that IFS is the most sensitive to the mesh. Mesh sensitivity is about the same for ANM and BPL. For some test cases with a very fine mesh, RK1 becomes faster than the other schemes.

Regarding the Meijer-G approximant, we wanted to explore the latter for temporal PDE problems. The problem of the 1D heat equation has been carried out by Meijer-G with the finite difference method. For this example, we have shown the efficiency of the Meijer-G for a series vector using the Gram-Schmidt procedure. Other examples, on the other hand, revealed some issues for specific domain areas and truncation orders. The Meijer-G sum was discovered to behave erratically in some areas. The cause is still being investigated. This could be due to the small value of t max in the Meijer-G sum, or to the fact that the Meijer-G function implementation does not always provide correct answers in some specific parameter cases. This is also mentioned in Viktor T. Toth

TOWARDS APPLICATION TO FLUID-STRUCTURE INTERACTION PROBLEMS

We recall that our long term future goal is to establish the necessary adaptations of this new class of time perturbation and resummation schemes to solve full Fluid-Structure Interaction (FSI) problems with free surfaces. The Navier-Stokes equations and capture interface methods constitute the basis of such problems. Since it is known that the time series for Navier-Stokes equations is divergent, resummation techniques of divergent series are required to obtain accurate numerical solutions for such models. Hence, FSI models based on the coupling of Navier-Stokes equations, followed by the most popular technique for capturing the interface which is the level set, are chosen. To serve our future purpose, progression difficulties sould be done in stages. In Chapter 3, the applications of the time perturbation appraches have been applied on some examples of the Navier-Stokes equations. In this chapter, we limit their study on the level set equation alone because there are many points to consider, test, and validate just for the separate analysis of the level set alone by time perturbation methods before its coupling with the Navier-Stokes equations. Therefore, this chapter aims to focus on time perturbation and resummation methods applied to the level set evolution equation subjected to a prescribed velocity field. We will first present a brief review of interface capture methods to clarify the motivation for choosing level set methods to be handled by the temporal perturbation techniques. We recall that the Meijer-G resummation will not be applied to the level set evolution equations since the numerical resolution of the level set has been written in the C programming language where the numerical access of the Meijer-G function G p,q m,n is not available.

Review and motivation

Two main approaches are broadly used in the literature to model the free surfaces namely, Lagrangian methods (e.g. front tracking [START_REF] Unverdi | A front-tracking method for viscous, incompressible, multi-fluid flows[END_REF]) and Eulerian methods [START_REF] Chang | A Level Set Formulation of Eulerian Interface Capturing Methods for Incompressible Fluid Flows[END_REF]- [START_REF] Scardovelli | Direct numerical simulation of free-surface and interfacial flow[END_REF].

In Lagrangian methods, interface location is explicitly known and therefore the interfacial boundary conditions can be easily implemented. However, a primary drawback of this method is that once the interface topology starts changing, "surgical" procedures are needed to smooth and distribute the interfacial elements used to represent the interface [START_REF] Raessi | A volume-of-fluid interfacial flow solver with advected normals[END_REF]. Procedures of the Lagrangian methods are problematic when interfaces merge or split [START_REF] Osher | Level set methods and dynamic implicit surfaces[END_REF], specifically in three dimensions (3D), [START_REF] Osher | Level set methods and dynamic implicit surfaces[END_REF]. On the other hand, in Eulerian methods, interface location is implicitly represented by a scalar function in general on an Eulerian grid, which acts as an indicator and, as such, allows interfaces to merge or break with relative ease. A fixed numerical mesh is used to evolve the interface by solving a transport equation involving this indicator function (for example, the level set function or the VOF function). The level set method, introduced in 1988 [129] is a popular Eulerian numerical approach to implicitly evolve the interface using a smooth (Lipschitz continuous) function (denoted by φ) defined over the entire physical domain. The initial interface is represented as the zero level set function (φ 0 ) and it propagates naturally with the evolution of this zero level set. Among the advantages of this method over others, one may cite its ability to handle topological changes such as breakup and merging, as well as its effective parallelization, the simple generalization to three dimensions without any complications in numerical processing, and its relatively easy implementation. It has been successfully used for various fields of applications [START_REF] Osher | Level set methods and dynamic implicit surfaces[END_REF], [START_REF] Sethian | Level set methods and fast marching methods: evolving interfaces in computational geometry, fluid mechanics, computer vision, and materials science[END_REF].

The first goal of this study is to analyze the behavior and efficiency of time perturbation and numerical resummation methods using techniques that can handle all the topological changes of the interface. Therefore, in this chapter, benchmarks level set problems, are studied by time perturbation and resummation methods. As during the resolution, the level set function rapidly loses its property, several methods to reinitialize this function are proposed in the literature [START_REF] Hysing | The Eikonal equation: Numerical efficiency vs. algorithmic complexity on quadrilateral grids[END_REF]. Therefore, the second goal of this chapter is related to a discussion of the reinitialization procedure in the context of solutions that are continuous in time.

Level set method

In this section, we introduce the principle of the level set method, the reinitialization technique, and the classical numerical resolution of those latter.

Principle of method

The level set method represents an interface Γ as the zero contours of a smooth sign distance function, φ(x, t), called the level set function. The latter is defined over a domain Ω (for example a two-dimensional domain 2D) in such a way that:

         φ(x, t) < 0 ∀ x ∈ Ω 1 φ(x, t) > 0 ∀ x ∈ Ω 2 φ(x, t) = 0 ∀ x ∈ Γ (4.1)
as shown in Fig. (4.1). Note that x = (x, y) since Ω is a two-dimensional computational domain. One of the most obvious observations is that the signed distance function is a regular function that allows to know immediately in which medium the point x = (x, y) is located.

The interface evolution in a given spatial time-dependant velocity field u(x, t) where u = (u, v), is controlled by solving the following transport equation:

∂φ ∂t + u • ∇φ = 0. (4.2)

Properties

The properties of this method are as follows:

• The location of the interface is directly known since the sign of the function φ(x, t) makes it possible to know immediately in which medium the point x = (x, y) is located, as shown in Fig. 4.1.

• The level set method permits the direct access to the local geometrical characteristics of the interface, namely the normal n and the curvature κ

n = ∇φ |∇φ| and κ = ∇n (4.3)
Thanks to the level set function, the normal n and the curvature of the interface κ are calculated in the whole domain. The normal vector n is oriented in the positive direction of the contours and its calculation is simplified using the distance property of the level set function.

• The level set function φ(x, t) is a regular function:

∀x ∈ Ω, |∇φ(x, t)| = 1 (4.4)
this is also the property of the signed distance function. When transporting a scalar field, this property is very interesting from a numerical point of view. It is always tricky to transport a discontinuous scalar for reasons of accuracy, stability, and diffusion of the numerical schemes. Therefore, when the level set function has a strong local curvature, the regularity property Eq.(4.4) may be lost.

Reinitialization technique

In this subsection, we present one of the major drawbacks of the level set method which is the loss of its distance property, and we describe a technique to overcome this problem.

The challenge of using the level set method is that the numerical solution of Eq.(4.2) rapidly loses its signed distance property during larger time computations. This means that the value of φ does not remain the shortest distance between the point and the interface: Eq.(4.4) will be no longer valid. To overcome this problem and maintain φ as a signed distance function, the latter is reinitialized after a certain time, after each iteration for example. Remark 4.1.1 The interface may move slightly during the reinitialization of the level set function. This is not desirable, but it is almost unavoidable. Because the interface moves slightly every time we use the reinitialization technique, we conclude that reinitialization should only be used when absolutely necessary. Too much reinitialization of the interface movement will have a significant impact on the final results, and too few reinitialization will ruin the good properties of the level-set function. It should also be noted that it is not necessary for all level sets to be signed distance functions, but the zero level set and a small area next to it are enough, as mentioned in [START_REF] Johansson | Implementation of a standard level set method for incompressible two-phase flow simulations[END_REF].

Several ways of reinitialization are proposed in the literature. One may cite the fast marching method, the fast sweeping method, the algebraic Newton method, the hyperbolic PDE, and also a "brute force" approach. The interested reader is referred to [START_REF] Hysing | The Eikonal equation: Numerical efficiency vs. algorithmic complexity on quadrilateral grids[END_REF] where the strength and weaknesses of each of these reinitialization methods are compared.

The most popular reinitialization technique is the PDE method proposed in [START_REF] Sussman | A Level Set approach for Computing Solutions to Incompressible Two-Phase Flow[END_REF] and is used in this study. It works as follow: a correction is introduced solving the so-called reinitialization equation [START_REF] Sussman | A Level Set approach for Computing Solutions to Incompressible Two-Phase Flow[END_REF] for few steps in a fictitious time τ (not a real physical time):

∂φ ∂τ + S(φ 0 ) (|∇φ| -1) = 0 (4.5)
where

S(φ 0 ) = φ 0 φ 2 0 + h (4.6)
is a smoothed sign function, φ 0 = φ(τ = 0) is the initial zero level set function and h is proportional to the grid size or smoothing distance. Classically, for each step size ∆t, several iterations of the reinitialization algorithm are needed to ensure the correct distribution of the level sets.

Classical numerical resolution

Various numerical experiences have shown that the discretization schemes of the transport equation Eq.(4.2) and the reinitialization equation Eq.(4.5) must have a sufficiently high order of convergence as well as a high precision to obtain suitable results [START_REF] Tanguy | Développement d'une méthode de suivi d'interface[END_REF].

Spatial discretization:

Since there is a significant sensitivity of level set methods to the approximation of convective terms (∇φ), a great deal of attention has been paid to spatial discretization. The literature review has shown that the majority of authors use the fifth-order of HJ-WENO scheme ("Hamilton-Jacobi Weighted Essentially Non-Oscillatory") [START_REF] Liu | Weighted Essentially Non-oscillatory Schemes[END_REF], [START_REF] Jiang | Efficient Implementation of Weighted ENO Schemes[END_REF] to achieve better results.

Temporal discretization: Concerning temporal discretization, it presents difficulties but less than those of spatial discretization. S. Tanguy [START_REF] Tanguy | Développement d'une méthode de suivi d'interface[END_REF] has shown that even schemes of order 2 such as the Bell, Collela, and Glaz [START_REF] Bell | A second-order projection method for the incompressible navier-stokes equations[END_REF] or Adams-Bashford schemes are inadequate. As mentioned in [START_REF] Shu | Efficient implementation of essentially non-oscillatory shock-capturing schemes[END_REF], [START_REF] Yousefzadeh | A Level-Set Immersed Boundary Method for Reactive Transport in Complex Topologies with Moving Interfaces[END_REF], the temporal discretization is performed either with the third or fourth order of the Total Variation Diminishing Runge Kutta (TVDRK) scheme. These schemes, proposed in [START_REF] Gottlieb | Total Variation Diminishing Runge-Kutta Schemes[END_REF], guarantee that no additional oscillations are created when the order of the temporal discretization is increased. In this chapter, we will use the fourth-order accuracy of TVDRK denoted by TVDRK4. Solutions obtained by TVDRK4 and the fifth-order of HJ-WENO, are considered as reference solutions for comparisons in the present work.

Solution with time perturbation and resummation methods

After presenting the classical resolution of the initial problem in the previous section using TVDRK4 for the temporal discretization and the finite difference scheme HJ-WENO for the spatial discretization, this section is dedicated to the development of time series solution of the level set evolution equation. For more clarification and details, we will present the recurrence formulas to obtain series terms for the the two following cases separately:

• Time-independent velocity fields u(x).

• Time-dependent velocity fields u(x, t).

In both cases, the velocity depends on space. The obtained reccurence formulas will be then discretized using the finite difference scheme HJ-WENO which also permit to validate the association of time perturbationresummation scheme to this finite difference scheme. A discussion about the use of the reinitialization technique in the context of continuous solution time is presented at the end of this section.

Recurrence formulas for time-independent velocity field

First, we consider the case where the velocity field u = (u, v) depends on space but not on time. Therefore, only the level set function is represented in terms of power series using time as perturbation parameter:

φ(x, t) = φ 0 (x) + N i=1 φ i (x) ti (4.7)
where N denotes the truncation order, and φ 0 is a known regular solution at an initial time t 0 .

Injecting Eq.(4.7) into the transport equation Eq.(4.2) and equating the respective power of t, yield to:

                     t0 : φ 1 + u • ∇φ 0 = 0 t1 : 2φ 2 + u • ∇φ 1 = 0 t2 : 3φ 3 + u • ∇φ 2 = 0 . . . tN-1 : N φ N + u • ∇φ N -1 = 0
We recall that u and φ i , ∀i ≥ 1 are space dependent. However, for the sake of simplicity of writing, we did not use φ i (x, t) and u(x) in this system of equations.

Reinitialization technique in the context of transients continuous solutions

Since our work is based on time perturbation methods, we will test three ways of using reinitialization within solution that are continuous in time:

1. The classical PDE reinitialization Eq.(4.5) only after the last step of the continuation procedure of ANM i.e. at the final physical time t, 2. The classical PDE reinitialization Eq.(4.5) at the end of each continuation step of ANM, i.e end of each t max , 3. Time perturbation of reinitialization procedure without using a modified version of Eq.(4.5).

The time perturbation reinitialization procedure is introduced here. We recall that the concept of a reinitialization technique is to maintain the level set function φ as a signed distance function, (i.e. |∇φ| = 1 must be respected during the resolution). We note that we have not solved the classical PDE reinitialization equation (Eq.(4.5)) by time perturbation approach because it contains a pseudo time τ . For this reason, our idea to test the reinitialization procedure by temporal perturbation is based on solving the following system (transport equation + property of the level set signed distance function Writing φ as time power series and replace it in the system (Eq.(4.12) and Eq.(4.13)), one can obtain the following recurrence formulas to obtain φ i series terms ∀i = 1, • • • , N : For the case of time-independent velocity field:

φ i = - 1 i u • ∇φ i-1 + i-1 k=0 ∇φ k • ∇φ i-1-k (4.15)
For the case of time-dependent velocity field:

φ i = - 1 i i-1 k=0 (u k + ∇φ k ) • ∇φ i-1-k (4.16)
Finally, note that if φ 0 is a signed distance function, then each series term φ i for i ≥ 1 may not remain so. This property is tracked by evaluating the value of |∇φ i |, ∀i ≥ 1, which is usually moves away from unity.

Numerical validation for time-independent velocity field

To study the behavior of time perturbation techniques and numerical resummation methods for the solutions of level set evolution equations, various academic test problems proposed in [START_REF] Tanguy | Développement d'une méthode de suivi d'interface[END_REF], [START_REF] Rider | Stretching and Tearing Interface Tracking Methods[END_REF]- [START_REF] Couderc | Développement d'un code de calcul pour la simulation d'écoulements de fluides non miscibles[END_REF] are studied in this section. We deal with two problems in detail where there is a large deformation of the interface. We impose a velocity field that depends on space but not on time. The so-called "Vortex flow" and "deformation field" simulations are the most frequently encountered in the literature, so we have logically carried them out. A test case called "Compression/decompression of a circle" is also studied.

Vortex flow

The "Vortex flow" test case, introduced by Bell, Collela, and Glaz [START_REF] Bell | A second-order projection method for the incompressible navier-stokes equations[END_REF], quantifies the performance of an interface tracking method when the interface stretches into a long ligament due to shear flow. This test case is rigorous because, as we will see, the ligament stretches indefinitely over time. The numerical schemes used must be robust enough to ensure that the code does not diverge. Therefore the latter should be very precise to capture the ligament as much as possible when its width approaches that of the mesh, preserving the initial surface. It is very interesting to test the time perturbation and numerical summation methods for the interface tracking method, as it is almost inevitable when simulating two-phase shear flows.

Consider a square domain Ω = [0, 1] 2 containing a circle of radius r 0 = 0.15 and origin (x 0 , y 0 ) = (0.5, 0.75). The initial interface is given by: φ 0 (x, y) = (x -x 0 ) 2 + (y -y 0 ) 2 -r 0 (4.17)

A velocity field u = (u, v) is defined by the following potential ψ:

ψ = 1 π sin 2 (πx) sin 2 (πy) (4.18)
The analytical expressions for the components of u = (u, v), are given by:

u(x, y) = - ∂ψ ∂y = -2

Numerical results of classical schemes

The two classical schemes (TVDRK4 + HJ-WENO) are applied here (see subsection 4.1.3). In this study, studies are conducted with the grid size 128 × 128 and for a discretetime step, ∆t = 10 -3 s for the TVDRK4 scheme to respect the CFL condition.

Numerical results of time perturbation and resummation methods

In this subsection, the effectiveness of the time perturbation and resummation methods is studied on the level set equation. After determining the final recursive relation that permits to calculate φ i series terms (Eq.(4.8)), the ability of ANM, BPL, and IFS to provide the same approximate solutions as the classical method is examined. These schemes are used to evaluate the solution of level set transport equation (Eq.(4.2)) until a final time t = 3s, for the following grid size: 128 × 128. First, the Euclidean norm of series terms at zeroth step (without continuation procedure) is plotted in Fig. (4.3). It can be observed that the norm increases with the truncation order. Then, the continuation is carried up to t = 3s, using a residual tolerance = 10 -4 . The series is computed up to a truncation order N = 20. First, the Euclidean norm of the level set function φ(x, y, t) 2 obtained by ANM, BPL and IFS are compared to the one obtained by TVDRK4 at time t = 3s in Fig. (4.4). No visible difference can be seen from this plot which is the first indicator of the global correctness of numerical time perturbation and resummation integration schemes. It is observed, in Fig. (4.5), that ANM, BPL, and IFS present the same numerical solution at a probe point of the domain, as the classic scheme TVDRK4. This is further supported by no visible distinction for the interface Γ (zero contour i.e {x = (x, y); φ(x, y, t) = 0}) between BPL (in red) and TVDKR4 (in black) in Fig. (4.6), for t = 0.5s, t = 1s and t = 3s. Same interfaces are obtained with IFS and ANM and are not presented in the figure because of the visualization difficulties among the 4 methods. This was verified by a computation of the relative error between the temporal perturbation methods which is of order 10 -15 . One can deduce that time perturbation methods alone and then followed by resummation methods can respect the evolution of the level set and represent well the important spatial variations in the evolution of the interface. The numerical result has shown the ability of the analyticity of the approximation to reproduce the large deformations of the initial circle.

Comparisons

Now we are interested to compare the range of validity (t max ) of these temporal approaches. The comparison is performed in Fig. (4.7) with the residual criteria for a tolerance = 10 -4 . It can be seen that BPL has the largest t max among the other methods and is significantly greater than the time step, ∆t = 10 -3 s, of TVDRK4. One can deduce that the range of validity of time perturbation and resummation schemes is much larger than the time step of TVDRK4. Therefore, they require relatively lesser steps than the classic scheme to reach the same final time. The classical time integration scheme TVDRK4 with ∆t = 10 -3 s needs 3000 steps to reach a fixed final time of t = 3s. The number of steps performed with time perturbationresummation schemes to reach t = 3s, for a truncation order of N = 15, is presented in Tab.(4.1). It is done for different choices of residual tolerance parameters using the same continuation criteria based on the residual. This table shows that the time perturbation method alone (ANM) and time perturbation-resummation methods (BPL, IFS) need much fewer steps than the classical scheme TVDRK4. Further, it is observed that IFS needs more steps than BPL. CPU time is now presented in Tab.(4.2). Those computation have been performed on an Intel(R) Xeon(R) CPU W3565 3.20GHz, using gcc compiler with optimization flag "-03". The reference CPU time is 6.3s for the classical time integration scheme. Total(s) TVDRK4 Total(s)ANM Total(s) BPL Total(s) IFS 6.3 2.5 5.6 10 Table 4.2 -The total CPU time required to reach t = 3s for the level set evolution equation using TVDRK4 with a discrete-time step ∆t = 10 -3 s and ANM, BPL, and IFS with a truncation order N = 15. The grid used is 128 × 128.

Computations using reinitialization technique

As mentioned in sections 4.1.2 and 4.2.3, the level set function rapidly loses its signed distance property over time. To avoid this, a reinitialization technique is applied. The example of Vortex flow is an ideal test case to study the efficiency of the reinitialization procedure. Note that the transport equation of the level set has been constructed such that each iso-contour (or called also level lines) in addition to the interface itself, is advected by the local velocity field. Such a methodology implies that the iso-contours do not remain locally at the same distance from each other when the tangential component of the gradient of the normal velocity to the iso-contours is non-zero. The principle is to correct the position of the iso-contours from the only valid iso-contour which is the interface. The purpose of this correction based on the solution of the transport equation is to reset the level set values from the calculated position of the interface, to restore the property of the signed distance function. For this reason, we will illustrate the iso-contours to validate the efficiency of the reinitialization algorithm. First, we start solving the classical PDE reinitialization Eq.(4.5) by TVDRK4.

TVDRK4 with the classical PDE reinitialization: Firstly, the level lines of the solution obtained by TVDRK4, before and after using PDE reinitialization technique (see Eq. [START_REF] Couderc | Développement d'un code de calcul pour la simulation d'écoulements de fluides non miscibles[END_REF]). Therefore, the result is significantly improved with the application of the reinitialization algorithm.

Time perturbation of reinitialization technique:

For this test, we apply time perturbation of the reinitialization procedure (see Eq.(4.15)). We noticed in Fig. (4.11) that just one iso-contour (in red) next to the interface (in black) is well-positioned but only if a very small ANM tolerance δ is used. One can always observe the bad behavior of others iso-contours in blue which may be due to the application of the reinitialization technique during the entire continuation process. As a result of this study, one can conclude that adapting ANM with the classic PDE reinitialization at the end of each continuation step of the level set evolution equation is the relatively better choice to obtain a good result. Therefore, for the following, only this choice will be considered.

After the comparison of the different behavior of the level lines between different techniques, now a comparison only for the zeroth contour (interface) between the following methods:

• TVDRK4 with PDE reinitialization, • The classical PDE reinitialization Eq.(4.5) at the end of each continuation step of ANM, i.e end of each t max ,

• The "reference" solution given by Rider and Kothe [START_REF] Rider | Stretching and Tearing Interface Tracking Methods[END_REF] by using the analytical particle method. (It is not an exact mathematical solution, since it is obtained numerically. Nevertheless, it is considered as an approximation of the solution as fine as possible.) is proposed in Fig. It can be seen that the solution given by ANM with PDE reinitialization end of each continuation step, is more accurate than the classic scheme TVDRK4 with reinitialization using a discrete-time step ∆t = 10 -3 s. The length of the coil obtained by ANM with PDE reinitialization end of each continuation step (in blue) is nearly the same as that of the "reference" solution (in red) than that of TVDRK4 (in black) in Fig. (4.12). However, the thickness of the plume is slightly larger than that of the "reference" solution. This is due to the problem of the conservation of the mass of the classical reinitialization (and not a numerical error of ANM).

Finally, as a conclusion of this study, one can deduce that adapting the PDE reinitialization end of every step of continuation of ANM is the better choice in terms of the interface (or zeroth contour) and the other iso-contours (level lines).

Compression-Decompression of a circle

The goal of this example [START_REF] Lakhlili | Modélisation et simulation numériques de l'érosion par méthode DDFV[END_REF] is to validate the choice of adapting the PDE reinitialization end of every step of continuation on another example than the "Vortex flow". This test case does not present any deformation of the initial interface. We will not present a full study for this example. Only the validation of the reinitialization end of each step of ANM continuation will be performed. The initial data given in Eq.(4.17) is used also for this problem. We consider the same domain as the first example Ω = [0, 1] 2 , containing a circle of radius r 0 = 0.4 and origin (x 0 , y 0 ) = (0.5, 0.5) successively compressed and decompressed by the velocity field u = (u, v), then by its opposite defined by:

u(x, y) = x 0 -x (4.21) v(x, y) = y 0 -y (4.22)
The analytical solution after compression for an instant t 1 denoted by φ c [START_REF] Lakhlili | Modélisation et simulation numériques de l'érosion par méthode DDFV[END_REF], is given by:

φ c (x, y, t) = (x -x 0 ) 2 + (y -y 0 ) 2 -(r 0 -t 1 ) (4.23)
The analytical solution after decompression for an instant t 2 denoted by φ d [START_REF] Lakhlili | Modélisation et simulation numériques de l'érosion par méthode DDFV[END_REF], is given by:

φ d (x, y, t) = (x -x 0 ) 2 + (y -y 0 ) 2 -(r 0 -t 1 + t 2 ) (4.

Deformation field

The third test case that we will deal with is the deformation field, that is also processed in [START_REF] Rider | Stretching and Tearing Interface Tracking Methods[END_REF] by using the analytical particle method. In this test case, we will also detect a large deformation of the initial interface and the topology. The same initial assumptions of the domain, and radius of the circle are the same as section 4.3.1. The initial data given in Eq.(4.17) is used also for this problem, but the complex velocity field u = (u, v) is defined by the following potential ψ:

ψ(x, y) = 1 4π sin 4π(x + 1 2 ) cos 4π(y + 1 2 ) (4.25)
The analytical expressions for the components of the velocity field u = (u, v) are given by

u(x, y) = - ∂ψ ∂y = sin 4π(x + 1 2 ) sin 4π(y + 1 2 ) (4.26) v(x, y) = ∂ψ ∂x = cos 4π(x + 1 2 ) cos 4π(y + 1 2 ) (4.27)
The fluid elements are deformed by the flow field. Level-set interface solutions for t = 1s and t = 2s obtained by BPL resummation procedure and TVDRK4 are in good agreement in Fig. (4.17). Note that BPL also reproduces the evolution of the interface topology as in Fig. (4.16). The same figures are obtained with ANM and IFS. This proves that these methods can follow this deformation and reproduce the solutions accurately. Finally as a conclusion of this study, one can deduce from these examples that the level set evolution equation with a time-independent velocity field has been successfully handled. Thus, the recurrence formula Eq.(4.8) is validated.

Numerical validation for time-dependent velocity field: "Time-Reversed Flow fields"

The goal of this test case is to validate the recurrence formula Eq.(4.11) in the case of time-dependent velocity field. This test has been also handled in [START_REF] Rider | Stretching and Tearing Interface Tracking Methods[END_REF]. The same initial assumptions as the previous sections are used for this problem, but the velocity field u(x, t) = u = (u, v) is multiplied by cos(ωt), therefore: u(x, y, t) = -2 sin 2 (πx) sin(πy) cos(πy) cos(ωt) (4.28)

v(x, y, t) = 2 sin(πx) cos(πx) sin 2 (πy) cos(ωt) (4.29) where ω denotes the angular frequency.

It should be noticed that any velocity field can be multiplied by cos πt T that made it periodic in time T , where T is the time period when the flow returns to its initial state [START_REF] Leveque | High-resolution conservative algorithms for advection in incompressible flow[END_REF]. The flow is time-reversed when the period T is equal to 2.

Calculation of velocity field series terms

A time perturbation of the prescribed time-dependent velocity field is required. The velocity coefficients u i = (u i , v i ), ∀i ≥ 1 in the recurrence formula Eq.(4.11) should be calculated in order to determine φ i For the sake of simplicity, let: A = -2 sin 2 (πx) sin(πy) cos(πy) (4.30)

Therefore, the first component of the velocity can be written in the following form:

u(x, y, t) = A cos(ωt) (4.31)
Let the change of variable D = cos(ωt) and E = sin(ωt) as proposed in [START_REF] Guillot | A Taylor series-based continuation method for solutions of dynamical systems[END_REF]. This leads to:

   D = cos(ωt) E = sin(ωt) =⇒    dD = -ωE dt dE = ωD dt (4.32)
After writing D and E as a series:

D = D 0 + N n=1 D n tn and E = E 0 + N n=1 E n tn , then
injecting them into the second part of Eq.(4.32), and finally equating terms with identical powers of t: we get for n ≥ 1:

   D n = -ω n E n-1 E n = ω n D n-1 (4.33)
The initial terms are D 0 = cos(ωt 0 ) and E 0 = sin(ωt 0 ) where t 0 denotes the initial time. Therefore,

u n = AD n ∀n = 1, • • • , N (4.34)
The second component of the velocity v n is calculated in the same way. Finally, after identifying u n = (u n , v n ) , series terms φ n can now be computed by Eq.(4.11).

We want to verify the ability of these methods to reproduce the solution continuously also for the case of the time-dependent velocity field.

Numerical results

The truncation order is set to N = 10, the frequency to ω = 0.125π and the residual tolerance to = 10 -5 . The interface given by the classic scheme TVDRK4 with ∆t = 10 -3 s, and ANM, BPL and IFS is indistinguishably visually in Fig. The error between ANM, BPL, IFS and the classical scheme TVDRK4 is around 10 -4 . Evaluating these errors clearly shows the efficiency of time perturbation-resummation methods also in the case of a time-dependent velocity field.

Finally as a conclusion of this study, one can deduce that the level set evolution equation with a time-dependent velocity field has been successfully handled. Thus, the recurrence formula Eq.(4.11) is well validated.

Conclusion

The numerical time perturbation and resummation methods ANM, BPL, and IFS are tested and compared to the classic integration scheme for solving level set evolution equations in terms of efficiency and robustness. As previously stated, the numerical schemes used for the level set must be robust and precise enough to prevent the resolution code from diverging. We recall that the Meijer-G resummation was not applied to the level set evolution equation because the numerical resolution of the level set was written in the C programming language, which does not support numerical access to the Meijer-G function.

It is found that the proposed methods capture very accurately the evolution of the interface locations for a given time-dependent and independent velocity fields (but depending on space). It has been well shown that BPL and IFS tend to improve the validity domain of the ANM. The validity domain of these three methods is much larger than the time step of the classical TVDRK4 scheme.

Regarding the reinitialization technique in the context of solutions that are continuous in time, we have to be careful how to perform it. Following various tests, we determine that the best choice is to adopt a classical reinitialization at the end of each ANM continuation step. This choice is effective and has the same advantages as the classical techniques. Regards the CPU time analysis, it is heartening to note that we have revealed that the CPU time required by ANM is approximately 3 times less than that of TVDRK4 for the same set of simulation parameters. Therefore, for this approach, the temporal perturbation methods seems to be more interesting than the classical schemes.

The testing of different reinitialization trials on level set evolution equations and the validation of numerical time perturbation and resummation techniques and in this chapter and on Navier-Stokes equations in Chapter 3, is an important step that will contribute to the process of their implementation to solve full Fluid-Structure Interaction (FSI) problems, or other coupled nonlinear problems.

CONCLUSIONS AND PERSPECTIVES

Conclusions

This thesis was related to the numerical resolution of time-dependent problems. Since the improvement of numerical methods for solving nonlinear evolutionary problems is currently a growing research field, the goal of this thesis was dedicated to the development of the numerical time perturbation and resummations methods for solving time-dependent problems. The main purpose was to study these schemes, as alternatives to time integration methods. The main idea of time perturbation approach is to represent the unknowns of the problem in the form of a polynomial by taking the time as the perturbation parameter. Therefore, an analytical approximation of the solution is then continuous and differentiable on a temporal validity domain to be determined. The most interesting property of this time perturbation approach is that the computed solutions and their derivatives are obtained continuously in time. However, one must be very careful when manipulating time series, since they can suffer from divergence. For this reason, resummation techniques of divergent series are used.

In a brief introduction, we provided an overview of the topic covered in this thesis, followed by a detailed description of the research study objectives. The work has been divided into four chapters. The first chapter focused on time perturbation methods ANM (Asymptotic Numerical Method) and summation algorithms of divergent series BPL (Borel-Padé Laplace), IFS (Inverse Factorial Series), and Meijer-G approximant (MG) for solving time-dependent problems. This chapter provided the necessary background of these methods, as well as some theoretical concepts and literature results. It also addressed some developments that we have proposed. The developments include a first use of the recent Meijer-G approximant for time-dependent problems, since the latter has been applied for summing divergent series in quantum mechanical problems. According to the literature [START_REF] Mera | Fast summation of divergent series and resurgent transseries from Meijer-G approximants[END_REF], MG is particularly designed for series with a branch cut. However, the encouraging results presented in [START_REF] Mera | Fast summation of divergent series and resurgent transseries from Meijer-G approximants[END_REF] motivated us to explore and develop the use of this technique to the sum of divergent time series. Thus, with this new summation method, we did not limit the analysis in this manuscript to BPL and IFS only. For this reason, we developed some points of this exploration. We introduced the MG algorithm to solve temporal problems, for even and odd truncation orders. We proposed a continuation technique based on the residual criterion, that allows us to obtain an entire solution rather than just one branch of the solution. As a result, the expressions of the derivatives, followed by the residual of MG, have been detailed. We also introduced a new Meijer-G version for the vector series. We also established the expressions of the first and second-order derivatives for BPL and IFS in an explicit way. Finally, we proposed to compare such exact approximations derivatives use with others suggestions in the numerical examples, in order to approximate better the derivatives and obtain the correct solutions. The second chapter focused on the numerical applications of these methods to nonlinear ODEs in solid mechanics. The third and fourth chapters are devoted to their applications on time-dependent PDEs in fluid mechanics. The conclusions that can be drawn from the numerical tests in Chapters 2, 3, and 4 will be divided into two main categories: in terms of efficiency feasibility, and computation effort.

We note first that increasing or decreasing the order of time perturbation-resummation schemes is very simple, one can just need to change the value of a parameter in the code. Concerning the calculation of the series terms using recurrence relations, previous researchs on perturbation methods has shown that rewriting nonlinear equations in quadratic form is recommendable. In this study, we then used this concept to propose simple recurrences relations.

In terms of efficiency and feasibility: After different numerical test cases, one can in principle state that the application of the continuation technique using the exact approximations of the derivatives of BPL, IFS, and MG sums, guarantees a good quality of the solution for different nonlinearity parameters of the problem. Moreover, the solution obtained by the latter is also efficient for long-time intervals. Another proposal of approximating the derivatives and calculating the residuals, than those formulas, has been investigated. The results, however, indicate that the latter may cause some shifts in the solution obtained for simulations over long time intervals. For this reason, we evaluated the residuals using the derivatives formulas that we have developed in Appendices. We start with the conclusion on ANM, BPL, and IFS before moving on to Meijer-G, which has not been applied to all of the test cases in the manuscript (Meijer-G function is not available in all scientific languages). ANM, BPL, and IFS have been successfully applied on ordinary differential equations, using examples of stiff and highly nonlinear problems where classical schemes suffer from many weaknesses. We have chosen interesting examples from the literature that are known to be excellent examples to evaluate the performance of new time integration techniques. Then, we showed that ANM, BPL, and IFS are well suited to long-term simulation. For some high nonlinear cases, classical schemes frequently introduce a small time lag at each period, whereas the latter the numerical methods based on perturbation reproduce more faithfully the limit cycle. Furthermore, these methods are suitable for stiff problems and can deal with cases where the system exhibits extremely fast dynamics. These approaches have also been shown to be effective when the considered problem exhibits for a wide range of degrees of freedom. We demonstrated that they are a viable alternative for solving ODEs.

The encouraging results on ODEs, have motivated us to apply the latter also on the time-dependent PDEs in particular for fluid mechanics problems. These latter have demonstrated their ability to provide appropriate solutions for other time-dependent PDEs such as the Burgers equation, the Heat equation, the Navier-Stokes equation and the level set evolution equation. They have proven their high levels of efficiency in terms of stability, robustness, and accuracy for solutions at high Reynolds numbers. The obtained results are consistent with the literature. The simulation is validated for large final times and for different parameters changes such as the Reynolds numbers, the viscosity, and the grid size. We also showed that these methods can be associated to different spatial discretization techniques (finite element method, finite difference method: HJ-WENO scheme, MAC method). A numerical study of two-dimensional (2D) incompressible flow in the Lid-driven cavity for Re = 12500 was carried out in detail in the present study. We have verified that the fluid behavior is well respected using these methods, by the velocity stream function plots. We recall that our long-term goal is to use time perturbation-resummations methods to simulate FSI problems where equations governing the fluid motion are coupled to equations governing the evolution of the free surface. The Navier-Stokes equations and the level set methods are the central elements of such problems. It is known that the resolution of the level set equation necessitates the use of high-order numerical schemes in order to capture the interface movements and, as a result, respect the mass conservation because this property is unfortunately not inherent in the method. However, it has been demonstrated through the execution of some case studies that combining the temporal perturbation-resummation schemes ANM, BPL, and IFS with the spatial integration schemes of nonlinear terms (HJ-WENO), provides satisfactory results. The perturbation and temporal resummation methods used for the level set have proven to be robust enough to ensure that the solution does not diverge. These latter capture very accurately the evolution of the interface for both time-dependent and time-independent velocity fields. We also emphasized the need to redistribute the level set in order to ensure an accurate calculation of the interface curvature using the reinitialization procedure. Regarding the popular reinitialization technique in the context of time-continuous solutions, we have seen that one has to be careful how to perform it. Following various numerical tests, we have determined that the best choice is to adopt a classical reinitialization at the end of each ANM continuation step. This choice has the same advantages as the classical techniques. If a difficulty is encountered, it is better to decrease the tolerance accuracy of the chosen method.

Regarding the Meijer-G approximant, we wanted to explore the latter for some timedependent ODEs and PDEs. For some cases, we found that Meijer-G provides a better approximation of the solution than BPL, ANM, and IFS by using a lower truncation order. It overcomes the fact that poles appears in Padé's coefficient in BPL and leads to numerical issues for the continuation process in some cases. The Meijer-G continuation is validated, allowing to obtain a complete continuous solution in time on all of the requested intervals, rather than just only one power series representation of the solution that is valid in a specific temporal range of the solution. The proposed continuation is based on the residual criterion. As a result, we presented and validated the Meijer-G sum derivative exact approximations. The solution obtained by MG with continuation, agreed with the exact solution and provided the correct limit cycle. For the heat equation, we were able to associate this method with a finite difference spatial discretization method. In addition, we have shown the efficiency of the Meijer-G approximant for a series vector using the Gram-Schmidt procedure. However, other examples revealed some issues for specific areas of the domain and truncation orders. It was found that the Meijer-G sum did not behave smoothly in some areas. The cause is still being investigated. This could be due to the small value of t max in the Meijer-G sum, or to incompatible parameters in the Meijer-G sum, or to the fact that the implementation of the Meijer-G function does not always give the right answers in some specific parameter cases. This is also mentioned in the technical notes of Viktor T. Toth (for Maple implementation) [104]. In all the document, we have plotted only the real part of the Meijer-G solution that interests us. We have noted that there are some cases where the physical interpretation of the obtained imaginary part is not so obvious and it is not our main goal. Our goal is to obtain an accurate approximation of the solution.

In terms of computation effort: We performed CPU time studies for one example of ODEs (Van der Pol oscillator) and all the studied examples of PDEs using ANM, BPL, and IFS. However, we did not perform CPU time analysis for MG because the study of its effectiveness and feasibility for time-dependent problems, is still in progress. It should be noted that for both small and large problems, IFS is less efficient in terms of computation time than ANM, BPL, and the classical scheme for certain large test cases. Furthermore, it is extremely sensitive to the mesh size since the CPU time increases significantly when we refine the mesh. Finally, one can say that it becomes slower with higher truncation orders. For small size problems, we have noticed that when only a moderate accuracy is required, the classical time integration scheme can outperform ANM and BPL in some cases. However, when high accuracy is required, BPL becomes the most attractive. This is especially true when the system is large. When we compared ANM and BPL, we found that BPL is similar to ANM or slightly faster for small size problems. However, for large problems with high number of degrees of freedom, the BPL algorithm improves ANM and is more than twice as fast. When simulating a high-dimensional problem, such as the Navier-Stokes equations, where the computation of serial terms takes a considerable amount of time, BPL can significantly reduce the computational cost. We found that when comparing the time perturbation and resummation schemes for the same comparison criterion and parameters, we revealed that: IFS is the most expensive in terms of time, BPL was faster than ANM for all examples with a larger t max requiring fewer steps to reach the final time. For solving the level set evolution equation, ANM and BPL was faster than the traditional TVDRK4 scheme. Compared to the classical RK1 scheme, RK1 was significantly faster than IFS in most cases. For the processed heat equation and the Navier-Stokes equation on an L-shaped domain, RK1 was significantly slower than ANM and BPL. For the Burgers equation, RK1 was slightly faster than BPL but much faster than ANM. Because the residual had incorrect values in the case of the LID-driven cavity for Re = 12500, only a comparison indication was presented rather than a true comparison based on the residual criterion. According to the studies, the calculation of serial terms is time consuming for some cases of Navier-Stokes problems. Therefore, we need to design a method to speed up the computation of the series terms for the time perturbation and summation methods of Navier-Stokes equations. With respect to grid and mesh size modification, we observed that IFS is the most sensitive to the mesh size. The mesh sensitivity is about the same for ANM and BPL. For some test cases with very fine mesh, RK1 becomes faster than the other schemes.

Perspectives

We end this manuscript by indicating some potential future research directions that may emerge from the work presented here, building on the results summarized in the previous chapters. The perspectives that can be drawn are divided into two ways of work: theoretical and numerical parts.

Theoretically

As we have already seen, the theoretical notion of the methods was skipped in this manuscript and only the numerical aspects of the temporal schemes are discussed. However, optimizations to the algorithm may be possible by considering some theoretical improvements.

• We recall that one of the conditions for the series being 1-summable is that is a Gevrey series of index one (condition for BPL algorithm presented in Chapter 1). Note that Maillet [START_REF] Maillet | Sur les séries divergentes et les équations différentielles[END_REF] has proven that formal series that are solutions of an algebraic ODE are Gevrey series. Malgrange [START_REF] Malgrange | Sur le théorème de Maillet[END_REF] has extended this observation to the analytic nonlinear ODEs. The theory for PDEs is not quite complete. Only some equations in fluid mechanics have been investigated. For this reason, an addition theoretical examination of the studied equation can be useful to determine the actual Gevrey index and evaluate it numerically using the series coefficients. For some numerical cases treated in this manuscript, the theoretical proof that the Gevrey index is 1 already existed. However, we did not present the Gevrey index studies for remaining test cases.

• We recall that the third step of BPL algorithm is to evaluate the Laplace transform of the Padé approximant in a direction linking the origin to the infinity in the direction θ (where θ is an angle), denoted by d θ = [0, e iθ ∞[, to return to the temporal space. Since the singular directions of the Borel transform are not known a priori, we have chosen θ = 0, therefore this direction is set as the real positive axis in this work; d 0 = R + . We could make a more detailed study of the equations to which we have applied the resummation. This would allow us to find possible singularities in the Borel space and to have more precision on the best direction of integration during the Laplace transformation.

• It is possible to investigate whether there is a relation between the appearance of the imaginary part of the solution obtained by the Meijer-G approximant (see Fig. • Since MG resummation is recent, the convergence properties are unknown. A convergence study of the algorithm at order N = 3 versus order N = 4 is presented in [START_REF] Antipin | Resummation in QFT with Meijer G-functions[END_REF]. If the theoretical convergence properties are more developed, it may be help to find an explicit expression of t max of MG or other error relation without going through the time-consuming residue criterion.

Numerically

From a numerical point of view, some directions that can arise from our work can be listed in the following:

• The computation time for evaluating the series terms of the Navier-Stokes equations is thought to be long. This is due to the numerical difficulties in dealing with big numbers in the series terms. Some optimization is required to accelerate this evaluation.

• The continuation technique of the time perturbation and numerical resummation schemes should be optimized. In some cases, the time step is quite expensive, in particular by considering the number of residual evaluations, especially for IFS.

A less expensive method of quality evaluation should be developed. This should increase the speed of the scheme. Therefore, an efficient continuation criterion that is not based on the residual is requested. By optimizing the algorithm, the efficiency can be increased considerably. The calculation of the residual consumes a significant amount of CPU time in the solution (IFS algorithm). Another technique must be used to evaluate the error.

• Another potential research direction arising from our work in Chapters 3 and 4 is the validation of these perturbation-resummation methods on the coupling of the level set with the Navier-Stokes equations. Many fundamental aspects of this coupling must be investigated. This new class of time integration schemes must be adapted to fluid-structure interaction problems with free surface and real test cases. Moreover, another well-known method for tracking the interface than level set is the volume of fluid method (VOF). It has the same advection equation as the level set function, it is much more difficult to solve it recisely due to the abrupt change in the interface transition region. However, it has a significant advantage over the level set method. It is the preservation of the mass conservation. Coupled techniques are expected to combine the advantages of the both methods, namely good interface representation and mass conservation. One can also see how to study these methods using time perturbation methods.

• In some cases, the residuals of the temporal perturbation and resummation methods had inappropriate values during the time, as we saw in the previous chapter for the LID-driven cavity example. It will be interesting to develop a residual corrector that is suitable for the temporal perturbation aspect (in case the solution quality deteriorates). As a result, coupling the time perturbation methods with a nonlinear corrector is preferable.

• Because the Meijer-G approximant was only used as the first application for time-dependent problems, no CPU time study for MG was presented in this manuscript. As a result, once the effectiveness study of MG is completed, evaluating the CPU time required by this method will be necessary. We notice that for a residual criterion, applying or calling a Meijer-G function G p,q m,n for each time can be timeconsuming.

• Modern computer algebra systems, such as Maple, include a Meijer-G function implementation. This implementation can be strictly symbolic, or it can include a numerical code. This function implementation is not yet available for all languages, such as C or Fortran, and a connection via Matlab, Python, Maple, or Mathematica is required. This link can be time-consuming sometimes.

• As we have seen, for some test cases and areas of domain, Meijer-G did not behave smoothly. This could be due to the small value of t max in the Meijer-G relation sum, or to the fact that the implementation of the Meijer-G function does not always give the right answers in some cases of specific parameters [104]; or that these approximants are affected by spurious poles; or the obtained hypergeometrics vectors are not positive to be able to apply the Euler Gamma function. Therefore, it would be interesting to study the cause.

• In the following article that is in progress ( https://arxiv.org/pdf/2007.00390.

pdf), the linear stability of the BPL integrator has been introduced. It has been shown that BPL is not A-stable as most of the explicit method, but more efficient in solving the stiff problems than classical schemes. Therefore, what about the A-stability study of Meijer-G and the size of the linear stability domain? Abstract: The development of numerical methods for solving nonlinear evolution problems is currently a growing research field. Therefore, the main goal of this thesis is to address several needs regarding the development of time perturbation methods and numerical summation of divergent series used as time integration schemes. We are interested in the Asymptotic Numerical Method, the Borel Padé Laplace summation, the Inverse factorial series, and the Meijer-G approximant. The most interesting property of these approaches is that the obtained solutions are continuous in time. These methods are applied to the ordinary and partial differential equations in mechanics. These include the heat equation, Burgers, the interface tracking equa-tion and the Navier-Stokes equations. These approaches have shown their efficiency over long time intervals on a variety of problems. Thus, numerical tests show that it is possible to reduce the computation time of classical methods for some problems while eliminating the constraints related to the classical scheme. For other problems such as the Navier-Stokes equations, we have noticed that in some cases an optimization must be proposed to improve these approaches. Moreover, a new insight is provided on the use of the Meijer-G approximants on time-dependent problems. Different developments related to this method have been validated and show interesting perspectives.

Figure 1 . 1 -

 11 Figure 1.1 -Illustration of the upper bound of the domain of validity t max .

Remark 1 . 1 . 1 Figure 1 . 2 -

 11112 Figure 1.2 -Example of ANM series solutions using different truncation orders N compared to the exact solution. (Case of a differential equation representing the motion of a simple pendulum).

3 Figure 1 . 3 -

 313 Figure 1.3 -Illustration of the continuation technique.

Figure 1 . 4 -

 14 Figure 1.4 -ANM series solutions for Euler equation using different truncation order N (without continuation), compared to the exact solution.

. 52 )Figure 1 . 5 -

 5215 Figure 1.5 -Approximated solutions for Euler equation using ANM and BPL with two truncation orders N = 15 and N = 20 (without continuation).

  (i) There is no uniqueness of Padé approximant. (ii) The Pole in Padé approximant on R + might appear. (See an example in Fig.(1.6

Figure 1 . 6 -

 16 Figure 1.6 -Example of the existence of a pole in Padé approximant on R + of the BPL solution.

63 ) 1 . 3 . 2

 63132 Theorem Let u( t) be a 1-summable series, A and C two constants implicated by the exponential decay of the prolongation P (ξ) of its Borel transformation and T = 1 C .

  Fig.(1.7)) can surmount the problem of poles in Padé approximants.

Figure 1 . 7 -

 17 Figure 1.7 -Figure extracted from[START_REF] Deeb | Intégrateurs temporels basés sur la resommation des séries divergentes[END_REF] where there is a pole in Padé approximant.

Figure 2 . 1 -

 21 Figure 2.1 -Approximated solution using time perturbation method ANM, time perturbation-resummation methods BPL, IFS and MG for N = 5 compared to the exact solution.

Figure 2 . 2 -Figure 2 . 3 -

 2223 Figure 2.2 -Approximated solution using time perturbation method ANM, time perturbation-resummation methods BPL and IFS for N = 10 and time perturbationresummation method MG for N = 5, compared to the exact solution.

Figure 2 . 4 -

 24 Figure 2.4 -Imaginary and real parts of the Meijer-G solution for N = 5, for a simple harmonic oscillator example.

Figure 2 . 5 -

 25 Figure 2.5 -(a): Approximated solution using time perturbation method ANM, time perturbation-resummation methods BPL, IFS, and MG for N = 5 compared to the exact solution. (b): Relative error between the solution obtained using time perturbation method ANM, time perturbation-resummation methods BPL, IFS, and MG for N = 5 and the exact solution:|reference solutionnumerical solution| |reference solution| .

Figure 2 . 6 -

 26 Figure 2.6 -Approximated solution using time perturbation method ANM and time perturbation-resummation method MG for different truncation orders compared to the exact solution.

  (a) µ = 0.1 (b) µ = 2 (c) µ = 10

Figure 2 . 7 -

 27 Figure 2.7 -Solutions (left) and phases portrait (right) for the free Van der Pol oscillator using RK4 with a discrete-time step ∆t = 10 -3 s, for different values of µ.

1 -

 1 Number of steps required to reach t = 10s for the free Van der Pol oscillator for µ = 2, using time perturbation-resummation method IFS with N = 15 and different values of the residual tolerance with "Series.D" and "Resummation.D".

Figure 2 . 8 -

 28 Figure 2.8 -Approximated solutions (displacement, velocity, acceleration) using time perturbation-resummation method IFS for = 10 -4 , N = 15, for the free Van der Pol oscillator for µ = 2. In black the reference RK4 solution using ∆t = 10 -3 s.

Figure 2 . 9 -

 29 Figure 2.9 -Approximated solutions (displacement, velocity, acceleration) using time perturbation-resummation method IFS for N = 15 and different values of the tolerance, for the free Van der Pol oscillator for µ = 2. In black the reference RK4 solution using ∆t = 10 -3 s.

4 (

 4 , ε=10-4 Series.D, ε=10-6 Resummation.D, ε=10 -, ε=10-4 Series.D, ε=10 -6 Resummation.D, ε=10 -4 (c) acceleration ẍ(t)

Figure 2 . 10 -

 210 Figure 2.10 -Approximated solutions (displacement, velocity, acceleration) using time perturbation-resummation method IFS for = 10 -4 , N = 15, for the free Van der Pol oscillator for µ = 2. In black the reference RK4 solution using ∆t = 10 -3 s.

Figure 2 . 11 -

 211 Figure 2.11 -Approximated solution (displacement, velocity, acceleration) using time perturbation-resummation method IFS for = 10 -4 , N = 15, for the free Van der Pol oscillator for µ = 1000. In black the reference RK4 solution using ∆t = 10 -3 s.

Figure 2 . 13 -

 213 Figure 2.13 -Approximated solution using time perturbation method ANM (series) for N = 15 and δ = 10 -6 , for the free Van der Pol oscillator for µ = 2. In black the reference RK4 solution using ∆t = 10 -3 s.

Figure 2 .

 2 Figure 2.14 -Approximated solution (left) and phase portrait (right) using time perturbation-resummation method BPL for = 10 -4 , N = 15 for the free Van der Pol oscillator for µ = 2. In black the reference RK4 solution using ∆t = 10 -3 s.

  (a) N = 5, µ = 2 (b) N = 10, µ = 2 (c) N = 35, µ = 2 (d) N = 8, µ = 2

Figure 2 . 15 -

 215 Figure 2.15 -Different examples of apparition of poles resulting from Padé approximant of BPL for the free Van der Pol oscillator for µ = 2.

  (2.15d)).A comparison of the solution obtained by the perturbation method ANM, resummation techniques BPL and IFS for N = 8 without continuation, is displayed in Fig.(2.16). As expected and per construction, no numerical poles are observed for IFS. A continuation procedure using IFS is performed using the residual criterion. The residual (Eq.(1.71)) is based on the first and second derivative of IFS presented in Appendix B. The continuation of IFS based on the residual is carried up to a final time t = 15s, using the residual tolerance = 10 -4 and N = 8. Evolution of the IFS solutions are plotted in Fig.(2.17). It can be observed that IFS presents the same numerical solution as RK4. Oscillations with IFS reach also a limit cycle as it can be appreciated in Fig.(2.17b).

Figure 2 .

 2 Figure 2.16 -Approximated solution using time perturbation method ANM (series), time perturbation-resummation BPL, and IFS for N = 8 without continuation, for the free Van der Pol oscillator for µ = 2. In black the reference RK4 solution using ∆t = 10 -3 s.

Figure 2 . 17 -

 217 Figure 2.17 -Approximated solution (left) and phase portrait (right) using time perturbation-resummation method IFS for = 10 -4 , N = 8 for the free Van der Pol oscillator. In black the reference RK4 solution using ∆t = 10 -3 s.

Figure 2 .

 2 Figure 2.18 -Approximated solution using time perturbation-resummation method BPL in the interval time [9900; 10000]s, for the free Van der Pol oscillator for µ = 2.

Figure 2 .

 2 Figure 2.19t max evolution for time perturbation method ANM, time perturbationresummation methods BPL and IFS, for N = 15, and = 10 -4 for the free Van der Pol oscillator for µ = 2.

3 -

 3 Number of steps required to reach t = 10s for the free Van der Pol oscillator for µ = 2, using RK4 adaptive step size with different values of the desire accuracy parameter ∆ 0 .

Fig.( 2

 2 Fig.(2.20) displays the results of the solution obtained by time perturbation method ANM and RK4 for a discrete-time step ∆t = 10 -4 s, for a final time t = 6000s. One can

Figure 2 . 20 -Figure 2 . 21 -

 220221 Figure 2.20 -Approximated solution using time perturbation method ANM for a truncation order N = 15 and RK4 with a discrete-time step ∆t = 10 -4 s, for the free Van der Pol oscillator with µ = 1000.

Figure 2 .

 2 Figure 2.22 -Zoom plot of the approximated velocity using time perturbation method ANM and RK4 with a discrete-time step ∆t = 10 -5 s, for the free Van der Pol oscillator with µ = 1000.

Figure 2 .

 2 Figure 2.23t max evolution for ANM(R), ANM(C), BPL and IFS for N = 25, for the free Van der Pol oscillator with µ = 1000. (R) stands for ranges of validity evaluated using a residual tolerance criterion = 10 -4 , while (C) stands for the criterion proposed by Cochelin [40] using δ = 10 -4 . (b) is a zoom plot of (a).

5 Figure 2 . 24 -

 5224 Figure 2.24 -Influence of the initial time t 0 for the time perturbed RHS: F (t) = A cos(ωt) without continuation procedure. Comparisons with the exact solution (plain black curve) are performed for A = 17, ω = 4 and a truncation order of N = 20.

Figure 2 .

 2 Figure 2.25 -Influence of the truncation order for the time perturbed RHS F (t) = A cos(ωt) without continuation procedure using the introduction of the new variables. Comparisons with the exact solution (plain black curve) are performed for different truncation orders N , for A = 17, ω = 4 and t 0 = 0.

Figure 2 . 26 -

 226 Figure 2.26 -Approximated solution using time perturbation method ANM (series) for different truncation orders N , in the case of the forced Van der Pol oscillator for µ = 2, ω = 4 and A = 17.

Figure 2 . 27 -

 227 Figure 2.27 -Approximated solution (left) and phase portrait (right) using time perturbation resummation method BPL for = 10 -4 , N = 15 in the case of the forced Van der Pol oscillator for µ = 2, ω = 4 and A = 17. In black the reference RK4 solution using ∆t = 10 -3 s.

  (a) A = 20, ω = 10, µ = 2 (b) A = 35, ω = 10, µ = 2 (c) A = 35, ω = 10, µ = 6 (d) A = 20, ω = 5, µ = 2

Figure 2 . 28 -

 228 Figure 2.28 -Phase portrait using time perturbation resummation method BPL for = 10 -4 , N = 15, in the case of the forced Van der Pol oscillator, for different values of A, ω and µ. In black the reference RK4 solution using ∆t = 10 -3 s.

Fig.( 2 Figure 2 . 29 -

 2229 Fig.(2.29) shows the exact solution compared to the solution obtained by the time perturbation method and the MATLAB solver ode45 with the zoom detail, for a final time equal to 2/R, which is t = 20000s. One can see that the solution obtained by the solver ode45 gives oscillatory motion, while the ANM solution is in very good agreement with the exact one. The same result as ANM is obtained by BPL and IFS.

Figure 2 . 30 -

 230 Figure2.30 -Approximated solutions of the hardening spring problem for a final time equal to 50T , by time perturbation method ANM and the third and fourth-order Runge-Kutta method (RK3/RK4) with a discrete-time step ∆t = T /32. Reference solution is presented in[START_REF] Wood | Stability properties of some algorithms for the solution of nonlinear dynamic vibration equations[END_REF] and is in black.

Figure 2 . 31 -

 231 Figure2.31 -Phase portraits of the hardening spring problem for a final time equal to 50T , by time perturbation method ANM and the third and fourth-order Runge-Kutta method (RK3/RK4) with a discrete-time, step ∆t = T /32. The reference solution is presented in[START_REF] Wood | Stability properties of some algorithms for the solution of nonlinear dynamic vibration equations[END_REF] and is in black. (b) is a zoomed plot of (a).

  Fig.(2.32) illustrates the phase portrait of Lorenz system and shows the efficiency of time perturbation method. Now, we consider x 0 = 5, y 0 = 5, z 0 = 4 and ρ = 27, Fig.(2.33) shows the ability of time perturbation methods to reproduce the two well-positioned attractors.

Figure 2 . 32 -

 232 Figure 2.32 -Phase portraits of Lorenz system with x 0 = 5, y 0 = 5, z 0 = 400 and ρ = 350 obtained by time perturbation method ANM and ode45.

Figure 2 . 33 -

 233 Figure 2.33 -Phase portraits of Lorenz system with x 0 = 5, y 0 = 5, z 0 = 4 and ρ = 27 obtained by time perturbation method ANM and ode45..

46 )Fig.( 2

 462 Fig.(2.34) describes the nonlinear single pendulum problem for two cases of initials conditions to illustrate two different motions of the single pendulum.

Figure 2 . 34 -

 234 Figure 2.34 -Figure extracted from [92] describing the nonlinear single pendulum problem.

Figure 2 . 35 -

 235 Figure 2.35 -Angles θ(t) of the oscillating nonlinear single pendulum problem for the following initials conditions θ(0) = 0.0 and θ(0) = 1.999999238456499, by time perturbation resummation method BPL and the third and fourth-order Runge-Kutta method (RK3/RK4) with different discrete-time step ∆t. The reference solution is in black[START_REF] Kim | An improved explicit time integration method for linear and nonlinear structural dynamics[END_REF],[START_REF] Fung | Solving initial value problems by differential quadrature method-part 2: second-and higher-order equations[END_REF].

Figure 2 . 36 -

 236 Figure2.36 -Angles θ(t) of the rotating nonlinear single pendulum problem for the following initial conditions θ(0) = 0.0 and θ(0) = 2.000000761543501, by time perturbation resummation method BPL and the third and fourth-order Runge-Kutta method (RK3/RK4) with different discrete-time step ∆t. The reference solution is in black[START_REF] Kim | An improved explicit time integration method for linear and nonlinear structural dynamics[END_REF],[START_REF] Fung | Solving initial value problems by differential quadrature method-part 2: second-and higher-order equations[END_REF].

  (3.1), the solution obtained by Meijer-G approximant without any continuation is in a good agreement with the exact one until t = 2.1s visually and through error calculation which is in the range of [10 -3 ; 10 -2 ].

Figure 3 . 1 -

 31 Figure 3.1 -Approximated solution for the heat equation using time perturbationresummation method Meijer-G (MG) for N = 5, without continuation at t = 2.1s. The exact solution is in black.

Figure 3 . 2 -

 32 Figure 3.2 -Approximated solution for the heat equation using time perturbationresummation method Meijer-G (MG) for N = 5, with and without continuation at t = 2.5s. The exact solution is in black.

Figure 3 . 3 -

 33 Figure 3.3 -Evolution of the norm of the series terms with the truncation order N at zeroth step (without continuation procedure), for the heat equation 2D. The grid used is 30 × 30.

Figure 3 . 4 -

 34 Figure 3.4 -Sequence of snapshots of the solution of the Diffusion of a Gaussian function obtained by time perturbation method ANM from t 0 = 0s to t = 0.3s.

Figure 3 . 5 -

 35 Figure 3.5 -Evolution of the solution of the Gaussian function for a probe-line from t 0 = 0s to t = 1s.

Figure 3 . 6 -

 36 Figure 3.6 -Approximated solutions using time perturbation method ANM probed along a line (presented in white color in the figure). The classic RK1 solution is in black.

  in FEniCS. The total CPU time evaluated in all the following tables of comparisons is computed via: CPU Total (s) =CPU Total measured in Fenics (s) -(Number of steps × N × 0.9× CPU time of the factorization and Solve (s)) + 0.9 × CPU time of the factorization and Solve(s) (3.24)

  (3.1), Tab.(3.2) and Tab.(3.3) and is depicted in Fig.(3.7). We recall that the classical method RK1 requires 110.05s of CPU time for the same set of parameters (final time, residual norm). BPL requires only 3.68s for an optimized truncation order N = 20 and an approximate CPU time of 15s for the other truncation orders. Based on Tab.(3.1), Tab.(3.2), Tab.(3.3) and Fig.(3.7

  Evolution of the CPU time with the truncation order. Evolution of t max average with the truncation order.

Figure 3 . 7 -

 37 Figure 3.7 -Comparison statistics between ANM, BPL, IFS and RK1 to reach the final time t = 1s for the heat equation 2D, for a fixed residual tolerance = 2.2 • 10 -3 and a viscosity ν = 1. The grid used is 30 × 30.

Figure 3 . 9 -

 39 Figure 3.9 -Evolution of the norm of the series terms with the truncation order N at zeroth step (without continuation procedure), for the Burgers equation 2D. The grid used is 20 × 20.

Figure 3 . 10 -

 310 Figure 3.10 -Comparison statistics between ANM, BPL, IFS and RK1 to reach t = 0.17s for the Burgers equation 2D, for a fixed residual tolerance = 7.3 • 10 -3 and a viscosity ν = 0.02. The grid used is 20 × 20.

Figure 3 . 11 -

 311 Figure 3.11 -Comparison statistics for several mesh sizes n x × n x , for ANM, BPL and IFS to reach t = 0.17s for the Burgers equation 2D, for a fixed truncation order N = 20. Three different mesh size are considered: 20 × 20, 40 × 40, and 80 × 80 (n x = 20, n x = 40 and n x = 80).

Remark 3 . 4 . 1

 341 [START_REF] Weniger | Summation of divergent power series by means of factorial series[END_REF] with the following Dirichlet boundary conditionu k+1 = g k+1 in ∂Ω (3.69)and the Neumann boundary condition∇p k • n = 0 in ∂Ω (3.70) Costin has shown that the formal solution of the Navier-Stokes equation is in general divergent, Gevrey of index 1 and Borel summable[START_REF] Costin | Divergent expansion, Borel summability and three-dimensional Navier-Stokes equation[END_REF].

72 )

 72 Combining Eq.(3.72) with Eq.(3.71) leads to:

Figure 3 . 12 -

 312 Figure 3.12 -Schematic of a 2D Poiseuille flow between two plates [115].

Figure 3 . 13 -

 313 Figure 3.13 -Approximated solutions (the x-component velocity u and the y-component velocity v) using time perturbation method ANM for N = 6 and δ = 10 -4 probed along a line (presented in white color in the figure). The exact solution is in black.

Figure 3 . 14 -

 314 Figure 3.14 -L2-norm of the relative error between the solutions (the x-component velocity u, the y-component velocity v) obtained using time perturbation method ANM for N = 6 and the exact solution.

Figure 3 . 15 -

 315 Figure 3.15 -Schematic of the L-shaped cavity and the problem considered in the current work.

Figure 3 . 16 -

 316 Figure 3.16 -Evolution of the norm of the series terms (velocity and pressure) with the truncation order N at t = 0.0002s for the Navier-Stokes equations on an L-shaped domain. The number of degrees of freedom is ndofs = 1147.

  (

Figure 3 . 17 -

 317 Figure 3.17 -Behavior of the flow in the L-shaped domain obtained by time perturbation method ANM for N = 5, and ν = 0.01. The same figure is obtained by BPL and IFS. The black arrows design the velocity vector in the domain.

Figure 3 . 18 -

 318 Figure 3.18 -Approximated solutions (velocity u) using time perturbation method ANM for N = 5 probed along a line (presented in white color in the figure) for t = 2s, and ν = 0.01. The same figure is obtained by BPL and IFS. The RK1 solution is in black.

Figure 3 . 19 -

 319 Figure 3.19 -Velocity streamlines for different values of the viscosity ν, obtained by time perturbation method ANM. Same results are obtained with time perturbationresummation methods BPL and IFS for t = 2s.

Figure 3 .

 3 Figure 3.21 -The MAC mesh [124].

  (a) ν = 0.01, Re = 100 (b) ν = 0.0025, Re = 400 (c) ν = 0.001, Re = 1000 (d) ν = 0.0002, Re = 5000 (e) ν = 0.0001, Re = 10000 (f) ν = 0.00008, Re = 12500

Figure 3 . 22 -

 322 Figure 3.22 -Velocity stream function plots for different Reynolds numbers Re, obtained by time perturbation method ANM. Same streamlines are obtained by BPL and IFS and are not presented in the figure.

Figure 3 . 23 -

 323 Figure 3.23 -Evolution of the v velocity and phase portrait for the probe point P1(x(i = 5), y(j = 5)) obtained by the classic scheme RK1 with a discrete-time step ∆t = 10 -4 s for a fixed Reynolds number Re = 12500. (b) is a zoom plot of (a).

Figure 3 . 24 -

 324 Figure 3.24 -Approximated solutions (velocity v and pressure p) using time perturbation method ANM for N = 10, ν = 0.00008 (Re = 12500) compared to the classic RK1 solution for different time intervals.

Without continuation :

 continuation Before starting the simulation of the continuation methods, we show the behavior of the velocity series terms u for an initial time t 0 = 400s, a Reynolds number fixed to Re = 12500 and a truncation order that varies from N = 0 to N = 35 in Fig.(D.1), Fig.(D.2) and Fig.(D.3) presented in Appendix D. One can deduce that the norm of the velocity series terms increases with the truncation order and reaches a magnitude equal to 10 38 for N = 35.

Figure 3 . 25 -Figure 3 .Figure 3 . 27 -

 3253327 Figure 3.25 -Approximated solutions (velocity v and pressure p) using time perturbation method ANM for N = 10, δ = 10 -32 compared to the classic RK1 solution for a time interval [400; 500]s.

  's technical notes (for Maple implementation) [104]. The obtained results open a perspective of the application of Meijer-G for temporal problems.Chapter 4

Figure 4 . 1 -

 41 Figure 4.1 -Illustration of a two-dimensional computational domain Ω, showing the evolution of a level set function, φ(x, t) with a zero contour representing the interface, Γ.

  perturbation method. The following notation |∇φ| = 1 stands for a modulus at each discretization point (∀x ∈ Ω, |∇φ(x, t)| = 1). As we want to use perturbation methods, it is proposed to use the following global information: ∇φ • ∇φ = 1 (4.14)

Fig.( 4

 4 .2) illustrates the evolution of the interface shapes at different times t computed by solving equation Eq.(4.2) with homogeneous Neumann boundary conditions without reinitialization and with TVDRK4. These numerical results clearly show the large deformation of the initial interface.

Figure 4 . 2 -

 42 Figure 4.2 -Evolution of level set interfaces shapes using TVDRK4 without reinitialization, from an initial time t = 0s to the final time t = 3s. The grid size is 128 × 128. The discrete time step is ∆t = 10 -3 s.

Figure 4 . 3 -

 43 Figure 4.3 -Evolution of the norm of the series terms with the truncation order N at zeroth step (without continuation procedure), for the level set evolution equation. The grid used is 128 × 128.

Figure 4 . 5 -

 45 Figure 4.5 -Solution for the probe point P1(x(i = 30), y(j = 40)) = (0.2362205, 0.3149606) (i.e φ(0.2362205, 0.3149606, t)) by time perturbation-resummation methods ANM, BPL and IFS, and the classic scheme TVDRK4 with a discrete-time step ∆t = 10 -3 s. (b) is a zoom plot of (a). The solution obtained by ANM, BPL, and IFS is continuous between the points of continuation indicated by symbols in Figure (b).

Figure 4 . 6 -

 46 Figure 4.6 -Level set interface by time perturbation-resummation method BPL (in red) and the classic scheme TVDRK4 (in black) with a discrete-time step ∆t = 10 -3 s. No visible distinction between the level set contours.

Figure 4 . 7 -

 47 Figure 4.7 -Comparison of the range of validity (t max ) of ANM, BPL (for a fixed Gauss number N G = 10) and IFS, for N = 15 and a residual tolerance = 10 -4 . The grid size is 128 × 128.

  (4.5)) with ∆t = 10 -3 s, are presented in Fig.(4.8). It can be observed that the first plot in Fig.(4.8a) shows disturbances of the level sets while those latter are correctly positioned on both sides of the interface in Fig.(4.8b) (according to the literature

Figure 4 . 9 -

 49 Figure 4.9 -Different level lines {0, 0.1, 0.2, 0.3} for solutions obtained by ANM with the classical PDE reinitialization at the end of the last continuation step at a final time t = 3s. ANM truncation order is N = 20. The grid size is 128 × 128. Indication on the colors of the level lines {0, 0.1, 0.2, 0.3}: {0} in black, {0.1} in red, {0.2} in blue, {0.3} in orange.

  (a) with a tolerance δ = 10 -6 . (b) with a tolerance δ = 10 -15 (c) with a tolerance δ < 10 -15

Figure 4 . 10 -

 410 Figure 4.10 -Different level lines {0, 0.1, 0.2, 0.3} for solutions obtained by ANM with PDE reinitialization at each continuation step, with different ANM tolerances δ, at a final time t = 3s. ANM truncation order is N = 20. The grid size is 128 × 128.

Figure 4 . 11 -

 411 Figure 4.11 -Different level lines {0, 0.1, 0.2} for solutions obtained by time perturbation of the reinitialization technique with a tolerance δ = 10 -30 , at a final time t = 3s. ANM truncation order is N = 20. The grid size is 128 × 128.

  (4.12).

  (a) grid size: 128 × 128 (b) grid size: 256 × 256

Figure 4 .

 4 Figure 4.12 -Level-set interfaces at a final time t = 3s. In blue: ANM + PDE reinitialization end of every step of continuation, N = 20, δ = 10 -6 . In red: the "reference" solution presented by Rider and Kothe[START_REF] Rider | Stretching and Tearing Interface Tracking Methods[END_REF]. In black: TVDRK4 with reinitialization with a discrete-time step ∆t = 10 -3 s.

  Fig.(4.13) illustrates the decompression of the interface shapes during the decompression, at different times t computed by solving equation Eq.(4.2) without reinitialization and with TVDRK4.

Figure 4 . 13 -

 413 Figure 4.13 -Evolution of the level set interfaces shapes using TVDRK4 without reinitialization, from an initial time t = 0s to the final time t = 0.5s. The grid size is 128 × 128. The discrete time step is ∆t = 10 -3 s.

Figure 4 . 14 -

 414 Figure 4.14 -Different level lines {0, 0.1, 0.2, 0.3, 0.4} for solutions obtained, at time t = 0.5s with a truncation order N = 20 and a grid size 128 × 128, by (a): ANM without reinitialization, (b): ANM with PDE reinitialization at each continuation step, with δ = 10 -6 , (c): ANM with PDE reinitialization at each continuation step, with δ = 10 -20 , (or a t max of order 10 -3 ). Same figure is obtained by the analytical solution. Indication on the colors of the level lines {0, 0.1, 0.2, 0.3}: {0} in black, {0.1} in red, {0.2} in blue, {0.3} in orange, {0.4} in purple.

Fig.( 4 . 15 )

 415 plots the level set interface using TVDRK4 for ∆t = 10 -3 s for the grid size 128 × 128. The deformed contours of the level set function appear stretched and ruptured. Furthermore, Fig.(4.16) shows the evolution of interface topology inside the whole domain during the time.

Figure 4 . 15 -

 415 Figure 4.15 -Evolution of the level set interfaces shapes using TVDRK4 without reinitialization using a discrete time step ∆t = 10 -3 s, from initial time t = 0s to final time t = 3s. The grid size is 128 × 128.

  (a) t = 0.5s (b) t = 3s

Figure 4 . 16 -

 416 Figure 4.16 -Interface topology in the whole domain using TVDRK4 without reinitialization using a discrete time step ∆t = 10 -3 s, at two different times t. The grid size is 128 × 128.

  (a) t = 1s (b) t = 2s

Figure 4 . 17 -

 417 Figure 4.17 -Level set interface computed by time perturbation resummation method BPL (in red) using N G = 10, N = 15, = 10 -4 and the classic scheme TVDRK4 (in black) with a discrete-time step ∆t = 10 -3 s. The grid size is 128 × 128.

  (

  4.18) for t = 3s.Local errors between the different couples of these methods is evaluated and some of them are displayed in Fig.(4.19), Fig.(4.20), and Fig.(4.21) respectively for t = 1s, t = 2s and t = 3s. The error between ANM, BPL, and IFS is approximately of 10 -14 order.

Figure 4 .

 4 Figure 4.18 -Level-set interface at final time t = 3s. In blue: ANM, in red: BPL, in green: IFS, in black: TVDRK4 with a discrete-time step ∆t = 10 -3 s. No visible distinction between the level set contours.

  (a) Error between TVDRK4 and ANM (b) Error between ANM and BPL (c) Error between BPL and IFS

Figure 4 . 19 -

 419 Figure 4.19 -Errors between the temporal methods for t = 1s.

  (a) Error between TVDRK4 and ANM (b) Error between ANM and BPL (c) Error between BPL and IFS

Figure 4 . 20 -

 420 Figure 4.20 -Errors between the temporal methods for t = 2s.

  (a) Error between TVDRK4 and ANM (b) Error between ANM and BPL (c) Error between BPL and IFS

Figure 4 . 21 -

 421 Figure 4.21 -Errors between the temporal methods for t = 3s.

( 2 . 4 )

 24 in Chapter 2) and the determination of the upper bound of the validity domain t max . If so, in what circumstances? Otherwise, determining the physical meaning of the imaginary part obtained by MG, related to the problem would be interesting.

Titre:

  Étude numérique de méthodes d'intégration temporelle basées sur des techniques de perturbation. Application à des problèmes de mécanique. Mot clés : Perturbation temporelle, Séries divergentes, Borel Padé Laplace, Séries factorielles inversées, Meijer-G approximant, Méthode asymptotique numérique, Équations de Navier-Stokes, Level-set. Résumé : Le développement des méthodes numériques pour la résolution des problèmes d'évolution non linéaires est actuellement un domaine de recherche en pleine expansion. Par conséquent, l'objectif principal de cette thèse est d'étudier et développer les méthodes de perturbation temporelle et de sommation numérique de séries divergentes utilisées comme schémas d'intégration temporelle. Nous nous intéressons à la méthode asymptotique numérique, à la sommation de Borel Padé Laplace, des séries factorielles inversées et des approximants de Meijer-G. La propriété la plus intéressante de ces approches est que les solutions obtenues sont continues en temps. Ces méthodes sont appliquées aux équations différentielles ordinaires et aux équations aux dérivées partielles de la mécanique. Celles-ci comprennent l'équation de la chaleur, de Burgers, l'équation de suivi d'interface et les équations de Navier-Stokes. Ces approches se sont avérées efficaces sur de longs intervalles de temps et pour de nombreux problèmes. Ainsi, les tests numériques montrent qu'il est possible de réduire le temps de calcul des méthodes classiques pour certains problèmes tout en éliminant les contraintes d'intégration et les faiblesses liées aux schémas classiques. Pour d'autres problèmes tels que les équations de Navier-Stokes, nous avons remarqué que dans certains cas une optimisation doit être proposée pour améliorer ces approches. De plus, un nouvel éclairage est apporté sur l'utilisation des approximants de Meijer-G sur des problèmes dépendant du temps. Les différents développements liés à cette méthode ont été validés et laissent entrevoir des perspectives intéressantes. Title: Numerical study of time integration methods based on perturbation techniques. Application to mechanical problems. Keywords: Time perturbation, Divergent series, Borel-Padé Laplace, Inverse factorial series, Meijer-G approximants, Asymptotic numerical method, Navier-Stokes equations, Level-set.

  1.1)): ||u t,N -u t,N -1 || 2 ||u t,N || 2 = tN ||u N || 2 ||u 1 t + • • • + u N tN || 2

	< δ	(1.16)
	N
	where u t,N =
	n=0

  max ]. The range of validity of the rational representation is usually larger than the polynomial one (i.e t max of Padé > t max of ANM ). Consequently, fewer steps of continuation are necessary.

2 u P adé,N ( t) 2 (1.19) with t in [0, t

Remark 1.3.5 With the new approximations of the first and second derivatives of

  Evaluate the validity domain t max such that t max is the last value of t for which Res(t) 2 < .Step 6: Evaluate a new starting point u 0 at a new time t 0 = t 0 + t max .Step 7: Go to Step 1 and repeat the algorithm until reaching the final time t end .

	Algorithm 3 Continuation technique for IFS
	Inputs:
	1. The initial conditions u 0 and u 1 .
	2. The truncation order N .
	3. The residual tolerance or the accuracy parameter .
	4. The final time t end or the domain in which we want to obtain the curve, for example:
	[0, t end ].
	Step 1: Compute series terms until a given order N .
	Step 2: Calculate the sum of IFS denoted by I(t).
	Step 3: Apply the first and second derivatives approximations of IFS to calculate İ(t)
	and Ï(t). (See Appendix B).
	Step 4: Evaluate the residual Res(t) given in Eq.(1.71).
	Step 5:
	2 denotes the Euclidean
	norm. For the next step, we redefine a new starting point u(t 1 ) where t 1 > t 0 , which
	corresponds to the end of the step where relation Eq.(1.72) holds.
	Therefore, the continuation procedure of IFS is summarized by the Algorithm 3.
	BPL
	and IFS in Appendices A, B, one can evaluate the residual and compare not only the IFS
	and BPL solutions to the exact or classical ones but also the velocity and acceleration
	(especially for nonlinear dynamics problems).
	Remark 1.3.6 In this manuscript, we have chosen a criterion based on the residual for
	the following two reasons:
	1. To ensure a good quality of the solution.
	2. To have the same criterion for all the methods and to compare the solutions ob-
	tained by the following schemes ANM, BPL, IFS with the same tolerance value of
	the residue.
	It should be noted that the continuation of ANM by the residual criterion is based on the
	same concept described for BPL and IFS.

u(t) -u 0

  = (α 11 + α 21 t + α 31 t 2 + α 41 t 3 + α 51 t 4 )tu * 1 + (α 22 + α 32 t + α 42 t 2 + α 52 t 3 )t 2 u * 2 + (α 33 + α 43 t + α 53 t 2 )t 3 u * 3 + (α 44 + α 54 t)t 4 u * 4 + α 55 t 5 u *

		.92)
	u 3 = α 31 u * 1 + α 32 u * 2 + α 33 u * 3 ,	(1.93)
	u 4 = α 41 u * 1 + α 42 u * 2 + α 43 u * 3 + α 44 u * 4 ,	(1.94)
	u 5 = α 51 u * 1 + α 52 u * 2 + α 53 u * 3 + α 54 u * 4 + α 55 u * 5 .	(1.95)
	This allows u(t) to be re-written in the new bases u * 1 , u * 2 , u * 3 , u * 4 , u * 5 by replacing these
	expansions in the polynomial form u(t) in Eq.(1.90):
	5	(1.96)
	Polynomials of degree 1, 2, 3, and 4 thus appear. The main point is to replace these
	polynomials by the sum of Meijer-G approximants in the following way:

t 4 u * 4 + α 55 t 5 u *

  (t) is the Meijer-G sum of α 11 + α 21 t + α 31 t 2 + α 41 t 3 + α 51 t 4 (i.e by applying the even Meijer-G algorithm to α 11 + α 21 t + α 31 t 2 + α 41 t 3 + α 51 t 4 ). (t) is the Meijer-G sum of α 22 + α 32 t + α 42 t 2 + α 52 t 3 (i.e by applying the odd Meijer-G algorithm to α 22 + α 32 t + α 42 t 2 + α 52 t 3 ).

	sense:
	u M G (t) -u(t) ≈ α 55 t 5 u * 5 + O(t 6 )	(1.98)
		5 (1.97)
	where:
	-S 1 M G -S 2 M G -S 3 M G (t) is the Meijer-G sum of α 33 + α 43 t + α 53 t 2 (i.e by applying the even
	Meijer-G algorithm to α 33 + α 43 t + α 53 t 2 ).
	-S 4 M G (t) is the Meijer-G sum of (α 44 + α 54 t) (i.e by applying the odd Meijer-G al-
	gorithm to α 44 + α 54 t).
	The two approximations of the curve are asymptotically equivalent in the following

Table 2 .

 2 2.3), while ensuring a high quality of the solution. 4 -Number of steps required to reach t = 10s for the free Van der Pol oscillator for µ = 2, using time perturbation and resummation methods with N = 15 and different values of the residual tolerance .

	Residual tolerance	ANM BPL IFS
	10 -2	24	21	55
	10 -3	30	26	79
	10 -4	36	33	99

Table 3 .

 3 1 -Statistics to reach the final time t = 1s for the heat equation 2D, using ANM with different truncation orders N , for a fixed residual tolerance = 2.2 • 10 -3 and a viscosity ν = 1. The grid used is 30 × 30.

		< t max >	Total (s) Series (s) Residual Series (s) #RES
	5	417 2.398 • 10 -3	23.69	4.11	15.67	9956
	10 285 3.509 • 10 -3	20.05	5.92	11.36	6840
	15 208 6.950 • 10 -3	16.27	5.68	8.54	4992
	20 164 6.098 • 10 -3	14.28	5.59	7.06	3936
	25 136 7.352 • 10 -3	13.47	5.63	6.19	3264

Table 3 .

 3 2 -Statistics required to reach the final time t = 1s for the heat equation 2D, using BPL, with different truncation orders N for a fixed residual tolerance = 2.2 • 10 -3 and a viscosity ν = 1. The grid used is 30 × 30. The number of Gauss used is N G = 6. For IFS, one can see in Tab.(3.3) that most of the total CPU time is spent for the residual loop. This time grows significantly with the truncation order. It can be observed also that IFS is time-consuming especially for high truncation orders N .

			< t max >	Total (s) Series (s) Residual BPL (s) #RES
	5	266 3.760 • 10 -3	21.49	2.66	16.22	6384
	10	79	1.266 • 10 -2	8.75	1.47	6.54	2291
	15	79	1.266 • 10 -2	9.94	2.00	7.21	2271
	20	25	4.000 • 10 -2	3.68	0.86	4.87	729
	25	70	1.429 • 10 -2	11.87	2.95	8.08	2030
	4. IFS:
	N Steps	< t max >	Total (s) Series (s) Residual IFS(s) #RES
	5	512 1.953 • 10 -3	146.87	5.21	136.64	12288
	10 293 3.413 • 10 -3	307.70	5.30	299.53	7032
	15 196 5.102 • 10 -3	482.92	5.36	475.65	4704
	20 137 7.299 • 10 -3	729.27	7.71	719.54	3288
	25	92	1.086 • 10 -2	704.04	3.96	689.68	2668

Table 3 .

 3

3 -Statistics required to reach the final time t = 1s for the heat equation 2D, using IFS, with different truncation orders N for a fixed residual tolerance = 2.2 • 10 -3

Table 3 .

 3 5 -Comparison of the total CPU time for several mesh sizes n x × n x , for ANM, BPL, IFS, and RK1 to reach the final time t = 1s for the heat equation 2D, for a fixed truncation order N = 20. Three different mesh size are considered: 30 × 30, 60 × 60, and 120 × 120 (n x = 30, n x = 60 and n x = 120). Comparison of < t max > for several mesh sizes n x × n x , for ANM, BPL and IFS to reach the final time t = 1s for the heat equation 2D, for a fixed truncation order N = 20. Three different mesh size are considered: 30×30, 60×60, and 120×120 (n x = 30, n

	x ndofs	Res 2	RK1	ANM	BPL	IFS
	30	961	2.2 • 10 -3 110.65	14.28	3.68	729.27
	60	3721 1.9 • 10 -3 136.81 125.54	72.11	2886.32
	120 14641 3.7 • 10 -2 238.13 1277.80 683.24 13010.75

x = 60 and n x = 120).

Table 3 .

 3

	present work

6 -Comparison of the computed values of u obtained by ANM, BPL, and IFS (denoted by the "present work") for ν = 0.02 at t = 0.625s, with different references in the literature. The grid size is 20 × 20.

Table 3 .

 3 8 -Statistics to reach t = 0.17s for the Burgers equations 2D, using ANM with different truncation orders N , for a fixed residual tolerance = 7.3 • 10 -3 and a viscosity ν = 0.02. The grid used is 20 × 20.

		< t max >	Total (s) Series (s) Residual Series (s) #RES
	5	528 3.220 • 10 -4	64.27	24.51	33.33	10560
	10 336 5.070 • 10 -4	60.05	34.31	21.87	6720
	15 245 6.950 • 10 -4	63.57	44.44	16.41	4900
	20 192 8.855 • 10 -4	70.01	54.35	13.36	3844
	25 171 9.987 • 10 -4	86.36	67.47	12.87	3660

Table 3 .

 3 The t max values are higher than those of ANM (see Tab.(3.8)). 9 -Statistics to reach t = 0.17s for the Burgers equations 2D, using BPL with different truncation orders N , for a fixed residual tolerance = 7.3 • 10 -3 and a viscosity ν = 0.02. The grid used is 20 × 20. The number of Gauss used is N G = 6. For IFS, one can see in Tab.(3.10) that most of the total CPU time is spent for the residual loop. This time grows significantly with the truncation order. It can be observed also that IFS is time-consuming especially for high truncation orders N .

	N Steps	<t max >	Total (s) Series (s) Residual BPL (s) #RES
	5	177 9.604 • 10 -4	27.51	7.13	17.66	3624
	10	95	1.789 • 10 -3	23.01	9.71	12.24	2285
	15	79	2.152 • 10 -3	26.43	14.23	11.31	1906
	20	85	2.000 • 10 -3	38.14	23.89	13.29	2045
	25	48	3.542 • 10 -3	28.97	19.13	8.12	1157

Remark 3.3.2 We also investigated the CPU time required to compute one-BPL solution for various truncation orders N , for this example. The results have shown that calculating one BPL solution is not time-consuming. We found that the CPU time required to compute one-BPL solution for N = 10 is approximately 1.3•10 -2 s and increases with the truncation order. The CPU time required to compute the Borel transforms, the Padé coefficient, the Gauss Laguerre quadrature step, and the BPL solution were measured. These four steps were proven to be extremely fast. The BPL derivative is also quickly evaluated. And, if one wants to know which step takes the most CPU time, one can say it is the step related to the evaluation of the Padé coefficients.

4. IFS:

Table 3 .

 3

10 -Statistics to reach t = 0.17s for the Burgers equations 2D, using IFS with different truncation orders N , for a fixed residual tolerance = 7.3 • 10 -3 and a viscosity ν = 0.02. The grid used is 20 × 20.

Table 3 .

 3 [START_REF] Berrahma-Chekroun | Resolution of the transient dynamic problem with arbitrary loading using the asymptotic method[END_REF] -Statistics to reach t = 2s for the Navier-Stokes equations on an L-shaped domain, using ANM with different truncation orders N , for a fixed momentum residual tolerance M = 10 -4 and a continuity residual tolerance c = 9.898 • 10 -2 . The viscosity is fixed to ν = 0.02.

	< t max >	Total (s) Series (s) Residual Pol (s) #RES
	5 4084 4.897 • 10 -4	406.01	123.85	271.59	133772
	10 4125 4.848 • 10 -4	574.25	301.17	261.79	119620
	15 4132 4.484 • 10 -4	840.57	543.44	285.22	119803
	20 4134 4.837 • 10 -4 1091.62	801.37	278.50	119786
	25 4134 4.837 • 10 -4 1474.29	1170.54	292.23	119786
					#RES
	5 1604 1.24 • 10 -3	162.57	47.47	110.77	32186
	10 1604 1.24 • 10 -3	253.15	115.45	133.40	32186
	15 1604 1.24 • 10 -3	358.53	201.10	154.10	32186
	20 1604 1.24 • 10 -3	495.32	313.53	177.46	32186
	25 1604 1.24 • 10 -3	625.85	420.84	200.01	32186

N Steps < t max > Total (s) Series (s) Residual BPL (s)

Table 3 .

 3 [START_REF] Braikat | A High Order Continuation Based On Time Power Series Expansion And Time Rational Representation For Solving Nonlinear Structural Dynamic Problems[END_REF] -Statistics to reach t = 2s for the Navier-Stokes equations on an L-shaped domain, using BPL with different truncation orders N , for a fixed momentum residual tolerance M = 10 -4 and a continuity residual tolerance c = 9.898 • 10 -2 . The viscosity is fixed to ν = 0.02.

	N Steps < t max > Total (s) Series (s) Residual IFS (s) #RES
	5 5233 3.82 • 10 -4 2760.86	348.22	2406.64	104840

Table

Table 3 .

 3 [START_REF] Daya | A numerical method for nonlinear eigenvalue problems application to vibrations of viscoelastic structures[END_REF] -The total CPU time required for the selected time interval [400 ; 405]s, for the LID driven cavity problem, using ANM, BPL and IFS, for N = 6. The Reynolds number is Re = 12500. The grid used is 64 × 64. The discrete-time step of RK1 used is ∆t = 10 -4 s. The computation study of ANM, BPL and IFS are only for a forced t max and are just presented as an indication.

	Total(s) RK1 Total(s) ANM Total(s) BPL Total(s) IFS
	32.942	612.271	179.430	890.484

  Solution norm ||φ(x, y, t)|| 2 computed by time perturbation-resummation methods ANM, BPL and IFS, and the classic scheme TVDRK4 with a discrete-time step ∆t = 10 -3 s. The grid size is 128 × 128.

		44.743 44.7435 44.744						TVDRK4 ANM BPL IFS
	SOLUTION NORM	44.741 44.7415 44.742 44.7425
		44.7405
		44.74
		44.7395	0	0.5	1	1.5	2	2.5	3
						TIME
	Figure 4.4 -

Table 4 .

 4 1 -Number of steps needed to reach t = 3s using ANM, BPL, and IFS for a truncation order N = 15, for different values of the residual tolerance . The grid used is 128 × 128. TVDRK4 with ∆t = 10 -3 s needs 3000 steps to reach the same final time.

	Residual tolerance ANM BPL IFS
	10 -2	108	81	103
	10 -3	109	86	105
	10 -4	110	90	109
	10 -5	114	94	116

Étude numérique de méthodes d'intégration temporelle basées sur des techniques de perturbation. : application à des problèmes de mécanique Cynthia Tayeh 2022 Étude numérique de méthodes d'intégration temporelle basées sur des techniques de perturbation. : application à des problèmes de mécanique Cynthia Tayeh 2022

Étude numérique de méthodes d'intégration temporelle basées sur des techniques de perturbation. : application à des problèmes de mécanique Cynthia Tayeh 2022

Étude numérique de méthodes d'intégration temporelle basées sur des techniques de perturbation. : application à des problèmes de mécanique Cynthia Tayeh

Towards application to Fluid-Structure Interaction problems

The recursive relation that permits to obtain the series terms φ i is obtained:

Recurrence formulas for time-dependent velocity field

Now, we consider the case where the velocity field u (x, t) is spatial time-dependent, as result it is necessary to look for a time polynomial approximation also for it. Therefore, consider the following two power series:

By substituting these power series into the transport equation Eq.(4.2) and equating the respective power of t, one can obtain the following recurrence formula that permits to calculate φ i series terms:

Therefore, series terms are given by the following recurrence formula:

We recall that u i and φ i , ∀i ≥ 1 are space dependent. However, for the sake of simplicity of writing, we did not use φ i (x, t) and u i (x) in this system of equations.

The terms u i and φ i are the unknowns of this recurrence formula. The velocity coefficients

∀i ≥ 1 should be calculated accordingly to the velocity data of each problem, in order to determine φ i terms (more details in the numerical example). Now, we seek to test numerically the three ways of reinitialization, in the context of solutions that are continuous in time. The latter are proposed in subsection 4.2.3 and recalled here:

1. The classical PDE reinitialization Eq.(4.5) only after the last step of the continuation procedure of ANM i.e. at the final physical time t, 2. The classical PDE reinitialization Eq.(4.5) at the end of each continuation step of ANM, i.e end of each t max , 3. Time perturbation of the reinitialization procedure without using a modified version of Eq.(4.5). We are looking for a way of these three proposals that yields to well-positioned iso-contours like those illustrated in Fig. (4.8b).

ANM with the classical PDE reinitialization at the end of the last continuation step:

Here, the classical PDE reinitialization (see Eq.(4.5)) only after the final step of the continuation procedure of ANM (i.e end of the final time t = 3s) is tested. As can be seen in Fig. (4.9), this does not lead to proper level lines as those presented in Fig. (4.8b). Thus, this test of reinitialization is not suitable.

APPENDIX A

Derivatives formulas for Borel-Padé-Laplace

The Borel sum denoted by Su(t) in a direction d = R + is given by:

As already mentioned in Chapter 1, we will note:

-S(t) = Su(t) the Borel sum.

-Ṡ(t) = dSu(t) dt the first derivative of Borel sum.

-S(t) = d 2 Su(t) dt 2 the second derivative of Borel sum. The first derivative of Borel sum is given by:

The change of variable ξ = rt, yields to:

Finally, the first derivative of Su(t) is

The second derivative of Borel sum is given by:

with

and

For I 1 :

Let the change of variable ξ = rt

Therefore

For I 2 :

Let the change of variable ξ = rt

Therefore, the second derivative of Borel sum is

Remark A.0.1 Note that the evaluation of the derivatives of the Borel sum must start from t 0 = 0. For this reason, for all the examples we have treated in this manuscript, we have numerically fixed the initial time with a small value very close to 0 (for example

Therefore, as a summary:

A brief recall of Gauss-Laguerre quadrature is presented: In numerical analysis Gauss-Laguerre quadrature is an extension of the Gaussian quadrature method for approximating the value of integrals of the following kind:

by the following sum

with w i are the quadrature coefficients (or weights). The x i points, or nodes are real, distinct, unique and are the roots of N G -Laguerre polynomials. Generally, one can also consider the integrals that have a singularity of the law of power

x α known to x = 0, for α > -1 with α a real number, that lead to integrals of forms:

This allows to estimate such integrals accurately for f (x) smooth or polynomial even in the case where α is not an integer.

It should be noticed that the nodes x i and the weighs w i depend on α and N G . Here are two different strategies for calculating first and second derivatives of BPL.

Strategy 1

Strategy 1 is an application of Gauss-Laguerre formula or the generalized formula with α = 0. Borel sum is given by:

The first derivative of Borel sum is given by:

The second derivative of Borel sum is given by:

Therefore, with α = 0, we have

For this first way, the Padé solution P (r i t) is calculated only once for the estimation of the Borel solution, the first and second derivatives. Note that we used this strategy in all the test cases in the manuscript.

Strategy 2

Strategy 2 is the second way is based on generalized Gauss-Laguerre. For α = 0, we have:

For α = 1, we have:

For α = 2, we have:

Therefore,

For this strategy, the development of Padé is not the same for Borel sum and its derivative since the roots and the weight depend on α. As a result, this strategy is more timeconsuming than the first one. However, it has been presented for the integrals that have a singularity of the law of power x α known to x = 0, for α > -1. It was not used in any of the test cases presented.

APPENDIX B

Derivatives formulas for Inverse factorial series

In this Appendix, we develop the first and second derivatives of Inverse factorial series I(t) of the second way of the subsection 1.3 (large order N ) in Chapter 1. We recall some definition of this subsection before starting to evaluate the derivatives. Suppose that u(t) is 1-summable in a direction d, with an angle θ with the positive half-axis. Let τ l = s l e iθ with l = (1, 2, . . . , N ), where {s 1 , s 2 , . . . , s N } is any sequence of complex numbers. Let

Note that if the series are summable in R + direction, so y = z and we will choose

We call Inverse factorial series the following :

where v n is the n-th term of IFS which is calculated using the following recursive algorithm based on

The terms v

2 , . . . , v (1) n , v

n+1 are the terms of the IFS denoted by: v 1 , v 2 , . . . , v n , v n+1 .

First derivative

The first derivative of Inverse factorial series is given by:

We will build a recursive algorithm of derivative to evaluate İ(t).

We have

This algorithm for n ≥ 1 is based on

with τ 0 = 0 and dv If the series are summable in the R + direction, so ẏ = ż (B.12)

The terms dv

, . . . , dv (1) n dt , dv

(1) n+1 dt . . . are the terms of the D IF S , (first derivative terms of IFS) denoted by v1 , v2 , . . . , vn , vn+1 .

Second derivative

The second derivative of Inverse factorial series is given by:

We will build a recursive algorithm of second derivative. Since the second derivative is complicated to write it directly, we will evaluate it using A, B, C such that

and and

The derivative of Eq.(B.9) is given by

APPENDIX C

Derivatives formulas for Meijer-G sum

The Meijer-G sum denoted by S M G (t) is given by:

Concerning the derivatives of arbitrary order h of the Meijer-G sum, first one needs to recall the following relationship for Meijer-G function.

We have ∀z = 0:

where w is a constant.

First derivative

Thus, the first derivative (order of derivative h = 1) of Meijer-G sum is:

Second derivative

The second derivative (order of derivative h = 2) of Meijer-G sum is:

Derivative of order h

The derivative of order h of Meijer-G sum is:

Appendix D

APPENDIX D

Velocity series terms for the Lid-driven cavity example at t 0 = 400s for a Reynolds number fixed to Re = 12500

The behavior of the velocity series terms u for the Lid-driven cavity example, for an initial time t 0 = 400s, a Reynolds number fixed to Re = 12500 and a truncation order that varies from N = 0 to N = 35 are shown in