2004-Fez conference on Differential Equations and Mechanics Electronic Journal of Differential Equations, Conference 11, 2004, pp. 53–59. ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu ftp ejde.math.txstate.edu (login: ftp)

ERROR ESTIMATE FOR THE CHARACTERISTIC METHOD

INVOLVING OLDROYD DERIVATIVE IN A TENSORIAL TRANSPORT PROBLEM

MOHAMMED BENSAADA, DRISS ESSELAOUI, PIERRE SARAMITO

Abstract. An optimal a priori error estimate O(h^{k+1} + Δt), result is presented for a tensor problem involving Oldroyd derivative when using a suitable characteristic method and a finite element method. To conclude, we present results of numerical tests which confirm the previous estimates. Our long time goal is to deal with the viscoelastic fluid flow problem.

1. Introduction

The Oldroyd derivative of a symmetric tensor σ is defined by

\frac{𝒟_{a}σ}{𝒟t} = \frac{∂σ}{∂t} + (u.∇)σ + σM_{a}(u) + M_{a}^{T} (u)σ,  (1.1)

where u is a given velocity field, M_{a}(u) = ((1 − a)∇u − (1 + a)∇u^{T} ) /2 and a ∈ [−1, 1] is the parameter of the Oldroyd derivative. In this paper, we suppose also, for simplicity, that u = 0 on ∂Ω. Problems involving the Oldroyd derivative appear in viscoelasticity (non-Newtonian polymer melt flow problems, see e.g. [6]), in turbulence modelling (R_{ij} −ϵ models) or in liquid crystals modelling. Let Ω be a bounded polygonal subset of ℝ^{d}, d = 1, 2 or 3, T > 0 be a time constant, γ a given symmetric tensor defined in Ω×]0, T [, σ_{0} given in Ω. and λ > 0. The aim of this paper is to study, as a preliminary result, the approximation of the linear transport equation involving the Oldroyd derivative: find σ defined in Ω×]0, T [ such that

\begin{matrix}λ\frac{𝒟_{a}σ}{𝒟t} + σ = γ in Ω×]0, T [, \\ σ(0) = σ_{0} in Ω. \end{matrix}  (1.2)

2000 Mathematics Subject Classification. 65M60, 65M25, 65M15, 76A10.

Key words and phrases. Finite element method; method of characteristics; error bound;

viscoelastic fluids.

\begin{matrix}c \\ ◯ \end{matrix}2004 Texas State University - San Marcos.

Published October 15, 2004.

Supported by grant STCI01/03 from the CNRST-CNRS.

and its discrete counterpart: find σ_{h}^{(n)} ∈ T_{h}, 1 ≤ n ≤ N, such that, for all τ_{h} ∈ T_{h},

(λ\frac{σ_{h}^{(n+1)} − R_{Δt}^{(n)} × σ_{h}^{(n)} ∘ X_{Δt}^{(n)} × (R_{Δt}^{(n)})^{T}}{Δt} + σ_{h}^{(n+1)}, τ_{h}) = (γ(t_{n+1}), τ_{h}) , (1.3)

and σ_{h}^{(0)} = π_{h}σ_{0}, where N ≥ 1, Δt = T/N, t_{n} = nΔt, X_{Δt}^{(n)}(x) = x − Δt u(x, t_{n+1}), R_{Δt}^{(n)}(x) = I −Δt M_{a}^{T} (u)(x, t_{n+1}), and T_{h} denotes the space of continuous piecewise polynomial symmetric tensors T_{h} = \{τ ∈ (C^{0}(\overline{Ω}) ∩ L^{2}(Ω))^{d×d}; τ_{|K} ∈ P_{k}, ∀K ∈ 𝒯_{h} }.

Here, (𝒯_{h} )_{h>0} denotes a suitable family of regular triangulations of \overline{Ω}. The Lagrange interpolation operator from L^{2} tensors into T_{h} is denoted by π_{h}, k ≥ 1.

The characteristic method [1] has been proposed for the numerical treatment of convected-dominated flows and transport equations. It is based on an approximation of the material derivative \frac{D}{Dt} = \frac{∂}{∂t} + (u.∇):

\frac{Dφ}{Dt}(x, t) ≈ \frac{φ(x, t) − φ(X(x, t; s), s)}{t − s}  (1.4)

The trajectory X(x, s; .) is defined for all (x, s) ∈ \overline{Ω} × [0, T] by

\begin{matrix}\frac{∂X}{∂t} (x, s; t) = u(X(x, s; t), t), t ∈]0, T [, \\ X(x, s; s) = x. \end{matrix}  (1.5)

In 1987, D. Esselaoui and M. Fortin [3] extended the characteristic method for the approximation of the Oldroyd derivative of a symmetric tensor. For all fixed s ∈ [0, T], these authors considered the following transformed problem: find σ̂(.,.;s) defined in Ω×]0, T [ such that

\begin{matrix}λ\frac{Dσ̂}{Dt} (x, t; s) + σ̂(x, t; s) = R(x, t; s)γ(x, t)R^{T} (x, t; s), (x, t) ∈ Ω×]0, T [, \\ σ̂(x, 0; s) = R(x, 0; s)σ_{0}(x)R^{T} (x, 0; s), x ∈ Ω, \end{matrix} (1.6)

where σ̂(x, t; s) = R(x, t; s)σ(x, t)R(x, t; s)^{T} and the tensor flow R(x, .; s) is defined for all (x, s) ∈ \overline{Ω} × [0, T] by

\begin{matrix}\frac{DR}{Dt}(x, t; s) = R(x, t; s)M_{a}^{T} (u)(x, t), t ∈]0, T [, \\ R(x, s; s) = I. \end{matrix}  (1.7)

The Oldroyd derivative has been replaced by a material derivative of a tensor, suitable for the characteristic method, as (1.3). A short computation shows that

\frac{D^{m}σ̂}{Dt^{m}} (x, t; s) = R^{m}(x, t; s)\frac{𝒟_{a}^{m}σ}{𝒟t^{m}} (x, t)(R^{m})^{T} (x, t; s), ∀m ∈ ℕ,  (1.8)

and thus problems (1.2) and (1.6) are equivalent. Numerical computations involving the scheme (1.6) has been already performed in [3] while the corresponding numerical analysis was not yet available. We show a 𝒪 (h^{k+1} + Δt) optimal estimate. Moreover, since the trajectories and the transformation are both approximated, the scheme (1.6) is of practical interest. To conclude, we present results of numerical tests that confirm the previous estimates.

2. Error estimate

Let ∥ · ∥ the L^{2}(Ω) norm, ∥ · ∥_{∞} the L^{∞}(Ω) one and | · |_{m,p,Ω} the W^{m,p}(Ω) semi-norm, for m ≥ 0 et p ∈ [1, ∞]. Also, C^{0,1}(\overline{Ω}) denotes the space of lipschitzian functions \overline{Ω}. For a Banach space Y, let us denote C(Y) the space C([0, T], Y). We suppose also that the data u ∈ C(C^{0,1}(\overline{Ω})): the existence and the continuity of x → X(x, s; t) follow then from the Cauchy-Lipschitz theorem. Let us denote finally X^{(n)}(x) = X(x, t_{n+1}; t_{n}) and R^{(n)}(x) = R(X^{(n)}(x), t_{n}; t_{n+1}). In this paper, C_{i}, i ∈ ℕ is a positive constant, independent of h and Δt.

For a Banach space Y, with norm ∥.∥_{Y} and 1 ≤ p < ∞ we introduce

l^{p}(0, T; Y) = {φ : (t_{1}, …, t_{N}) → Y; ∥φ∥_{l^{p}(0,T;Y)} = ( ∑_{n=1}^{N} ∥φ(t_{i})∥_{Y}^{p} Δt)^{1/p} < ∞},

l^{∞}(0, T; Y) = {φ : (t_{1}, …, t_{N}) → Y; ∥φ∥_{l^{∞}(0,T;Y)} = max_{1≤i≤N} ∥φ(t_{i})∥_{Y} < ∞}.

Theorem 2.1 (Error estimate). Let σ and σ_{h} the solutions of (1.2) and (1.3), respectively. Suppose that u is in (C(C^{0,1}(\overline{Ω})) ∩ W^{1,∞}(W^{1,∞}(Ω))^{d} and σ is in (W^{1,∞}(H^{m+1}(Ω)) ∩ W^{2,∞}(L^{2}(Ω)))^{d×d}, with m ≥ 0. Then, there exist three positive constants Δt_{0}, h_{0} and c, independent of h and Δt, such that, if Δt ≤ Δt_{0} and h ≤ h_{0}, we have

∥σ − σ_{h}∥_{l^{∞}(0,T;L^{2}(Ω))} ≤ c (h^{r+1} + Δt) ,  (2.1)
∥σ − σ_{h}∥_{l^{2}(0,T;L^{2}(Ω))} ≤ c (h^{r+1} + Δt) ,  (2.2)

where r = min(k, m).

Proof. Let us introduce σ̃_{h}(t) = π_{h}σ(t), t ∈ [0, T], σ̃_{h}^{(n)} = σ̃_{h}(t_{n}) and ε_{h}^{(n)} = σ_{h}^{(n)} − σ̃_{h}^{(n)}. By a development and using (1.3), we get

(λ\frac{ε_{h}^{(n+1)} − R_{Δt}^{(n)} × ε_{h}^{(n)} ∘ X_{Δt}^{(n)} × (R_{Δt}^{(n)})^{T}}{Δt} + ε_{h}^{(n+1)}, ε_{h}^{(n+1)}) = (ρ_{h}^{(n+1)}, ε_{h}^{(n+1)}),

(2.3)

where

ρ_{h}^{(n+1)} = γ(t_{n+1}) − σ̃_{h}^{(n+1)} − λ\frac{σ̃_{h}^{(n+1)} − R_{Δt}^{(n)} × σ̃_{h}^{(n)} ∘ X_{Δt}^{(n)} × (R_{Δt}^{(n)})^{T}}{Δt} . (2.4)

From the Cauchy-Schwartz inequality and using the identity ab ≤ (a^{2} + b^{2})/2,

(R_{Δt}^{(n)}ε_{h}^{(n)} ∘ X_{Δt}^{(n)}(R_{Δt}^{(n)})^{T} , ε_{h}^{(n+1)})

≤ \frac{1}{2}(∥R_{Δt}^{(n)}ε_{h}^{(n)} ∘ X_{Δt}^{(n)} (R_{Δt}^{(n)})^{T} ∥^{2} + ∥ε_{h}^{(n+1)}∥^{2})

and

(ρ_{h}^{(n+1)}, ε_{h}^{(n+1)}) ≤ \frac{1}{2} (\frac{1}{2}∥ρ_{h}^{(n+1)}∥^{2} + 2∥ε_{h}^{(n+1)}∥^{2}) .

Then, (2.3) becomes

∥ε_{h}^{(n+1)}∥^{2} ≤ ∥R_{Δt}^{(n)} × ε_{h}^{(n)} ∘ X_{Δt}^{(n)} × (R_{Δt}^{(n)})^{T} ∥^{2} + \frac{Δt}{2λ} ∥ρ_{h}^{(n+1)}∥^{2} .

From lemma 3.1 (paragraph 2), for Δt small enough, we get

∥ε_{h}^{(n+1)}∥^{2} ≤ (1 + C_{0}Δt) ∥ε_{h}^{(n)}∥^{2} + \frac{Δt}{2λ} ∥ρ_{h}^{(n+1)}∥^{2} .

From the discrete Gronwall lemma, and using n ≤ N = T/Δt,

∥ε_{h}^{(n)}∥^{2} ≤ \frac{e^{C_{0}T}}{2λ} ∥ρ_{h}∥_{l^{2}(L^{2})}^{2}  (2.5)

where ρ_{h} = (ρ_{h}^{(1)}, …, ρ_{h}^{(N)}). It still remains to bound the right-hand side. From (1.6) with s = t = t_{n+1}, we get

λ\frac{Dσ̂}{Dt} (x, t_{n+1}; t_{n+1}) + σ(x, t_{n+1}) − γ(x, t_{n+1}) = 0.  (2.6)

Adding (2.6) to the expression (2.4) of ρ_{h}^{(n+1)},

ρ_{h}^{(n+1)} = (σ − σ̃_{h})(t_{n+1})

+ λ\{\frac{Dσ̂}{Dt} (x, t_{n+1}; t_{n+1}) − \frac{σ̃_{h}^{(n+1)} − R_{Δt}^{(n)} × σ̃_{h}^{(n)} ∘ X_{Δt}^{(n)} × (R_{Δt}^{(n)})^{T}}{Δt} \}.

Let us introduce the splitting ρ_{h}^{(n+1)} = ζ_{h}^{(n+1)} + η_{h}^{(n+1)}, where

ζ_{h}^{(n+1)}

= (σ − σ̃_{h})(t_{n+1}) + \frac{λ}{Δt} \{(σ − σ̃_{h})(t_{n+1}) − R^{(n)} × (σ − σ̃_{h})(X^{(n)}, t_{n}) × (R^{(n)})^{T} \} ,

and

η_{h}^{(n+1)}

= λ\{\frac{Dσ̂}{Dt} (x, t_{n+1}; t_{n+1}) − \frac{σ(t_{n+1}) − R^{(n)} × σ (X^{(n)}, t_{n}) × (R^{(n)})^{T}}{Δt} \}

− \frac{λ}{Δt}\{R^{(n)} × σ̃_{h} (X^{(n)}, t_{n}) × (R^{(n)})^{T} − R_{Δt}^{(n)} × σ̃_{h} (X_{Δt}^{(n)}, t_{n}) × (R_{Δt}^{(n)})^{T} \}.

On the one hand, using a classical interpolation result [2], we have ∥φ − π_{h}φ∥ ≤ C_{4}h^{r+1}|φ|_{r+1,2,Ω}, for all φ ∈ H^{r+1}(Ω). Then, lemma 3.2 (paragraph 2) yields ∥ζ_{h}^{(n+1)}∥ = 𝒪(h^{r+1}). On the other hand, the lemma 3.2 and 3.3 (paragraph 2) and the continuity of π_{h} in H^{1}(Ω) give ∥η_{h}^{(n+1)}∥ = 𝒪(Δt). Thus, the result yields by reporting ∥ρ_{h}^{(n+1)}∥ = 𝒪(h^{r+1} + Δt) in (2.5). □

3. Auxiliary lemma

Lemma 3.1 (Estimate on the approximate transformation). There exist two positive constants C_{0} and Δt_{0} such that, if Δt ≤ Δt_{0}, then

∥R_{Δt}^{(n)} × τ ∘ X_{Δt}^{(n)} × (R_{Δt}^{(n)})^{T} ∥^{2} ≤ (1 + C_{0}Δt)∥τ∥^{2}, ∀τ ∈ L^{2}(Ω).

Proof. Let d(x) = dist(x, ∂Ω) for x ∈ Ω. Since u = 0 on ∂Ω, we have

|X_{Δt}^{(n)}(x) − x| = |u(x, t_{n+1})| Δt ≤ ∥u∥_{L^{∞}(W^{1,∞})} d(x)Δt.

This leads to X_{Δt}^{(n)}(x) ∈ Ω for Δt < 1/|u|_{L^{∞}(W^{1,∞})}. Since X_{Δt}^{(n)}(∂Ω) = ∂Ω, and

from the continuity of u(., t_{n+1}), we get X_{Δt}^{(n)}(Ω) = Ω. Let J_{Δt}^{(n)} = det (R_{Δt}^{(n)}). Let

|.| denotes the matrix norm in ℝ^{d×d}. We get

∫_{Ω} |R_{Δt}^{(n)} × τ ∘ X_{Δt}^{(n)} × (R_{Δt}^{(n)})^{T} |^{2} dx
≤ ∥R_{Δt}^{(n)}∥_{∞}^{4} ∫_{Ω} |τ ∘ X_{Δt}^{(n)}|^{2} dx
≤ ∥R_{Δt}^{(n)}∥_{∞}^{4} ∫_{X_{Δt}^{(n)}(Ω)} |τ(y)|^{2} J_{Δt}^{(n)}(y) dy ≤ ∥R_{Δt}^{(n)}∥_{∞}^{4}∥J_{Δt}^{(n)}∥_{∞}^{2}∥τ∥^{2}.

From the definition of R_{Δt}^{(n)}, we get ∥R_{Δt}^{(n)}∥_{∞} ≤ 1 + Δt∥u∥_{L^{∞}(W^{1,∞})} and, for Δt small enough, there exists a positive constant c depending only of u such that ∥J_{Δt}^{(n)}∥_{∞} ≤ 1 + cΔt. Then we obtain the result. □

Lemma 3.2 (Time approximation). There exists K_{1}, C_{1} and C_{2} such that, if Δt ≤ K_{1}/∥u∥_{L^{∞}(W^{1,∞})}, then

(i) for all τ ∈ W^{1,∞}(L^{2}(Ω))^{d×d}

∥τ(., t_{n+1}) − R^{(n)} × τ(X^{(n)}, t_{n}) × (R^{(n)})^{T} ∥ ≤ C_{1}Δt∥\frac{𝒟_{a}τ}{𝒟t} ∥_{L^{∞}(L^{2})}

(ii) for all τ ∈ W^{2,∞}(L^{2}(Ω))^{d×d} we have

∥\frac{Dτ̂}{Dt}(., t_{n+1}; t_{n+1}) − \frac{τ(., t_{n+1}) − R^{(n)} × τ(X^{(n)}, t_{n}) × (R^{(n)})^{T}}{Δt} ∥

≤ C_{2}Δt ∥\frac{𝒟_{a}^{2}τ}{𝒟t^{2}} ∥_{L^{∞}(L^{2})}

Proof. Let f(t) = τ̂(X(x, t_{n+1}, t), t; t_{n+1}). From one hand, remark that f(t_{n+1}) = τ(x, t_{n+1}) and f(t_{n}) = R^{(n)} × τ(X^{(n)}, t_{n}) × (R^{(n)})^{T}. From other hand, by the property of the material derivative f^{(m)}(t) = \frac{D^{m}τ̂}{Dt^{m}}(X(x, t_{n+1}, t), t; t_{n+1}), m ≥ 0. Then, the result yields from the error estimate for the Taylor interpolation polynoms, from (1.8) and from the properties of the solution R(x, .; t_{n+1}) of the linear differential equation (1.7). □

Lemma 3.3 (Trajectories and transformation approximations). If u is in the space W^{1,∞}(W^{1,∞})^{d} and τ in H^{1}(Ω)^{d×d} then there exists C_{3} such that

∥R^{(n)} × τ ∘ X^{(n)} × (R^{(n)})^{T} − R_{Δt}^{(n)} × τ ∘ X_{Δt}^{(n)} × (R_{Δt}^{(n)})^{T} ∥
≤ C_{3}Δt^{2}∥u∥_{W^{1,∞}(W^{1,∞})}∥τ∥_{1,2,Ω}

Proof. Let us consider the following splitting

R^{(n)} × τ ∘ X^{(n)} × (R^{(n)})^{T} − R_{Δt}^{(n)} × τ ∘ X_{Δt}^{(n)} × (R_{Δt}^{(n)})^{T}

= (R^{(n)} − R_{Δt}^{(n)}) × τ ∘ X^{(n)} × (R^{(n)})^{T}

+ R_{Δt}^{(n)} × τ ∘ X^{(n)} × (R^{(n)} − R_{Δt}^{(n)})^{T} + R_{Δt}^{(n)} × (τ ∘ X^{(n)} − τ ∘ X_{Δt}^{(n)}) × (R_{Δt}^{(n)})^{T} .

From the Taylor expansion of f_{1}(t) = R(X(x, t_{n+1}, t), t, t_{n+1}), we get R^{(n)} −R_{Δt}^{(n)} = 𝒪(Δt^{2}). Then, from the expansion of f_{2}(t) = X(x, t_{n+1}, t), we get X^{(n)} − X_{Δt}^{(n)} = 𝒪(Δt^{2}). Finally, from the development of

g(θ) = τ (θX^{(n)}(x) + (1 − θ)X_{Δt}^{(n)}(x)) ,

we have τ ∘ X^{(n)} − τ ∘ X_{Δt}^{(n)} = 𝒪(Δt^{2}). □

4. Numerical experiments

In this paragraph, we present numerical results for d = 2 and k = 1. We choose u = (−x_{2}, x_{1}) Ω =] − 1/2, 1/2[^{2}, T = 2π, γ = 0 and the initial condition σ_{0} is chosen such that the solution of (1.2) is given by

σ(x, t) = \frac{1}{2} exp \{−\frac{t}{λ} − \frac{(x_{1} − x_{1,c}(t))^{2} + (x_{2} − x_{2,c}(t))^{2}}{r_{0}^{2}} \}

× (\begin{matrix}1 + cos(2t) & sin(2t) \\ sin(2t) & 1 − cos(2t) \end{matrix})

where x_{1,c}(t) = x̄_{1,c} cos(t)−x̄_{2,c} sin(t) et x_{2,c}(t) = x̄_{1,c} sin(t)+x̄_{2,c} cos(t) with r_{0} > 0 and (x̄_{1,c}, x̄_{2,c}) ∈ ℝ^{2}. The numerical tests are performed for λ = 1, r_{0} = 1/10 and (x̄_{1,c}, x̄_{2,c}) = (1/4, 0).

The scheme (1.3) is based on the free software rheolef [7]. The computation of the scalar product (R_{Δ}^{(n)}σ_{h}^{(n)} ∘ X_{Δ}^{(n)}(R_{Δ}^{(n)})^{T} , τ_{h}) in the approximate problem (1.3) is not an obvious task (see also [4, 5]). In each triangle, we use the six point fourth order Gauss quadrature formulae.

Since our domain is a square, we split each edge in M segments, and obtain M^{2} small squares of length 1/M. Then we split each of them in two triangles: the step of the mesh is then h = \sqrt{2}/M. Fig. 1.a represents ∥σ_{h} − π_{h}σ∥_{l^{2}(L^{2})} as a function of Δt for four triangulations given by h = \sqrt{2}/2^{i}, 3 ≤ i ≤ 6. We observe that when Δt → 0, the error starts to decrease and then tends to a constant that behaves as h^{2}. A simultaneous choice of (h, Δt) such that Δt = 𝒪(h^{2}) is convenient : let us choose (h, Δt) = (\sqrt{2}/2^{i}, 2π/4^{i}), 3 ≤ i ≤ 6. Fig. 1.b and 1.c presents the corresponding error for the l^{2}(L^{2}) and l^{∞}(L^{2}) norms as a function of h and Δt, respectively. These tests confirm the optimality of our estimates (2.1)-(2.2).

References

[1] M. Bercovier and O. Pironneau. Characteristics and the finite element method. In T. Kawai, editor, Proc. 4th Int. Symp. on finite element methods in flow problems, pages 67–73. Amsterdam: North-Holland, 1982.

[2] S. Brenner and R. Scott. The mathematical theory of finite element methods. Springer-Verlag, 2002.

[3] M. Fortin and D. Esselaoui. A finite element procedure for viscoelastic flows. Int. J. Numer. Methods Fluids, 7:1035–1052, 1987.

[4] K. W. Morton, A. Priestley, and E. Süli. Stability of the Lagrange-Galerkin method with non-exact intergration. Math. Model. Numer. Anal., 22:625–653, 1988.

[5] A. Priestley. Exact projections and the Lagrange-Galerkin method: a realistic alternative to quadrature. J. Comput. Phys., 112(2):316–333, 1994.

[6] P. Saramito. Operator splitting for viscoelastic fluid with a differential constitutive law. C. R. Acad. Sci. Paris, Série II, t. 319, No. 3:267–270, 1994.

[7] P. Saramito, N. Roquet, and E. Jocelyn. rheolef: A finite element environment.

   http://www-lmc.imag.fr/lmc-edp/Pierre.Saramito/rheolef, 2003.

\sqrt{2}/64^{\sqrt{2}/32}^{\sqrt{2}/16}^{h = \sqrt{2}/8}^{∥σ_{h} − π_{h}σ∥_{l^{2}(L^{2}(Ω))}} (a) Δt^{10^{−4}}  ^{10^{−3}}  ^{10^{−2}}  ^{10^{−1}}  ^{1}^{10^{−3}}^{10^{−2}}^{10^{−1}}^{1} ^{^{∥σ_{h} − π_{h}σ∥_{l^{∞}(L^{2}(Ω))}}^{∥σ_{h} − π_{h}σ∥_{l^{2}(L^{2}(Ω))}}^{2} (b) h^{10^{−2}}  ^{10^{−1}}  ^{1}^{10^{−3}}^{10^{−2}}^{10^{−1}}^{1}}

∥σ_{h} − π_{h}σ∥_{l^{∞}(L^{2}(Ω))}^{∥σ_{h} − π_{h}σ∥_{l^{2}(L^{2}(Ω))}}^{1} (c) Δt^{10^{−4}}  ^{10^{−3}}  ^{10^{−2}}  ^{10^{−1}}  ^{1}^{10^{−3}}^{10^{−2}}^{10^{−1}}^{1}

Figure 1. Method convergence versus h and Δt.

Mohammed Bensaada Laboratoire des Sciences de l’Ingénieur, Analyse Numérique et Optimisation (SIANO), Faculté des Sciences, Université Ibn Tofail, B.P.133, 14000-Kénitra, Maroc

E-mail address: m bensaada@hotmail.com

Driss Esselaoui Laboratoire des Sciences de l’Ingénieur, Analyse Numérique et Optimisation (SIANO), Faculté des Sciences, Université Ibn Tofail, B.P.133, 14000-Kénitra, Maroc

E-mail address: desselaoui@yahoo.fr   Fax:  212 37 37 27 70

Pierre Saramito LMC-IMAG, B.P. 53, 38041 Grenoble cedex 9, France

E-mail address: Pierre.Saramito@imag.fr   Fax:  33 4 76 63 12 63