We consider the singularly perturbed boundary value problem:

where $\epsilon$$\epsilon$epsi\varepsilon$\epsilon$ is a small positive parameter, $a(x)>0$$a(x)>0$a(x) > 0a(x)>0$a(x)>0$ for all $x\in [0,1]$$x\in [0,1]$x in[0,1]x \in[0,1]$x\in [0,1]$, and the functions $a$$a$aa$a$ and $f$$f$ff$f$ are sufficiently smooth. The solution of (P1) has a boundary layer at $x=1$$x=1$x=1x=1$x=1$.

It has long been recognized that difficulties can arise when certain "centered" finite-difference and finite-element methods are applied to (P1) when the diffusion coefficient $\epsilon$$\epsilon$epsi\varepsilon$\epsilon$ is small. In particular, such schemes when applied to (P1) on a uniform grid have an inherent formal cell Reynolds number limitation. Namely, with a uniform mesh length $h$$h$hh$h$ and $a(x)$$a(x)$a(x)a(x)$a(x)$ constant, one finds that the cell Reynolds number $\frac{ah}{\epsilon }$$\frac{ah}{\epsilon }$(ah)/(epsi)\frac{a h}{\varepsilon}$\frac{ah}{\epsilon }$ must be bounded by some constant depending on the scheme in order to avoid spurious oscillations or gross inaccuracies. For small $\epsilon$$\epsilon$epsi\varepsilon$\epsilon$ this requires a prohibitive number of grid points and so alternative approaches have been developed. One approach is to use a nonuniform mesh (which must be appropriately chosen) which is very fine "in the boundary layer" and coarser elsewhere. Another approach has been to devise schemes which have no formal cell Reynolds number limitation. Schemes of this type have been constructed by using uncentered ("upwind") differencing for the first derivative term, or, more generally, by adding an "artificial viscosity" to the diffusion coefficient $\epsilon$$\epsilon$epsi\varepsilon$\epsilon$, e.g., [10].

When the discretization is made by finite differences, Temam introduced in [11] the concept of Incremental Unknowns (IU in short). The idea, which stems from dynamical systems approach, consists in writing the approximate solution $u_i$$u_i$u_(i)u_{i}$u_i$ in the form $u_i=y_i+z_i$$u_i=y_i+z_i$u_(i)=y_(i)+z_(i)u_{i}=y_{i}+z_{i}$u_i=y_i+z_i$, where $z$$z$zz$z$ is a small increment. Passing from the nodal unknowns $u_i$$u_i$u_(i)u_{i}$u_i$ to the IUs $(y_i,z_i)$$\left(y_i,z_i\right)$(y_(i),z_(i))\left(y_{i}, z_{i}\right)$(y_i,z_i)$ amounts to a linear change of variables, that is to say, in the language of linear algebra, to the construction of a preconditioner. Many numerical simulations have shown the efficiency of such induced preconditioners.

Numerical solution of a problem such as (P1) using Incremental Unknowns has been considered in [6] and [7] but the IU's that have been used in these articles were connected to the Laplacian only: they were induced only by the discretization matrix of the Laplacian.

In [3], the authors propose a construction of IUs that are more adapted to the problem in the sense that they take into account the convection term in the construction of the IUs, thus leading to the use of an adapted interpolator and of a hierarchichal preconditioner.

We propose in this paper a different approach for the construction of adapted incremental unknowns for (P1): we first make a change of variable (assuming that there exists a grid function, see below) and then discretize the problem. This change of variable allows us to work on a uniform grid, which makes the calculations (in particular via Taylor's expansions) easier; the effects of the grid will then be on the coefficients of the differential operators.

Let $g:[0,1]\to [0,1]$$g:[0,1]\to [0,1]$g:[0,1]rarr[0,1]g:[0,1] \rightarrow[0,1]$g:[0,1]\to [0,1]$ such that $g(0)=0,g(1)=1,\mathrm{\exists }{g}^{\prime }(y)$$g(0)=0,g(1)=1,\mathrm{\exists }{g}^{\prime }(y)$g(0)=0,g(1)=1,EEg^(')(y)g(0)=0, g(1)=1, \exists g^{\prime}(y)$g(0)=0,g(1)=1,\mathrm{\exists }{g}^{\prime }(y)$ for all $x\in [0,1]$$x\in [0,1]$x in[0,1]x \in[0,1]$x\in [0,1]$ and $0<J(y):=g^{\prime }(y)<M,\mathrm{\forall }y\in (0,1)$$0<J(y):=g^{\prime }(y)<M,\mathrm{\forall }y\in (0,1)$0 < J(y):=g^(')(y) < M,AA y in(0,1)0<J(y):=g^{\prime}(y)<M, \forall y \in(0,1)$0<J(y):=g^{\prime }(y)<M,\mathrm{\forall }y\in (0,1)$. Furthermore, in order to obtain a finer resolution near the boundary layer, we assume that $J(1)=0$$J(1)=0$J(1)=0J(1)=0$J(1)=0$.

With the change of variable $x=g(y)$$x=g(y)$x=g(y)x=g(y)$x=g(y)$ we obtain from (P1), the following problem:

Let $\varphi (y)=\frac{1}{J(y)}{v}^{\prime }(y)$$\varphi (y)=\frac{1}{J(y)}{v}^{\prime }(y)$phi(y)=(1)/(J(y))v^(')(y)\phi(y)=\frac{1}{J(y)} v^{\prime}(y)$\varphi (y)=\frac{1}{J(y)}{v}^{\prime }(y)$. By integration of this equation we obtain:

which can also be written,

Let us consider the Taylor expansion of the function $\varphi$$\varphi$phi\phi$\varphi$ :

$\varphi (s)=\varphi (y)+(s-y){\int }_y^{y+h}J(s)\mathrm{d}s{\varphi }^{\prime }(y)+\frac{1}{2}(s-y{)}^2{\varphi }^{\prime \prime }(y+{\theta }^{-}(s-y))$$\varphi (s)=\varphi (y)+(s-y){\int }_y^{y+h} J(s)\mathrm{d}s{\varphi }^{\prime }(y)+\frac{1}{2}(s-y{)}^2{\varphi }^{\prime \prime }\left(y+{\theta }^{-}(s-y)\right)$phi(s)=phi(y)+(s-y)int_(y)^(y+h)J(s)dsphi^(')(y)+(1)/(2)(s-y)^(2)phi^('')(y+theta^(-)(s-y))\phi(s)=\phi(y)+(s-y) \int_{y}^{y+h} J(s) \mathrm{d} s \phi^{\prime}(y)+\frac{1}{2}(s-y)^{2} \phi^{\prime \prime}\left(y+\theta^{-}(s-y)\right)$\varphi (s)=\varphi (y)+(s-y){\int }_y^{y+h}J(s)\mathrm{d}s{\varphi }^{\prime }(y)+\frac{1}{2}(s-y{)}^2{\varphi }^{\prime \prime }(y+{\theta }^{-}(s-y))$,

Injecting these expressions of $\varphi$$\varphi$phi\phi$\varphi$ in (1), we find from (2) and (3):

If we denote

we obtain

and

which is a system of two equations with two unknowns $\varphi (y)$$\varphi (y)$phi(y)\phi(y)$\varphi (y)$ and ${\varphi }^{\prime }(y)$${\varphi }^{\prime }(y)$phi^(')(y)\phi^{\prime}(y)${\varphi }^{\prime }(y)$ whose determinant will be denoted by $r(y)$$r(y)$r(y)r(y)$r(y)$.

Proposition 1. The functions $a^{\pm },b^{+}$$a^{\pm },b^{+}$a^(+-),b^(+)a^{ \pm}, b^{+}$a^{\pm },b^{+}$are non-negative functions and the function $r$$r$rr$r$ is a negative one.

Since $r<0$$r<0$r < 0r<0$r<0$ the system has the solution,

and

We consider the following approximation for ${\varphi }^{\prime }(y)$${\varphi }^{\prime }(y)$phi^(')(y)\phi^{\prime}(y)${\varphi }^{\prime }(y)$ and $\varphi (y)$$\varphi (y)$phi(y)\phi(y)$\varphi (y)$ :

and

which is a backward type approximation.
If we substitute in (P2), we obtain the approximation:

Let $h=\frac{1}{2N}$$h=\frac{1}{2N}$h=(1)/(2N)h=\frac{1}{2 N}$h=\frac{1}{2N}$, and for $j=0,1,2,\dots ,2N,y_j=jh$$j=0,1,2,\dots ,2N,y_j=jh$j=0,1,2,dots,2N,y_(j)=jhj=0,1,2, \ldots, 2 N, y_{j}=j h$j=0,1,2,\dots ,2N,y_j=jh$. For $y=y_j$$y=y_j$y=y_(j)y=y_{j}$y=y_j$ in (1) we have for $j=1,\dots ,2N-1$$j=1,\dots ,2N-1$j=1,dots,2N-1j=1, \ldots, 2 N-1$j=1,\dots ,2N-1$ :

where in general $f_j=f\left(y_j\right)$$f_j=f\left(y_j\right)$f_(j)=f(y_(j))f_{j}=f\left(y_{j}\right)$f_j=f\left(y_j\right)$.
Since $v\left(y_j\right):=u\left(g\left(y_j\right)\right)=u\left(x_j\right)$$v\left(y_j\right):=u\left(g\left(y_j\right)\right)=u\left(x_j\right)$v(y_(j)):=u(g(y_(j)))=u(x_(j))v\left(y_{j}\right):=u\left(g\left(y_{j}\right)\right)=u\left(x_{j}\right)$v\left(y_j\right):=u\left(g\left(y_j\right)\right)=u\left(x_j\right)$, for $j=0,1,\dots ,2N$$j=0,1,\dots ,2N$j=0,1,dots,2Nj=0,1, \ldots, 2 N$j=0,1,\dots ,2N$, we can write (7):

Remark 1. If $g(y):=y$$g(y):=y$g(y):=yg(y):=y$g(y):=y$ then we obtain the classical upwind scheme.
Let ${\alpha }_j=\frac{\epsilon a_j^{-}(x_{j+1}-x_{j-1})}{2r_j{J}_j}$${\alpha }_j=\frac{\epsilon a_j^{-}\left(x_{j+1}-x_{j-1}\right)}{2r_j{J}_j}$alpha_(j)=(epsia_(j)^(-)(x_(j+1)-x_(j-1)))/(2r_(j)J_(j))\alpha_{j}=\frac{\varepsilon a_{j}^{-}\left(x_{j+1}-x_{j-1}\right)}{2 r_{j} J_{j}}${\alpha }_j=\frac{\epsilon a_j^{-}(x_{j+1}-x_{j-1})}{2r_j{J}_j}$ and ${\beta }_j=(\frac{a_j}{a_j^{-}}-\frac{\epsilon a_j^{+}}{r_j{J}_j})\frac{(x_{j+1}-x_{j-1})}{2}$${\beta }_j=\left(\frac{a_j}{a_j^{-}}-\frac{\epsilon a_j^{+}}{r_j{J}_j}\right)\frac{\left(x_{j+1}-x_{j-1}\right)}{2}$beta_(j)=((a_(j))/(a_(j)^(-))-(epsia_(j)^(+))/(r_(j)J_(j)))((x_(j+1)-x_(j-1)))/(2)\beta_{j}=\left(\frac{a_{j}}{a_{j}^{-}}-\frac{\varepsilon a_{j}^{+}}{r_{j} J_{j}}\right) \frac{\left(x_{j+1}-x_{j-1}\right)}{2}${\beta }_j=(\frac{a_j}{a_j^{-}}-\frac{\epsilon a_j^{+}}{r_j{J}_j})\frac{(x_{j+1}-x_{j-1})}{2}$ for $j=1,\dots ,2N-1$$j=1,\dots ,2N-1$j=1,dots,2N-1j=1, \ldots, 2 N-1$j=1,\dots ,2N-1$. We obtain a finite dimensional linear system which can be written as

By using trapezoidal rule:

we have
so,

As usual when an IU method is implemented, two different kinds of unknowns must be distinguished: those associated with the coarse grid components which are on $G_c$$G_c$G_(c)G_{c}$G_c$, and whose indices are even and those associated with the complementary points (odd indices) which are on $G_f\mathrm{\setminus }{G}_c$$G_f\mathrm{\setminus }{G}_c$G_(f)\\G_(c)G_{f} \backslash G_{c}$G_f\mathrm{\setminus }{G}_c$;