Nabila -Hassan Nagid

Ridouan Belhadj -Mohamed

Amattouch
email: amattouch36@yahoo.fr

Accelerated Restricted Additive Schwarz Method for Nonlinear Partial Differential Equations

Keywords: nonlinear PDEs, domain decomposition, restricted additive Schwarz (RAS), additive Schwarz (AS), vector epsilon algorithm, Bratu problem EDPs non linéaires, décomposition de domaine, Schwarz additif restreint epsilon algorithme vectoriel, problème de Bratu

In this paper, the Restricted additive Schwarz (RAS) method is applied to solve a nonlinear partial differential equations (PDEs). To accelerate the RAS iterations, we propose to apply vector ε-algorithm. Some convergence analysis of the proposed method is presented, and applied succeffully to Bratu problem. The obtained results show the efficiency of the proposed approach. Moreover, the algorithm yields much faster convergence than the classical Schwarz iterations.

RÉSUMÉ.

Dans ce papier, on applique la méthode de Schwarz additif restreint (RAS), pour résoudre un problème aux limites non linéaire. On applique ensuite une méthode d'accélération des itérations de Schwarz basée sur l'epsilon algorithme vectoriel (ε-algorithme). Nous présentons par la suite une analyse de la convergence de la méthode que nous validons sur l'exemple du problème non linéaire de Bratu. Les résultats obtenus montrent l'efficacité de l'approche proposée, de plus ces résultats montrent également que l'algorithme converge plus rapidement que la méthode de Schwarz classique.

Introduction

In scientific computing, the domain decomposition methods are now commonly used when solving large linear or nonlinear systems arising from the discretization of partial differential equations PDEs [START_REF] Quarteroni | Domain decomposition methods for PDEs[END_REF][START_REF] Smith | Domain decomposition: parallel multilevel methods for elliptic partial differential equations[END_REF][START_REF] Jacques | Introduction au méthodes de décomposition de domaine[END_REF]. The first models of these methods have been established by H.A.Schwarz, the idea is to decompose a large problem into a series of smaller sub-problems, and therefore more easily resolved. There are several variant of the Schwarz method, for example additive Schwarz method (AS), and the restricted additive schwarz method (RAS) [START_REF] Formmer | An algebraic convergence theory for restricted additive schwarz methods using weighted max norms[END_REF][START_REF] Cai | A restricted additive Schwarz preconditioner for general sparse linear systems[END_REF][START_REF] Efstathiou | Why restricted additive Schwarz converges faster than additive Schwarz[END_REF][START_REF] Nabben | Convergence theory of restricted multiplicative Schwarz methods[END_REF].

In numerical analysis, the ε-algorithm is a nonlinear algorithm for accelerating the convergence of numerical sequences. This algorithm was proposed by Peter Wynn to calculate the Shanks transformation [START_REF] Brezinski | Variations on Richardson's method and acceleration[END_REF][START_REF] Brezinski | Accélération de la convergence en analyse numérique[END_REF]. There are different variants of the ε-algorithm that can be used with vector sequences (vector ε-algorithm, topological or scalar applied to each component of the vector) [START_REF] Salam | On the vector εalgorithm for solving linear systems of equations, Numerical Algorithms[END_REF][START_REF] Brezinski | Méthodes d'accélération de la convergence en analyse numérique[END_REF][START_REF] Brezinski | Variations on Richardson's method and acceleration[END_REF][START_REF] Brezinski | Accélération de la convergence en analyse numérique[END_REF][START_REF] Roland | Méthodes d'accélération de convergence en analyse numérique et en statistique[END_REF].

In litterature, there exist several works in which we apply the acceleration method to Schwarz iterations, for example, in [START_REF] Duminil | Nonlinear Schwarz iterations with reduced rank extrapolation[END_REF], the authors propose an algorithm based on polynomial methods especially the reduced rank extrapolation method (RRE), in [START_REF] Linel | Méthodes de décomposition de domaine en temps et en espace pour la résolution de systèmes d'EDOs non-linéaires[END_REF], the authors accelerate Schwarz iterations for the ordinary differential equations ODEs, using epsilon algorithm. There exist many other authors that have proposed different ideas for accelerating domain decomposition methods [START_REF] Lube | Acceleration of a Domain Decomposition Method for Advection-Diffusion Problems[END_REF]24,[START_REF] Garbey | Acceleration of the Schwarz method for elliptic problems[END_REF][START_REF] Garbey | On some Aitken-like acceleration of the Schwarz method[END_REF][START_REF] Linel | TROMEUR-DERVOUT " Aitken-schwarz and schur complement methods for time domain decomposition[END_REF]. The purpose of this paper is to accelerate the nonlinear iterative Schwarz, using the vector ε-algorithm for PDEs. In particular we treat the nonlinear reaction diffusion problems. We apply this algorithm to the sequences of vectors produced by AS, and RAS methods, and we show experimentally that the proposed algorithm can provide faster convergence measured both in number of iterations and in CPU Times

Iterative Schwarz Methods

Linear Schwarz iterations

We consider the following problem

L(u) = f in Ω Bu = g on ∂Ω [1]
where L is a linear operator, B is a boundary operator and Ω is a bounded domain of

R d (d = 1, 2, ..).
H.A.Schwarz proposed an iterative method for the solution of classical boundary value problems. There are several variants of Schwarz algorithms, additive, multiplicative, and several hybrid types, a number of them are discussed in detail in [START_REF] Dryja | Schwarz analysis of iterative substructuring algorithms for elliptic problems in three dimensions[END_REF][START_REF] Lottes | Hybrid multigrid/Schwarz algorithms for the spectral element method[END_REF].

In the present work, we have considered the additive Schwarz method .

Let consider these notations, Ω as a union of nonoverlapping domains Ω j , j = 1, .., p and

Γ j = ∂Ω j ∩ ∂Ω, Γ ij = ∂Ω i ∩ Ω j . τ is the Richardson parameter (0 < τ ≤ 1/p).
The additive Schwarz algorithm in the Richardson version is written as follows:

F or n = 0, ...

F or each j = 1, ..., p solve    L(v n+1,j ) = f in Ω j Bv n+1,j = g on ∂Γ j v n+1,j = u n on ∂Ω j \Γ j Compute w n+1 = v n+1,1 + ... + v n+1,p U pdate u n+1 = (1 -pτ )u n + τ w n+1
The discretization of problem 1 leads to a linear system of equations of the form

Au = f [2]
where A is the discretization matrix by a numerical methods (Finite element, Finite Difference , or finite volume). We use the same notation f after discretization. a stationary iterative method for 2 is given by

u n+1 = u n + M -1 (f -Au n ) [3]
with a given initial approximation u 0 to the solution of 2. Algebraic domain decomposition methods group the unknowns into subsets, u j = R j u, j = 1, ..., p, where R j are rectangular restriction matrices. Coefficient matrices for subdomain problems are defined by A j = R j AR T j . The additive Schwarz (AS) preconditioner , and the restricted additive Schwarz (RAS) preconditioner are defined by:

M -1 AS = p j=1 R T j A -1 j R j , C M -1 RAS = p j=1 RT j A -1 j R j [4]
where the RT j correspond to a non-overlapping decomposition, and it consists of zeroes and ones, in such a way that p j=1

RT j R j = I.
The additive Schwarz method constructs the sequence of approximations {u n } n∈N by setting :

u n+1 = u n + p j=1 R T j A -1 j R j (f -Au n ) n = 0, 1, ... [5]
The restricted additive Schwarz (RAS) algorithm is given by :

u n+1 = u n + p j=1 RT j A -1 j R j (f -Au n ) n = 0, 1, ... [6]

Nonlinear Schwarz iterations

Consider now the problem 1 with a nonlinear operator L. After discretization, we obtain an algebric non linear system

F (u) = 0 [7]
we transform this problem to a fixed point form

G(u) = u [8]
where F and G are two mappings from R n → R n . using the same notation as before, we define G on each subdomain Ω j , j = 1, 2, ..., p, as follows :

G j (X) = R j G( R T j (X)) [9]
the corresponding nonlinear restricted Schwarz method is defined by

u n+1 = u n + p j=1 RT j G j (R j (u n ))n = 0, 1, ..., [10]
we also consider for the solution of 5 the Schwarz-Newton methods, where in each subdomain, the non linear problem is solved by a Newton, see [START_REF] Oana | Résolution des équations d'Euler[END_REF][START_REF] Cai | Inexact Newton methods with restricted additive Schwarz based nonlinear elimination for problems with high local nonlinearity[END_REF].

The Vector ε-Algorithm

Definitions and properties

The εalgorithm is a nonlinear extrapolation method for accelerating the convergence of sequences, one can say also that this is a generalization of Aitken method [START_REF] Brezinski | Méthodes d'accélération de la convergence en analyse numérique[END_REF][START_REF] Brezinski | Accélération de la convergence en analyse numérique[END_REF]. There exist several versions of the epsilon algorithm( topological, scalar, and vector epsilon algorithm) [START_REF] Salam | On the vector εalgorithm for solving linear systems of equations, Numerical Algorithms[END_REF][START_REF] Brezinski | Méthodes d'accélération de la convergence en analyse numérique[END_REF][START_REF] Brezinski | Accélération de la convergence en analyse numérique[END_REF].

In this work, we are only interested in the vector form. We consider thereafter the fundamental algebraic results in the theory of the vector ε-algorithm.

Before that, we recall some results concerning the Aitken's process.

Definition 3.1 Let U = (u n ) n∈N be a sequence of real numbers, the ∆ 2 Aitken process consists to transform the sequence

(u n ) in a new sequence ε (n)
2 defined by :

ε n 2 = u n+2 u n -u 2 n+1 ∆ 2 u n Theorem 3.1 If the sequence U = (u n ) n∈N satisfies the condition lim n→∞ (u n+1 -u)/(u n -u) = lim n→∞ ∆u n+1 /∆u n = ρ = 1
then the sequence ε n 2 converges to u faster than u n+1 . Proof. see for example [START_REF] Brezinski | Accélération de la convergence en analyse numérique[END_REF].

Writing ε

(n)
2 based on determinants:

ε (n) 2 = u n u n+1 ∆u n ∆u n+1 1 1 ∆u n ∆u n+1
we seek the conditions on (u n ) that verifies ε

(n) 2
= u for n > N (N is a given rank), therefore, we have the following theorems Theorem 3.2 A necessary and sufficient condition to have ε

(n) 2 = u ∀n > N , is that the sequence (u n ) verifies a 0 (u n -u) + a 1 (u n+1 -u) = 0 ∀n > N ,with a 0 + a 1 = 0
This theorem can be generalized to high order using a non-linear acceleration method, the Shanks transformation [START_REF] Salam | On the vector εalgorithm for solving linear systems of equations, Numerical Algorithms[END_REF][START_REF] Brezinski | Accélération de la convergence en analyse numérique[END_REF].

This transformation called e n k (U ) is built such that e n k (U ) = u ∀n > N, and it consists in computing the quantities e n k (U ) as follows

e n k (U ) = k i=0 a (n,k) i u n+i ∀n > N with k i=0 a (n,k) i = 1
from these equations, it is easy to obtain a determinantal formula for e n k (U )

e n k (U ) = u n ... u n+k ∆u n ... ∆u n+k . . . . . . ∆u n+k-1 ∆u n+2k-1 1 ... 1 ∆u n ... ∆u n+k . . . . . . ∆u n+k-1 ∆u n+2k-1 with k i=0 a (n,k) i = 1
where ∆ is the difference operator ∆u n = u n+1 -u n .

the previous results leads to the following theorem.

Theorem 3.3 If for a fixed k, the sequence U is such that there exists u ∈ R and a 0 , ..., a k ∈ R with

k i=0 a i = 0 satisfying k i=0 a i (u n+i -u) = 0 ∀n > N then e n k (U ) = ε (n) 2k = u ∀n > N.
The proof of Theorem 3.2 and Theorem 3.3 are given for example in [START_REF] Salam | On the vector εalgorithm for solving linear systems of equations, Numerical Algorithms[END_REF][START_REF] Brezinski | Accélération de la convergence en analyse numérique[END_REF][START_REF] Roland | Méthodes d'accélération de convergence en analyse numérique et en statistique[END_REF].

Remark 1 A recursive rule for computing the quantities e n k (U ) of shanks transformation has been given by Wynn [START_REF] Brezinski | Variations on Richardson's method and acceleration[END_REF][START_REF] Brezinski | Accélération de la convergence en analyse numérique[END_REF]. These quantities can be computed by the following ε-algorithm:

ε (n) -1 = 0 ε (n) 0 = U n n = 0, 1, ... ε (n) k+1 = ε (n+1) k-1 + (∆ε (n) k ) -1
n, k = 0, 1, ... where the inverse of a vector y is defined by:

y -1 = y y 2 2
.

using the theorem 3.3, Gekeler [START_REF] Gekeler | On the solution of systems of equations by the epsilon algorithm of Wynn[END_REF] has proved that vector epsilon algorithm provides a direct method for solving the linear systems of equations.

Theorem 3.4 If we apply the vector ε-algorithm to the sequence {u n } produced by

u n+1 = Au n + b
with a given u 0 and A is a real square matrix such that I -A is invertible, then we have

ε (n) 2m = u f or n = 0, 1, ...
where u = (I -A) -1 b and m is the degree of minimal polynomial of A for the vector u 0 -u.

Proof Let p(t) = m i=0 a i t i the minimal polynomial of A for the vector u 0 -u, the definition of the minimal polynomial of a matrix for a vector is

m i=0 (a i A i ) (u 0 -u) = 0 the matrix I -A is invertible, therefore p(1) = m i=0 a i = 0 on the other hand, we have u = Au + b so u n+1 -u = A(u n -u) and u k -u = A k (u 0 -u) ∀k0
replacing in p(t) we have

A n m i=0 a i (u i -u) = m i=0 a i (u n+i -u) = 0 ∀n so using Theorem 3.3 we prove that ε (n) 2m = u ∀n ≥ 0

Vector ε-algorithm applied to AS/RAS for linear systems

We consider the following problem :

L(u) = f inΩ Bu = g on∂Ω [11]
In the case where the operator L is linear, a discretization of the equation 11 leads to a linear system of equations of the form

Au = f [12]
The additive Schwarz methods allows to compute the sequence of approximations {u n } n∈N by setting :

u n+1 = u n + M -1 AS (f -Au n ) [13] with M -1 AS = p j=1 R T j A -1 j R j i.e.,         u n+1,1 u n+1,2 . . . u n+1,p         =         u n,1 u n,2 . . . u n,p         +         A 11 A 22 0 . 0 . . A pp         -1         r n,1 r n,2 . . . r n,p        
we can write

u n+1 = B u n + F [14]
where

B = I -M -1 AS A Theorem 3.5 Suppose A and M -1
AS a real square matrices that have the same size such that C = M -1 AS A is non singular. If we apply the vector ε-algorithm to the adiitive Schwarz sequence 14, then ε n 2m = u, where u is the solution of the linear system Cu = F, where m the degree of the minimal polynomial of B .

Proof The solution u is a fixed point of the operator

u → u + M -1 AS (f -Au)
Let P = M AS -A is the difference between A and M AS when 13 converge, it converges to the solution of the preconditioned system

M -1 AS Au = M -1 AS f
by setting y n = u n -u we obtain

u n+1 = u n + M -1 AS (f -Au n ) = (I -M -1 AS (M AS -P ))u n + M -1 AS f = M -1 AS P u n + M -1 AS Au = M -1 AS P u n + M -1 AS (M AS -P )u = u + M -1 AS P (u n -u)
the equivalent system becomes y n+1 = M -1 AS P y n = By n using theorem 3.4, we have ε n 2m = u,where m is the degree of the minimal polynomial of B and we have

P m (B)(u 0 -u) = k n=0 γ n B n (u 0 -u) = k n=0 γ n (u n -u) = 0
where γ n are the coefficients of the polynomial P d such that P d (1) = 1.

Vector ε-algorithm applied to AS/RAS for nonlinear systems

We consider now the nonlinear reaction diffusion problem defined by:

Lu -f (u) = g in Ω [15]
We can write Lu = G(u) in the sense of D (Ω) [START_REF] Dinca | A fixed point method for the p(.)-Laplacien[END_REF] The corresponding discretized problem can be written as follows:

AU = G(U ) ∈ R p [17]
where A is the matrix of the discretized operator L, obtained by the finite elements on a regular grid, G : R p -→ R p is the nonlinear function, and U is the vector containing the approximation of the solution of the continuous problem to grid points. We remark , that if we put F = A -1 G(.), then the problem 17 is equivalent to U = F (U ).

Let solve the following problem:

f ind x ∈ R p such that x = F (x) [18]
where F : R p -→ R p is differentiable in the sense of Frechet, in a neighborhood of x, m is the degree of minimal polynomial of F (x) for the vector x n -x and r is the multiplicity of the root (λ = 0) for this minimal polynomial.

Knowing x 0 we set u 0 = x n and we solve for k = 1, ..., 2m-r the following iterative problem u k = F (u k-1 ).

To calculate ε (r) 2(m-r) we apply the ε-algorithm to the vectors u 0 , ..., u 2m-r . then we take

x n+1 = ε (r) 2(m-r) .
The application of the epsilon algorithm to the nonlinear RAS provides a method of resolution with quadratic convergence, see [START_REF] Brezinski | Accélération de la convergence en analyse numérique[END_REF].

Theorem 3.6 Let u n+1 = u n + p j=1 RT j F j (R j (u n )) n = 0, 1, ..., [19]
where F : R p -→ R p is defined on each subdomain

F j (X) = R j F ( R T j (X)) [20]
and

G(u) := F(u) -u = p j=1 RT j F j (u) -u = 0 [21]
If F is differentiable in the sense of Frechet in a neighborhood of u, and I -F (u) is invertible, then there exists a neighborhood V of u such that

∀ x 0 ∈ V x n+1 -u = o( x n -u 2 ) n = 0, 1, ...
Proof If F is differentiable in the sense of Frechet in a neighborhood of u, we have :

u k+1 -u = F (u)(u k -u) + o( u k -u 2 )
Let p(t) = m i=0 a i t i the minimal polynomial of F (u) for the vector u n -u.

Since I -F (u) is invertible, we have p(1) = m i=0 a i = 0. We have :

u 1 -u = F (u)(u 0 -u) + o( u 0 -u 2 ) and u k -u = [F (u)] k (u 0 -u) + o( u 0 -u 2 ) Replacing in the minimal polynomial m i=0 a i [F (u)] i (x n -u) = m i=0 a i (u i -u) + o( u 0 -u 2 ) = 0 therefore u 0 = x n ,

The algorithm (Epsilon-RAS)

In case of convergence , lim n-→∞ u n = u.

1) choose a starting approximation x 0.

2) set u 0 = x n at the iteration n, and

u k+1 = u k + p j=1 k = 0, ....2m -r.
3) apply the epsilon algorithm to the vectors u 0 , ..., u 2m-r to calculate ε (r) 2(m-r) . 4) compute x n+1 such that

x n+1 = ε (r) 2(m-r) (r = 0 if the partial Frechet derivative of F is invertible)

Numerical Experiments

In this section, we compare the performance of Schwarz iterations with those accelerated with the vector epsilon algorithm in terms of number of iterations and CPU time.

We treat two different applications, the first one in the linear case and the second one in the nonlinear case.

We have implemented the finite element discretization in two spatial dimensions and use the software freefem++ to compute numerical results.

Application to the Bratu problem

We consider now the following nonlinear reaction diffusion problem

-∆u + λe u = f [22]
The domain is the unit square Ω = [0, 1] × [0, 1] decomposed uniformly into p non overlapping subdomains.

Using a finite element discretization, we obtain the following nonlinear system of equations

AX + λe X -b = 0 in Ω [23]
Using FreeFem++, the solution of the problem on all the domain Ω is presented in Fig .1 Let p the number of subdomains, Fig. 2 shows the solution for p=4. The following results reported in Table 1 and in Fig. 3, show the L ∞ Error between the approximate solution and the exact solution of the nonlinear problem for each algorithm. The RAS algorithm reaches the precision 10 -7 for p=36, while the Epsilon-RAS algorithm reaches the same precision only for p=25.

When we compare CPU-Time, it can be observed that additive nonlinear Schwarz and restricted nonlinear Schwarz take too long to converge than their accelerated forms. This shows that the Epsilon-algorithm performs very well for accelerate Schwarz method. The results are reported in Table2 and in Fig. 4 Subdomains(P ) CPU-TimeAS CPU-TimeEpsilon -AS CPU-TimeRAS CPU-

Conclusion

We have proposed an accelerated form of Schwarz iterations for nonlinear problems (AS-RAS) using the vector epsilon algorithm .

Comparing CPU-Time and the error, we show that this accelerated method is fast and it has a better accuracy than the direct classical Schwarz method.

As perspective of the present work, we can generalize the acceleration method for non stationary PDEs, and apply it to a real modelling case.

  using theorems 3.2 and 3.3 we get ε (r) 2(m-r) = u + o( x n -u 2 ) .

Fig 1 :

 1 Fig 1 : Bratu problem: Solution on allΩ

Fig 2 :

 2 Fig 2 : Bratu problem: Solution for p=4

  Subdomains(P ) ErrorAS ErrorEpsilon -AS ErrorRAS ErrorEpsilon -RAS P = 4 0.8542 × 10 -5 0.6899 × 10 -5 0.7974 × 10 -5 0.6123 × 10 -5 P = 16 0.5365 × 10 -5 0.9410 × 10 -6 0.4123 × 10 -5 0.7745 × 10 -6 P = 25 0.4036 × 10 -5 0.5123 × 10 -7 0.7741 × 10 -6 0.3012 × 10 -7 P = 36 0.9023 × 10 -7 -0.7841 × 10 -7 -

Fig 3 :

 3 Fig 3 : Bratu problem: Convergence for different number of subdomains p

Fig 4 :

 4 Fig 4 : Bratu problem: Convergence for different number of subdomains p

Table 1 .

 1 Bratu problem: Error norm versus p

Table 2 .

 2 Bratu problem. CPU-Time versus p

	TimeEpsilon -RAS


H Biyadi -M

Er-Riani -K Chetehouna -M

Jarroudi El
email: eljarroudi@hotmail.com

Finite element modeling of exothermic reaction with convective boundary conditions Sous-titre

Keywords: Finite element method, exothermic reaction, Biot number, Arrhenius kinetics. MOTS-CLÉS : méthode des éléments finis, réaction exothermique, nombre de Biot, la cinétique d'Arrhenius

In this work, we used a finite element method for the boundary value problem governing the ignition of a solid reactant undergoing slow oxidation in an arbitrary domain. We also investigate the effect of various parameters on the temperature distribution profiles.

RÉSUMÉ. Dans ce travail, nous avons utilisé une méthode des éléments finis pour le problème de la valeur limite régissant l'allumage d'un réactif solide qui subit une oxydation lente dans un domaine quelconque. Nous étudions aussi l'effet de différents paramètres sur les profils de distribution de température.

Introduction

The problem of self-ignition of energetic materials is of great importance in the field of safe use and storage of explosives and pyrotechnics. The self -ignition occurs when an explosive substance is heated to a temperature sufficient for the decomposition process begins to produce significant exothermic effects. Then there is a release of heat, accompanied by a rise in temperature and velocity of thermal decomposition. This causes a thermal runaway generating an explosion that can sometimes produce a devastating detonation wave. These phenomena were first studied analytically and Frank-Kamenetskii [START_REF] El | Diffusion and heat transfer in chemical kinetics[END_REF] has developed the quasi-stationary theory of thermal explosion for determining the critical conditions for self-ignition of explosives.

Many problems of practical importance in ignition theory involve the definition and determination of either the critical ambient temperature or the critical initial temperature. These phenomena were first subject to analytical studies and Frank-Kamenetskii [I] developed the theory of quasi-stationary thermal explosion, which identifies the critical conditions for self-ignition of explosive substances.

After the pioneer work of Frank-Kamenetskii [START_REF] El | Diffusion and heat transfer in chemical kinetics[END_REF], various aspects of these phenomena were studied intensively by many authors, including [START_REF] Boddington | Thermal theory of spontaneous ignition: Criticality in bodies of arbitary shape[END_REF][START_REF] Boddington | Criteria for thermal explosions with and without reactant consumption[END_REF][START_REF] Balakrishnan | Critical values for some nonclass A geometries in thermal ignition theory[END_REF][START_REF] Makinde | Strongly exothermic explosions in a cylindrical pipe:a case study of series summation technique[END_REF][START_REF] Chetehouna | On the Critical behaviour of exothermic explosions in class A geometries[END_REF]. The book [START_REF] Bebernes | Mathematical problems from combustion theory[END_REF] contains a part of these investigations and a large bibliography.

When reactant consumption is neglected, the equation for the heat balance in a bounded domain Ω, can be written as:

k∆T + C 0 QA exp(- E RT ) = c ∂T ∂t in Ω [1]
with the boundary condition

k ∂T ∂n + H(T -T a ) = 0 on ∂Ω [2]
Where k is the thermal conductivity, c is the thermal capacity, C 0 is the initial concentration of the reactant species assumed constant, Q is the heat of reaction.

According to the Arrhenius law, the dependence of the rate constant on the absolute temperature has the form A exp(-E RT ) , where E is the activation energy, R the universal gas constant. T a is the ambient temperature, H is the heat transfer coefficient. ∂T ∂n is the normal derivative of the temperature, n is the outward normal vector on the boundary . Following Balakrishnan et al [START_REF] Balakrishnan | Critical values for some nonclass A geometries in thermal ignition theory[END_REF], we introduce the following dimensionless variables in Eq. 1:

u = ( E RT 2 a )(T -T a ), ε =( RTa E ), Bi = Ha0 k , δ = QσAEa 2 0 kRT 2 a exp(-E RTa ).
and obtain the dimensionless governing equation together with the corresponding boundary conditions as:

∆u + δ exp( u 1 + εu ) = 0 on Ω [3] ∂u ∂n + Biu = 0 on ∂Ω [4]
δ is the Frank-Kamenetskii parameter, Bi is the Biot number and a 0 is some characteristic length scale of the problem.

One of the commonly used formulation is the exponential approximation of Frank-Kamenetskii, in which ε is supposed so small to be neglected. So the problem (3-4) is reduced to the following nonlinear eigenvalue problem:

∆u + δ exp(u) = 0 x ∈ Ω ∂u ∂n + Biu = 0 x ∈ ∂Ω [5]
The problem ( 5) have been studied both analytically and numerically, for simple class A geometries (infinite cylinder, infinite slab, and a sphere) and also for some non class A geometries [START_REF] Balakrishnan | Critical values for some nonclass A geometries in thermal ignition theory[END_REF]. Thermal ignition in all these geometries can be formulated using a single variable [START_REF] El | Diffusion and heat transfer in chemical kinetics[END_REF][START_REF] Balakrishnan | Critical values for some nonclass A geometries in thermal ignition theory[END_REF][START_REF] Chambré | On the solution of the Poisson-Boltzmann equation with application to the theory of thermal explosions[END_REF][START_REF] Scatter | Path following, critical points and shape factors[END_REF].

Thermal stability analysis in the theory of thermal ignition, with or without consideration of consumption and diffusion of the reactant, has generally been performed using the Frank-Kamenetskii number as a parameter of bifurcation, ε is kept fixed. But this number dissimulates the role of the most practically significant control variable, namely the ambient temperature. Several authors [START_REF] Burnell | Determination of critical ambient temperatures for thermal ignition[END_REF][START_REF] Gray | Critical initial conditions for thermal ignition[END_REF][START_REF] Gray | Critical initial conditions for spontaneous thermal ignition[END_REF][START_REF] Sexton | Thermal ignition in rectangular and triangular regions[END_REF] have assessed thermal ignition models with regard to multiplicity and stability of these steady states using the dimensionless ambient temperature.

In terms of the following dimensionless quantities:

u = RT E , U = RTa E , λ = Ra 2 0 QC0A KE
the governing steady state of heat balance are reduced to:

∆u + λ exp( -1 u ) = 0 x ∈ Ω ∂u ∂n + Bi(u -U ) = 0 x ∈ ∂Ω [6]
The aim is to determine the critical value of U , say U crit , beyond which multiple solutions of the problem (6) occur.

In this paper we first show how to use the finite element method to find an approximation of the solution of the problem (6). This can be done using a fixed point iteration method that we will describe below. Then, we deal with the study of the stability and bifurcation using the block Lanczos algorithm.

Computational method 2.1. Variational Formulation :

In this section we first derive a variational formulation of the boundary value problem (6) subject to Robin boundary condition.

We multiply the equation (6.1) by a test function v H 1 (Ω) multiply, integrate and use Green's formula to obtain:

Ω ∇u∇v dx+Bi ∂Ω uvdS+ Ω λ exp(- 1 u ) vdx = BiU ∂Ω vdS , ∀v ∈ H 1 (Ω)[7]
It is obvious that the left term of this expression is not a bilinear form. Therefore the Lax-Milgram theorem cannot be applied here.

To solve this problem, we must linearize the problem in order to bring the solution to a series of linear problem converging to the solution of the initial nonlinear problem. We have several approaches to deal with such problems in particular Picard method or Newton method.

Picard's method :

We briefly recall the principle of Picard's method (or successive iterations), which is a variant of the fixed-point method.

Let the matrix system:

Ku = f [8]
as a result of the finite element discretization of a non-linear problem, where K is the stiffness matrix, u the unknown vector and f the 'load vector'. K may be linear or not.

The residue R, which is defined by: R = f -Ku is a measure of the distance from the solution.

The methods of solving nonlinear systems are all iterative. At every step, they calculate du , the solution increment defined by: du = u i+1 -u i , where i is the index of the iteration.

One seeks the convergence of the system by successive approximations. This convergence can be measured by a standard of the increment or a norm of R. du residue, the solution increment can also be seen as a descent direction, down to the zero residue, i.e. to the solution.

The method of Picard is a fixed point method [START_REF] Reddy | The Finite Element Method in Heat Transfer and Fluid Dynamics[END_REF]. It is defined by:

K(u i )u i+1 = f (u i ) given u 0 [9]
which can also be written as:

K(u i )du = R(u i )
for the solution increment du. an accuracy of ε is required the process may be terminated by setting a criterion:

max u i+ -u i ≺ ε
With . is the Euclidian norm.

Rather than being recalculated at each iteration, the matrix K can be maintained constant over a number of sub-iterations.

At convergence, and if the solution exists, then it satisfies [START_REF] Chambré | On the solution of the Poisson-Boltzmann equation with application to the theory of thermal explosions[END_REF]. The advantage of this method lies in the existence of a large enough radius of convergence that allows us not to worry too much about the initial approximation. In return, its convergence rate remains often too slow for practical applications of interest. We can then use a relaxation method speeding a little bit the scheme [START_REF] Scatter | Path following, critical points and shape factors[END_REF]. The latter is defined as follows:

K(u i )u * = f (u i ) , u 0 given u i+1 = αu i+1 + (1 -α)u * , 0 ≤ α ≤ 1 [10]

Method of solution

After the choice of appropriate finite elements by a method adapted from discretization in space, the system of nonlinear equations resulting system ( 7) is solved using techniques such as Newton's method or the method of Picard.

The discrete problem arises as usual by restricting V = H 1 (Ω) to a discrete space V h according to some mesh with some element type. Similarly, we let u d be the discrete solution, so the discrete nonlinear problem is then written as: find u d such that:

F (u d , v) = L(v).
with u d = N j=1 α j Φ j . Since F is a nonlinear function of u d , the variational statement gives rise to a system of nonlinear algebraic equations. Picard's method can be used in alternative way for solving a nonlinear PDE problem.

We simply use a known, previous solution in the nonlinear terms so that these terms become linear in the unknown u d .For our particular problem, we use a known, previous solution in the coefficient e -1 u . More precisely, given a solution u k from iteration k, we seek a new (hopefully improved) solution u k+1 in iteration k + 1 such that u k+1 solves the linear problem,

∆u k+1 + λ exp( -1 u k ) = 0 x ∈ Ω ∂u k+1 ∂n + Bi(u k+1 -U ) = 0 x ∈ ∂Ω [11]
The iterations require an initial guess u 0 . The hope is that u k -→ u d as k → ∞, and that u k+1 is sufficiently close to the exact solution u d of the discrete problem after just a few iterations.

We can easily formulate a variational problem for u k+1 from the last equation. Equivalently, we can approximate e -1 u to obtain the same linear variational problem. Therefore, the problem consists of seeking u k+1 V h such that :

a(u k+1 , v) = L(v) , ∀v ∈ V h [12] whith : a(u k+1 , v) = Ω ∇u k+1 ∇v dx + Bi ∂Ω u k+1 vdS L(v) = BiU ∂Ω vdS -λ Ω v exp( -1 u k ) dx
The iterations can be stopped when ε = u k+1 -u k ≺ tol where tol is a small tolerance, or when the number of iterations exceed some critical limit. The latter case will pick up divergence of the method or unacceptable slow convergence.

Note that the existence and uniqueness of the solution of problem ( 12) is established through the Lax-Milgram theorem [START_REF] Daners | Robin boundary value problems on arbitrary domains, Trans[END_REF].

Thus, Picard's linearization reads as follows:

Find u k+1 such that :

Ω ∇u k+1 ∇v dx + Bi ∂Ω u k+1 vdS = BiU ∂Ω vdS -λ Ω v exp( -1 u k ) dx , given u 0

Numerical results

In this section we first present the spatial discretization using finite element method with P1 continuous piecewise linear functions. Then we develop all the steps of the FreeFem++ code to solve the problem by using the technique of mesh adaptation and at the end we present some numerical results.

To discretize (12), let T h denote a regular, quasi uniform triangulation of Ω with triangles of maximum size h ≺ 1. Let V h = v h ∈ C 0 (_Ω), v h/T ∈ P 1 (T ), ∀T ∈ T h denote a finite dimensional subspace of H 1 (Ω) where P 1 is the set of polynomials of degrees ≤ 1 .

Let Ω be the rectangle [0, 1] × [0, 1] , the triangulation T h of Ω is automatically generated by using a FreeFem command (Figure1) . In Figures 234, we present the evolution of the temperature profiles according to different parameters (λ, U , Bi). These solutions are in agreement with those obtaind by Balakrishnan et al. [START_REF] Balakrishnan | Critical values for some nonclass A geometries in thermal ignition theory[END_REF] by means of the shape factor method based on finite difference approximation.

Conclusion

Analysis has been carried out for steady state of exothermic chemical reaction of combustible material, in an arbitrary domain, with convective boundary conditions. The nonlinear PDE governing the problem are solved numerically using a FEM approach based on Picard's method. We have presented a numerical approach with FreeFem software, and we established the adequacy of the chosen finite element discretization by comparing the results with those of Balakrishnan et al. [START_REF] Balakrishnan | Critical values for some nonclass A geometries in thermal ignition theory[END_REF]. Thermal stability of such problems remains to be investigated. This can be undertaken using methods such as block Lanczos algorithm or Hermite-Padé approximations.

Figure 1 :

 1 Figure 1 : mesh Th.

Figure 2 :

 2 Figure 2 : Temperature distribution profile for λ = 1, U = 1, Bi = 1.

Figure 3 :

 3 Figure 3 : Temperature distribution profile for λ =10, U = 10, Bi = 1000.

Figure 4 :

 4 Figure 4 : Temperature distribution profile for λ= 9, U = 1, Bi = 1.