High-order compact schemes and incremental unknowns for 1D convection diffusion problems

Uncategorized

Abstract

Authors

A.C. Muresan (Tiberiu Popoviciu Institute of Numerical Analysis, Romanian Academy)
Adrian Cristian Mureşan

Alina Curseu
Tiberiu Popoviciu Institute of Numerical Analysis, Romanian Academy

Keywords

?

Paper coordinates

A.C. Mureṣan and A. Curṣeu, High-order compact schemes and incremental unknowns for 1D convection diffusion problems, in Proceeding of the ‘Journees Jeunes Numericiens-Colloque en l’honeur de R. Temam’, editor: editions Scientifiques Poitou-Charentes 2001, 77-89.

PDF

About this paper

Journal

Proceeding of the ‘Journees Jeunes Numericiens-Colloque en l’honeur de R. Temam’

Publisher Name
DOI
Print ISSN
Online ISSN

google scholar link

??

Paper (preprint) in HTML form

High-Order Compact Schemes and Incremental Unknowns for 1D Convection Diffusion Problems *

A.C. Mureşan and Alina Curşeu
Abstract

Induced hierarchical preconditioner are investigated for linear 1D convectiondiffusion problems. In the context of high-order compact schemes we give a priori estimates and a uniform estimate of the condition number.

1 Introduction

The utilization of incremental unknowns (IU in short) with multilevel finite differences was proposed by R.Temam in [13] for the integration of elliptic partial differential equations, instead of the usual nodal unknowns. The idea, which stems from dynamical systems approach, consists in writing the approximate solution uiu_{i} in the form ui=yi+ziu_{i}=y_{i}+z_{i}, where zz is a small increment. Passing from the nodal unknowns uiu_{i} to the IUs (yi,zi)\left(y_{i},z_{i}\right) 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.

The incremental unknowns play the role of the small structures. Unfortunately they are not always "small", in certain practical cases. Indeed this situation arises typically when the number of points is not large enough and the discretized function has strong gradients (this is the case of convection diffusion problems). A small step size (or a large number of grid points) can be a limitation for practical applications since the dimension of the system will be increased artificially. The use of high-order schemes can be a solution to the above problem. Exist two main classes of high-order schemes: explicit schemes and compact schemes. Explicit schemes directly compute the numerical derivative by employing large computational stencils for accuracy. High-order schemes achieved in this manner always require non-compact stencils that utilize grid points located beyond those directly adjacent to the node about which the differences are taken. Compact scheme, proposed by Kreiss and Oliger[6] and which was later improved upon by Lele[7], use smaller stencils but requires a scalar tridiagonal or pentadiagonal matrix inversion. Another idea (proposed,

00footnotetext: *2000 MSC Subject Classification: 65L10, 65L12, 65F35.

for convection diffusion problems, by MacKinnon, Carey, Johnson and Langerman [8],[9],[10],[11]) to obtain high-order compact schemes is to operate on the differential equation as an auxiliary relation to obtain expressions for higherorder derivatives in the truncation error. We will use this idea (also presented in [12]) to approximate the solution of the following (1D) convection-diffusion equation:

{u′′(x)+c(x)u(x)=f(x),x(0,1)u(0)=u(1)=0\left\{\begin{array}[]{c}-u^{\prime\prime}(x)+c(x)u^{\prime}(x)=f(x),x\in(0,1)\\ u(0)=u(1)=0\end{array}\right.

where fL2(0,1)f\in L^{2}(0,1) and c(x)L(0,1)c(x)\in L^{\infty}(0,1).
This equation often appears in the description of transport phenomena. The magnitude of c(x)\mathrm{c}(\mathrm{x}) determines the ratio of the convection to diffusion. Of course the 1D problems are not difficult from a computational point of view; we consider them because the algebra is more transparent than for higher dimensions.

Numerical solution of a problem such as (1) using Incremental Unknowns has been considered in [3] and [4] but the IUs that have been used in these articles were connected to the second derivative only.

In [1], the authors propose a construction of IUs that are more adapted to the problem in the sense that they take into account and the convection term in the construction of the IUs.

We consider approximations to (1) on a uniform mesh xi=ih(h=1/2Nx_{i}=ih(h=1/2N, NN a positive integer, i=0,1,,2N)i=0,1,\ldots,2N) having the form

βi(uiui1)+αi(uiui+1)=hfi;i=1,,2N.\beta_{i}\left(u_{i}-u_{i-1}\right)+\alpha_{i}\left(u_{i}-u_{i+1}\right)=hf_{i};i=1,\ldots,2N. (2)

Here uiu(xi),fi=f(xi)u_{i}\simeq u\left(x_{i}\right),f_{i}=f\left(x_{i}\right) and αi\alpha_{i} and βi\beta_{i} are real numbers whose definition depends on the discretization scheme; we have indeed, for example the two following possibilities: the firs one consists in approaching the first derivative by a (centered type) three points scheme and the second one consists in using an uncentered type scheme for the approximation of uu^{\prime} at the grid points. For the moment we do not explicite the discretization scheme.

As usual when an IU method is implemented, two different kinds of unknowns must be distinguished: those associated with the coarse grid components and which are on GcG_{c} (coarse grid), and whose indices are even and those associated with complementary points (odd indices) which are on Gf\Gc(GfG_{f}\backslash G_{c}\left(G_{f}\right. is the fine grid),

.o×o×oo×o..o\times o\times o\ldots o\times o.

Fig.1: Ω=(0,1),x\Omega=(0,1),x : points in Gc,o:G_{c},o: points in Gf\GcG_{f}\backslash G_{c},
If we write the system (2) at the complementary points, we obtain

(α2i+1+β2i+1)u2i+1β2i+1u2iα2i+1u2i+2=hf2i+1\left(\alpha_{2i+1}+\beta_{2i+1}\right)u_{2i+1}-\beta_{2i+1}u_{2i}-\alpha_{2i+1}u_{2i+2}=hf_{2i+1}

Hence, assuming that α2i+1+β2i+10\alpha_{2i+1}+\beta_{2i+1}\neq 0, we have

u2i+1=1α2i+1+β2i+1(β2i+1u2i+α2i+1u2i+2)+1α2i+1+β2i+1hf2i+1u_{2i+1}=\frac{1}{\alpha_{2i+1}+\beta_{2i+1}}\left(\beta_{2i+1}u_{2i}+\alpha_{2i+1}u_{2i+2}\right)+\frac{1}{\alpha_{2i+1}+\beta_{2i+1}}hf_{2i+1}

We note that u2i+1u_{2i+1} is expressed as the sum of a convex combination of u2iu_{2i} and u2i+2u_{2i+2} which is nothing but a bilinear interpolation scheme, and a correction term whose order is connected to the order of the interpolation scheme. If we set

z2i+1=u2i+11α2i+1+β2i+1(β2i+1u2i+α2i+1u2i+2)z_{2i+1}=u_{2i+1}-\frac{1}{\alpha_{2i+1}+\beta_{2i+1}}\left(\beta_{2i+1}u_{2i}+\alpha_{2i+1}u_{2i+2}\right) (3)

then the system (1), at the complementary points, is reduced to

z2i+1=1α2i+1+β2i+1hf2i+1z_{2i+1}=\frac{1}{\alpha_{2i+1}+\beta_{2i+1}}hf_{2i+1}

The numbers z2i+1z_{2i+1} are the incremental unknowns attached to the (discrete) problem (2). They depend closely on the scheme used for the discretization.

We can obviously repeat recursively the process described above and if the coarsest grid is reduced to one point only, the preconditioned matrix becomes diagonal.

From the point of view of the matriceal framework, this construction can be summarized by the determination of two matrices SS and Tt{}^{t}T under and upper triangular respectively such that TtAS{}^{t}TAS is bloc diagonal, A being the discretization matrix.

We first consider two grid levels. The discretization matrix AA is written with the hierarchical ordering (considering first the coarse grid unknowns and then the complementary ones) in the form

A~=(Λ1B1B2Λ2)\widetilde{A}=\left(\begin{array}[]{ll}\Lambda_{1}&B_{1}\\ B_{2}&\Lambda_{2}\end{array}\right)

where Λi,i=1,2\Lambda_{i},i=1,2 are invertible diagonal matrices.

Construction of SS

We want to construct a matrix SS of the form:

S=(I0G1I)S=\left(\begin{array}[]{rr}I&0\\ G_{1}&I\end{array}\right)

and such that A~S\widetilde{A}S is upper triangular. We have

A~S=(Λ1B1B2Λ2)(I0G1I)=(Λ1+B1G1B1B2+Λ2G1Λ2)\widetilde{A}S=\left(\begin{array}[]{ll}\Lambda_{1}&B_{1}\\ B_{2}&\Lambda_{2}\end{array}\right)\left(\begin{array}[]{rr}I&0\\ G_{1}&I\end{array}\right)=\left(\begin{array}[]{ll}\Lambda_{1}+B_{1}G_{1}&B_{1}\\ B_{2}+\Lambda_{2}G_{1}&\Lambda_{2}\end{array}\right)

Therefore the under-matrix G1G_{1} satisfies G1=Λ21B2G_{1}=-\Lambda_{2}^{-1}B_{2}, hence

S=(I0Λ21B2I)S=\left(\begin{array}[]{rr}I&0\\ -\Lambda_{2}^{-1}B_{2}&I\end{array}\right)

Construction of Tt{}^{t}T

We now want to construct a matrix Tt{}^{t}T of the form:

Tt=(IG20I){}^{t}T=\left(\begin{array}[]{rr}I&G_{2}\\ 0&I\end{array}\right)

and such that TtA~S{}^{t}T\widetilde{A}S is bloc diagonal. We have

A^=TtA~S=(Λ1+B1G1B1+G2Λ20Λ2)\hat{A}={}^{t}T\tilde{A}S=\left(\begin{array}[]{rr}\Lambda_{1}+B_{1}G_{1}&B_{1}+G_{2}\Lambda_{2}\\ 0&\Lambda_{2}\end{array}\right)

and then G2G_{2} must satisfy B1+G2Λ2=0B_{1}+G_{2}\Lambda_{2}=0.
Thus

Tt=(IB1Λ210I){}^{t}T=\left(\begin{array}[]{rr}I&-B_{1}\Lambda_{2}^{-1}\\ 0&I\end{array}\right)

and then A^\hat{A} can be written in the form

A^=TtA~S=(Λ1+B1G100Λ2)\hat{A}={}^{t}T\widetilde{A}S=\left(\begin{array}[]{rr}\Lambda_{1}+B_{1}G_{1}&0\\ 0&\Lambda_{2}\end{array}\right)

We note that since the linear system is non-symmetric, these IUs lead to a non-symmetric hierarchical preconditioner.

The first diagonal bloc of A^\hat{A} is still tridiagonal and we can repeat recursively the reduction procedure described above by using dd levels of IUs.

2 A Priori Estimates

We can obtain (see [1])a priori energy estimates that show that the (induced) IUs are indeed small structures as expected.

We multiply (2) by δiui,δi\delta_{i}u_{i},\delta_{i} to be fixed later, and summing these expressions on all indices ii, we obtain

i=12N1βiδi(uiui1)uii=12N1αiδi(ui+1ui)ui=i=12N1hfiδiui\sum_{i=1}^{2N-1}\beta_{i}\delta_{i}\left(u_{i}-u_{i-1}\right)u_{i}-\sum_{i=1}^{2N-1}\alpha_{i}\delta_{i}\left(u_{i+1}-u_{i}\right)u_{i}=\sum_{i=1}^{2N-1}hf_{i}\delta_{i}u_{i}

Hence,

i=02N1(βi+1δi+1ui+1αiδiui)(ui+1ui)=i=12N1hfiδiui.\sum_{i=0}^{2N-1}\left(\beta_{i+1}\delta_{i+1}u_{i+1}-\alpha_{i}\delta_{i}u_{i}\right)\left(u_{i+1}-u_{i}\right)=\sum_{i=1}^{2N-1}hf_{i}\delta_{i}u_{i}. (4)

We now choose δi\delta_{i} such that

δ0=1,δi+1=αiβi+1δi\delta_{0}=1,\delta_{i+1}=\frac{\alpha_{i}}{\beta_{i+1}}\delta_{i}

and we make the following hypothesis [1]:
Hypothesis(H)
i. There exists two strictly positive real numbers α\alpha and β\beta such that

αδiβi.\alpha\leq\delta_{i}\leq\beta\forall i.

ii. αi\alpha_{i} and βi\beta_{i} are strictly positive numbers, i=0,,2N\forall i=0,\ldots,2N.
iii. There exists a strictly positive real number γ\gamma such that

αiγhi\alpha_{i}\geq\frac{\gamma}{h}\forall i

We shall see later that (H) is not too restrictive for practical cases.
From (4) we infer

i=02N1αiδi(ui+1ui)2=i=12N1hfiδiui\sum_{i=0}^{2N-1}\alpha_{i}\delta_{i}\left(u_{i+1}-u_{i}\right)^{2}=\sum_{i=1}^{2N-1}hf_{i}\delta_{i}u_{i} (5)

We recall the discrete Poincaré inequality (see e.g. [2]):
Lemma 1 Let ui,i=0,,2Nu_{i},i=0,\ldots,2N be a sequence of real numbers such that u0=u2N=0u_{0}=u_{2N}=0. Then we have

i=12N1hui2i=02N11h(ui+1ui)2\sum_{i=1}^{2N-1}hu_{i}^{2}\leq\sum_{i=0}^{2N-1}\frac{1}{h}\left(u_{i+1}-u_{i}\right)^{2}

Using (H), Cauchy-Schwartz inequality and Poincaré inequality, we bound the right-hand side of (5) by

i=12N1hfiδiuiβi=12N1hfiuiβi=02N1hui2i=12N1hfi2βi=02N11h(ui+1ui)2i=12N1hfi2\sum_{i=1}^{2N-1}hf_{i}\delta_{i}u_{i}\leq\beta\sum_{i=1}^{2N-1}hf_{i}u_{i}\leq\beta\sqrt{\sum_{i=0}^{2N-1}hu_{i}^{2}}\sqrt{\sum_{i=1}^{2N-1}hf_{i}^{2}}\leq\beta\sqrt{\sum_{i=0}^{2N-1}\frac{1}{h}\left(u_{i+1}-u_{i}\right)^{2}}\sqrt{\sum_{i=1}^{2N-1}hf_{i}^{2}} (6)

Also, using (H) we can bound the left-hand side of (5) by

i=02N1αiδi(ui+1ui)2αi=02N1αi(ui+1ui)2αγi=02N11h(ui+1ui)2.\sum_{i=0}^{2N-1}\alpha_{i}\delta_{i}\left(u_{i+1}-u_{i}\right)^{2}\geq\alpha\sum_{i=0}^{2N-1}\alpha_{i}\left(u_{i+1}-u_{i}\right)^{2}\geq\alpha\gamma\sum_{i=0}^{2N-1}\frac{1}{h}\left(u_{i+1}-u_{i}\right)^{2}. (7)

From (6) and (7) we have:

i=02N11h(ui+1ui)2β2α2γ2i=12N1hfi2\sum_{i=0}^{2N-1}\frac{1}{h}\left(u_{i+1}-u_{i}\right)^{2}\leq\frac{\beta^{2}}{\alpha^{2}\gamma^{2}}\sum_{i=1}^{2N-1}hf_{i}^{2} (8)

We now introduce the incremental unknowns defined by (3). Since

i=02N11h(ui+1ui)2=i=0N11h(u2i+1u2i)2+i=0N11h(u2i+2u2i+1)2\sum_{i=0}^{2N-1}\frac{1}{h}\left(u_{i+1}-u_{i}\right)^{2}=\sum_{i=0}^{N-1}\frac{1}{h}\left(u_{2i+1}-u_{2i}\right)^{2}+\sum_{i=0}^{N-1}\frac{1}{h}\left(u_{2i+2}-u_{2i+1}\right)^{2}

we find, setting y2i=u2iy_{2i}=u_{2i}

i=0N11h(z2i+1+1α2i+1+β2i+1(β2i+1y2i+α2i+1y2i+2)y2i)2\sum_{i=0}^{N-1}\frac{1}{h}\left(z_{2i+1}+\frac{1}{\alpha_{2i+1}+\beta_{2i+1}}\left(\beta_{2i+1}y_{2i}+\alpha_{2i+1}y_{2i+2}\right)-y_{2i}\right)^{2}
+i=0N11h(y2i+2z2i+11α2i+1+β2i+1(β2i+1y2i+α2i+1y2i+2))2β2α2γ2i=12N1hfi2+\sum_{i=0}^{N-1}\frac{1}{h}\left(y_{2i+2}-z_{2i+1}-\frac{1}{\alpha_{2i+1}+\beta_{2i+1}}\left(\beta_{2i+1}y_{2i}+\alpha_{2i+1}y_{2i+2}\right)\right)^{2}\leq\frac{\beta^{2}}{\alpha^{2}\gamma^{2}}\sum_{i=1}^{2N-1}hf_{i}^{2}

which yields

i=0N1(z2i+1+α2i+1α2i+1+β2i+1(y2i+2y2i))2++i=0N1(z2i+1β2i+1α2i+1+β2i+1(y2i+2y2i))2β2α2γ2i=12N1h2fi2\begin{gathered}\sum_{i=0}^{N-1}\left(z_{2i+1}+\frac{\alpha_{2i+1}}{\alpha_{2i+1}+\beta_{2i+1}}\left(y_{2i+2}-y_{2i}\right)\right)^{2}+\\ +\sum_{i=0}^{N-1}\left(z_{2i+1}-\frac{\beta_{2i+1}}{\alpha_{2i+1}+\beta_{2i+1}}\left(y_{2i+2}-y_{2i}\right)\right)^{2}\leq\frac{\beta^{2}}{\alpha^{2}\gamma^{2}}\sum_{i=1}^{2N-1}h^{2}f_{i}^{2}\end{gathered}

Replacing fi=f(xi)f_{i}=f\left(x_{i}\right) by fi=12hxi1xi+1f(x)𝑑xf_{i}=\frac{1}{2h}\int_{x_{i-1}}^{x_{i+1}}f(x)dx and using Hölder inequality we

i=0N1(z2i+1+α2i+1α2i+1+β2i+1(y2i+2y2i))2++i=0N1(z2i+1β2i+1α2i+1+β2i+1(y2i+2y2i))2β24α2γ2i=12N1[xi1xi+1f(x)𝑑x]2β24α2γ2i=12N1(xi1xi+1f2(x)𝑑x)(xi1xi+1𝑑x)=β22α2γ2hi=12N1(xi1xi+1f2(x)𝑑x)β2α2γ2hf2\begin{gathered}\sum_{i=0}^{N-1}\left(z_{2i+1}+\frac{\alpha_{2i+1}}{\alpha_{2i+1}+\beta_{2i+1}}\left(y_{2i+2}-y_{2i}\right)\right)^{2}+\\ +\sum_{i=0}^{N-1}\left(z_{2i+1}-\frac{\beta_{2i+1}}{\alpha_{2i+1}+\beta_{2i+1}}\left(y_{2i+2}-y_{2i}\right)\right)^{2}\leq\frac{\beta^{2}}{4\alpha^{2}\gamma^{2}}\sum_{i=1}^{2N-1}\left[\int_{x_{i-1}}^{x_{i+1}}f(x)dx\right]^{2}\leq\\ \leq\frac{\beta^{2}}{4\alpha^{2}\gamma^{2}}\sum_{i=1}^{2N-1}\left(\int_{x_{i-1}}^{x_{i+1}}f^{2}(x)dx\right)\left(\int_{x_{i-1}}^{x_{i+1}}dx\right)=\frac{\beta^{2}}{2\alpha^{2}\gamma^{2}}h\sum_{i=1}^{2N-1}\left(\int_{x_{i-1}}^{x_{i+1}}f^{2}(x)dx\right)\leq\frac{\beta^{2}}{\alpha^{2}\gamma^{2}}h\|f\|^{2}\end{gathered}

where f=|f|L2(0,1)\|f\|=|f|_{L^{2}(0,1)}.
We now develop the left-hand side of this inequality and we obtain:

A=2i=0N1z2i+12+i=0N1α2i+12+β2i+12(α2i+1+β2i+1)2(y2i+2y2i)2+2i=0N1α2i+1β2i+1α2i+1+β2i+1z2i+1(y2i+2y2i).A=2\sum_{i=0}^{N-1}z_{2i+1}^{2}+\sum_{i=0}^{N-1}\frac{\alpha_{2i+1}^{2}+\beta_{2i+1}^{2}}{\left(\alpha_{2i+1}+\beta_{2i+1}\right)^{2}}\left(y_{2i+2}-y_{2i}\right)^{2}+2\sum_{i=0}^{N-1}\frac{\alpha_{2i+1}-\beta_{2i+1}}{\alpha_{2i+1}+\beta_{2i+1}}z_{2i+1}\left(y_{2i+2}-y_{2i}\right).

Using Young inequality, we then find

2i=0N1α2i+1β2i+1α2i+1+β2i+1z2i+1(y2i+2y2i)εi=0N1z2i+121εi=0N1(α2i+1β2i+1)2(α2i+1+β2i+1)2(y2i+2y2i)2.2\sum_{i=0}^{N-1}\frac{\alpha_{2i+1}-\beta_{2i+1}}{\alpha_{2i+1}+\beta_{2i+1}}z_{2i+1}\left(y_{2i+2}-y_{2i}\right)\geq-\varepsilon\sum_{i=0}^{N-1}z_{2i+1}^{2}-\frac{1}{\varepsilon}\sum_{i=0}^{N-1}\frac{\left(\alpha_{2i+1}-\beta_{2i+1}\right)^{2}}{\left(\alpha_{2i+1}+\beta_{2i+1}\right)^{2}}\left(y_{2i+2}-y_{2i}\right)^{2}.

Here ε\varepsilon is a strictly positive real number which will be fixed later. Hence

A(2ε)i=0N1z2i+12+i=0N1(α2i+12+β2i+12(α2i+1+β2i+1)21ε(α2i+1β2i+1)2(α2i+1+β2i+1)2)(y2i+2y2i)2A\geq(2-\varepsilon)\sum_{i=0}^{N-1}z_{2i+1}^{2}+\sum_{i=0}^{N-1}\left(\frac{\alpha_{2i+1}^{2}+\beta_{2i+1}^{2}}{\left(\alpha_{2i+1}+\beta_{2i+1}\right)^{2}}-\frac{1}{\varepsilon}\frac{\left(\alpha_{2i+1}-\beta_{2i+1}\right)^{2}}{\left(\alpha_{2i+1}+\beta_{2i+1}\right)^{2}}\right)\left(y_{2i+2}-y_{2i}\right)^{2}

Since αi>0\alpha_{i}>0 and βi>0,i\beta_{i}>0,\forall i, we have

1(α2i+1+β2i+1)2(α2i+12+β2i+121ε(α2i+1β2i+1)2)(11ε)α2i+12+β2i+12(α2i+1+β2i+1)2.\frac{1}{\left(\alpha_{2i+1}+\beta_{2i+1}\right)^{2}}\left(\alpha_{2i+1}^{2}+\beta_{2i+1}^{2}-\frac{1}{\varepsilon}\left(\alpha_{2i+1}-\beta_{2i+1}\right)^{2}\right)\geq\left(1-\frac{1}{\varepsilon}\right)\frac{\alpha_{2i+1}^{2}+\beta_{2i+1}^{2}}{\left(\alpha_{2i+1}+\beta_{2i+1}\right)^{2}}.

We set ξi=αiβi\xi_{i}=\frac{\alpha_{i}}{\beta_{i}} and we introduce the function ϕ(ξ)=1+ξ2(1+ξ)2\phi(\xi)=\frac{1+\xi^{2}}{(1+\xi)^{2}} for ξ>0\xi>0. It is easy to check that ϕ(ξ)12,ξ>0\phi(\xi)\geq\frac{1}{2},\forall\xi>0. Hence

α2i+12+β2i+12(α2i+1+β2i+1)2=ϕ(ξ2i+1)12,i=0,,N\frac{\alpha_{2i+1}^{2}+\beta_{2i+1}^{2}}{\left(\alpha_{2i+1}+\beta_{2i+1}\right)^{2}}=\phi\left(\xi_{2i+1}\right)\geq\frac{1}{2},\forall i=0,\ldots,N

and we then have

1(α2i+1+β2i+1)2(α2i+12+β2i+121ε(α2i+1β2i+1)2)12(11ε).\frac{1}{\left(\alpha_{2i+1}+\beta_{2i+1}\right)^{2}}\left(\alpha_{2i+1}^{2}+\beta_{2i+1}^{2}-\frac{1}{\varepsilon}\left(\alpha_{2i+1}-\beta_{2i+1}\right)^{2}\right)\geq\frac{1}{2}\left(1-\frac{1}{\varepsilon}\right).

Thus, taking ε=32\varepsilon=\frac{3}{2}, we obtain, after some simplifications

A16i=0N1z2i+12+16i=0N1(y2i+2y2i)2A\geq\frac{1}{6}\sum_{i=0}^{N-1}z_{2i+1}^{2}+\frac{1}{6}\sum_{i=0}^{N-1}\left(y_{2i+2}-y_{2i}\right)^{2}

We have then the following result:
Proposition 2 [1]Under the hypothesis (H)(H), the incremental unknowns defined by (3) satisfy the following a priori estimates:

i=0N1z2i+126β2α2β2hi=0N1(y2i+2y2i)26β2α2β2h.\begin{gathered}\sum_{i=0}^{N-1}z_{2i+1}^{2}\leq\frac{6\beta^{2}}{\alpha^{2}\beta^{2}}h\\ \sum_{i=0}^{N-1}\left(y_{2i+2}-y_{2i}\right)^{2}\leq\frac{6\beta^{2}}{\alpha^{2}\beta^{2}}h.\end{gathered}

Proposition is valid for general definition of the IUs given in (3), under Hypothesis (H). In particular, the scheme used for the discretization of the convective term is not specified.
i. Centered Convection-Diffusion IUs

We have

αi=1h(1cih2),βi=1h(1+cih2).\alpha_{i}=\frac{1}{h}\left(1-\frac{c_{i}h}{2}\right),\beta_{i}=\frac{1}{h}\left(1+\frac{c_{i}h}{2}\right).

Hence, if there exists a strictly positive real number γ\gamma such that

1maxi|ci|h2γ1-\max_{i}\frac{\left|c_{i}\right|h}{2}\geq\gamma

then the asumptions i., ii. and iii. of (H) are satisfied. This condition can be satisfied by taking hh small enough.
ii. Uncentered Convection-Diffusion IUs

  • If ci0c_{i}\geq 0 then

αi=αi+=1h,βi=βi+=1h+ci.\alpha_{i}=\alpha_{i}^{+}=\frac{1}{h},\quad\beta_{i}=\beta_{i}^{+}=\frac{1}{h}+c_{i}.

Hence, αi\alpha_{i} and βi\beta_{i} are stricly positive real numbers.

  • If ci0c_{i}\leq 0 then

αi=αi=1hci,βi=βi=1h\alpha_{i}=\alpha_{i}^{-}=\frac{1}{h}-c_{i},\quad\beta_{i}=\beta_{i}^{-}=\frac{1}{h}

Here again, αi\alpha_{i} and βi\beta_{i} are strictly positive real numbers.
In conclusion, and without any asumption, the asumptions i., ii. and iii. of (H) are satisfied.

3 High-order compact schemes (HOC)

The real advantage of (HOC) lies not in the increased accuracy (although this is sometimes important), but rather in the fact that problems requiring fine grid can be solved on coarse grids. This implies that less memory and thus less time is required to solve the same problem to the same accuracy. Since "time is money", (HOC) schemes directly reduce the expense of approximating a differential equation numerically.

High-order compact (HOC) schemes of the type studied here increase the accuracy of the standard central difference approximation from O(h2)O\left(h^{2}\right) to O(h4)O\left(h^{4}\right) by including compact approximation to the leading truncation error terms. The idea is to operate on the differential equation as an auxiliary relation to obtain expressions for higher-order derivatives in the truncation error.

We define δxnui,n=1,2\delta_{x}^{n}u_{i},n=1,2 to be the standard central difference operator for the nn - th derivative of uu at point ii on a uniform grid of mesh size hh,

δxui=ui+1ui12h,δx2ui=ui+12ui+ui1h2\delta_{x}u_{i}=\frac{u_{i+1}-u_{i-1}}{2h},\delta_{x}^{2}u_{i}=\frac{u_{i+1}-2u_{i}+u_{i-1}}{h^{2}}

Central differencing (1) yields

δx2ui+ciδxuiτi=fi-\delta_{x}^{2}u_{i}+c_{i}\delta_{x}u_{i}-\tau_{i}=f_{i} (9)

where τi\tau_{i} is the local truncation error at node ii,

τi=h212[2cd3udx3d4udx4]i+O(h4)\tau_{i}=\frac{h^{2}}{12}\left[2c\frac{d^{3}u}{dx^{3}}-\frac{d^{4}u}{dx^{4}}\right]_{i}+O\left(h^{4}\right) (10)

We seek to approximate the leading term in (10) and include it in the difference formulation to yield an O(h4)O\left(h^{4}\right) method. Assuming the solution is sufficiently regular, we may accomplish this by differentiating (1) to yield

d3udx3|i=[cd2udx2+dcdxdudxdfdx]i\left.\frac{d^{3}u}{dx^{3}}\right|_{i}=\left[c\frac{d^{2}u}{dx^{2}}+\frac{dc}{dx}\frac{du}{dx}-\frac{df}{dx}\right]_{i}

which can be approximated compactly as

d3udx3|i=ciδx2ui+δxciδxuiδxfi+O(h2)\left.\frac{d^{3}u}{dx^{3}}\right|_{i}=c_{i}\delta_{x}^{2}u_{i}+\delta_{x}c_{i}\delta_{x}u_{i}-\delta_{x}f_{i}+O\left(h^{2}\right) (11)

and similarly

d4udx4|i=[cd3udx3+2dcdxd2udx2+d2cdx2dudxd2fdx2]i\displaystyle\left.\frac{d^{4}u}{dx^{4}}\right|_{i}=\left[c\frac{d^{3}u}{dx^{3}}+2\frac{dc}{dx}\frac{d^{2}u}{dx^{2}}+\frac{d^{2}c}{dx^{2}}\frac{du}{dx}-\frac{d^{2}f}{dx^{2}}\right]_{i} (12)
=cid3udx3|i+2δxciδx2ui+δx2ciδxciδx2fi+O(h2)\displaystyle=\left.c_{i}\frac{d^{3}u}{dx^{3}}\right|_{i}+2\delta_{x}c_{i}\delta_{x}^{2}u_{i}+\delta_{x}^{2}c_{i}\delta_{x}c_{i}-\delta_{x}^{2}f_{i}+O\left(h^{2}\right)

Relations (11) and (12) can be combined with (10) to yield the new truncation error expression:

τi=h212[(ci22δxci)δx2ui+(ciδxciδx2ci)δxuiciδxfi+δx2fi]+O(h4)\tau_{i}=\frac{h^{2}}{12}\left[\left(c_{i}^{2}-2\delta_{x}c_{i}\right)\delta_{x}^{2}u_{i}+\left(c_{i}\delta_{x}c_{i}-\delta_{x}^{2}c_{i}\right)\delta_{x}u_{i}-c_{i}\delta_{x}f_{i}+\delta_{x}^{2}f_{i}\right]+O\left(h^{4}\right)

which can be use to increase the accuracy of the approximation scheme. The resulting high-order compact scheme is

Aiδx2ui+Ciδxui=Fi+O(h4)-A_{i}\delta_{x}^{2}u_{i}+C_{i}\delta_{x}u_{i}=F_{i}+O\left(h^{4}\right) (13)

where

Ai=1+h212(ci22δxci)\displaystyle A_{i}=1+\frac{h^{2}}{12}\left(c_{i}^{2}-2\delta_{x}c_{i}\right) (14)
Ci=ci+h212(δx2ciciδxci)\displaystyle C_{i}=c_{i}+\frac{h^{2}}{12}\left(\delta_{x}^{2}c_{i}-c_{i}\delta_{x}c_{i}\right) (15)
Fi=fi+h212(δx2ficiδxfi)\displaystyle F_{i}=f_{i}+\frac{h^{2}}{12}\left(\delta_{x}^{2}f_{i}-c_{i}\delta_{x}f_{i}\right) (16)

We can write (13) in the form (2) where:

αi=AihCi2,βi=Aih+Ci2,fi:=Fi for i=1,,2N\alpha_{i}=\frac{A_{i}}{h}-\frac{C_{i}}{2},\beta_{i}=\frac{A_{i}}{h}+\frac{C_{i}}{2},f_{i}:=F_{i}\text{ for }i=1,\ldots,2N

and we can define IU in this case. If we consider 1D convection diffusion problem with constant coefficients then:

Ai=1+h2c212,Ci=c and Fi=fi for i=1,,2NA_{i}=1+\frac{h^{2}c^{2}}{12},C_{i}=c\text{ and }F_{i}=f_{i}\text{ for }i=1,\ldots,2N

In this case we can see that, without any asumptions, the hypothesis (H) are satisfied.

4 Condition Number Analysis

The condition number is defined by

κ=max|λk|min|λk|\kappa=\frac{\max\left|\lambda_{k}\right|}{\min\left|\lambda_{k}\right|} (17)

where {λk}\left\{\lambda_{k}\right\} is the set of eigenvalues of the real matrix we wish to solve.
The tridiagonal matrices generated by compact 1D numerical schemes lend themselves to theoretical eigenvalue analysis, and our hope is that we can draw some conclusions from such analysis that might apply to higher dimensions.

We will consider 1D convection diffusion problem with constant coefficients. For a tridiagonal Toeplitz matrix which has the form

TN(c,d,e)=(decdecd)T_{N}(c,d,e)=\left(\begin{array}[]{cccccc}d&e&&&&\\ c&d&\cdot&&&\\ &\cdot&\cdot&\cdot&&\\ &&\cdot&\cdot&\cdot&\\ &&&\cdot&\cdot&e\\ &&&&c&d\end{array}\right)

the eigenvalues are known explicitly (see [5]),

d+2(ce)1/2cos(kπ/(N+1)),k=1,,Nd+2(ce)^{1/2}\cos(k\pi/(N+1)),k=1,\ldots,N

Using this result, we can obtain condition number for each of the 3 schemes (see [12]):

κHOC=c2h2+12+(c2h2+6ch+12)(c2h26ch+12)c2h2+12(c2h2+6ch+12)(c2h26ch+12),\displaystyle\kappa_{HOC}=\frac{c^{2}h^{2}+12+\sqrt{\left(c^{2}h^{2}+6ch+12\right)\left(c^{2}h^{2}-6ch+12\right)}}{c^{2}h^{2}+12-\sqrt{\left(c^{2}h^{2}+6ch+12\right)\left(c^{2}h^{2}-6ch+12\right)}}, (18)
κCDS={2+4c2h22+4c2h2,|ch|2|c|2,|ch|>2,\displaystyle\kappa_{CDS}=\left\{\begin{array}[]{c}\frac{2+\sqrt{4-c^{2}h^{2}}}{2+\sqrt{4-c^{2}h^{2}}},|ch|\leq 2\\ \frac{|c|}{2},|ch|>2\end{array},\right. (19)
κCDS={2+4c2h22+4c2h2,|ch|2|c|2,|ch|>2,\displaystyle\kappa_{CDS}=\left\{\begin{array}[]{c}\frac{2+\sqrt{4-c^{2}h^{2}}}{2+\sqrt{4-c^{2}h^{2}}},|ch|\leq 2\\ \frac{|c|}{2},|ch|>2\end{array},\right. (19)
κUDS=|ch|+2+|ch|+2|ch|+2|ch|+2\displaystyle\kappa_{UDS}=\frac{|ch|+2+\sqrt{|ch|+2}}{|ch|+2-\sqrt{|ch|+2}} (20)

For HOC scheme, as ch0ch\rightarrow 0, the condition number behaves as c2h2c^{-2}h^{-2} and for large chch, the condition number behaves as c2h2c^{2}h^{2}. This means that for large NN and |ch||ch| very large or very small, the condition number is extremely large and may pose some difficulty for an iterative solver. This is one instance where CDS and UDS appear to compare favorably with HOC, except of course for the fact that for ch>2ch>2, CDS is oscillatory and for large chch, the UDS models an overly diffusive problem.

Introducing IU, after using all the levels (what means LL if h=12Lh=\frac{1}{2^{L}} ) we will obtain a diagonal matrix which have the elements,

α1=c2h2+126h2,αn=(c2h2+6ch+1212h2)2n1+(c2h26ch+1212h2)2n1i=0n2[(c2h2+6ch+1212h2)2i+(c2h26ch+1212h2)2i],n=2,,L,\alpha_{1}=\frac{c^{2}h^{2}+12}{6h^{2}},\alpha_{n}=\frac{\left(\frac{c^{2}h^{2}+6ch+12}{12h^{2}}\right)^{2^{n-1}}+\left(\frac{c^{2}h^{2}-6ch+12}{12h^{2}}\right)^{2^{n-1}}}{\prod_{i=0}^{n-2}\left[\left(\frac{c^{2}h^{2}+6ch+12}{12h^{2}}\right)^{2^{i}}+\left(\frac{c^{2}h^{2}-6ch+12}{12h^{2}}\right)^{2^{i}}\right]},n=2,\ldots,L,

and α1>α2>>αL\alpha_{1}>\alpha_{2}>\ldots>\alpha_{L}. In conclusion, the condition number will be:

κIUHOC=2(12+c2h2)i=0L2[(c2h2+6ch+12)2i+(c2h26ch+12)2i](c2h2+6ch+12)2L1+(c2h26ch+12)2L1\kappa_{IUHOC}=2\frac{\left(12+c^{2}h^{2}\right)\prod_{i=0}^{L-2}\left[\left(c^{2}h^{2}+6ch+12\right)^{2^{i}}+\left(c^{2}h^{2}-6ch+12\right)^{2^{i}}\right]}{\left(c^{2}h^{2}+6ch+12\right)^{2^{L-1}}+\left(c^{2}h^{2}-6ch+12\right)^{2^{L-1}}} (21)

In the next table we have a comparison between condition number of HOC scheme and condition number of HOC scheme preconditioned by IU method.

ch=2\mathrm{ch}=2 ch=5\mathrm{ch}=5 ch =10=10 ch =20=20
h c k c k c k c k
18\frac{1}{8} 16 4.1426 40 3.3550 80 8.0960 160 16.215
1.5960 1.4133 2.5549 4.5692
1.3322 1.2330 1.8357 2.8620
116\frac{1}{16} 32 4.6936 80 3.6955 160 10.639 320 31.394
1.7267 1.4915 3.1834 8.3567
1.3718 1.2524 2.1057 4.6964
1.3333 1.2333 1.8664 3.3772
132\frac{1}{32} 64 4.8523 160 3.7903 320 11.526 640 40.698
1.7646 1.5135 3.4032 10.680
1.3844 1.2584 2.2051 5.8448
1.3341 1.2335 1.8887 3.8699
1.3333 1.2333 1.8667 3.4329

A more acurate bound for kIUHOCk_{IUHOC} can be given as follows:
Proposition 3 Let kIUHOCk_{IUHOC} be the condition number given by (21), where h=1/2L(LN)h=1/2^{L}(L\in N). Then

kIUHOC2L1=12h,k_{IUHOC}\leqslant 2^{L-1}=\frac{1}{2h},

for each L=2,3,L=2,3,\ldots.
Proof. Let

κIUHOC=2(12+c2h2)i=0L2[(c2h2+6ch+12)2i+(c2h26ch+12)2i](c2h2+6ch+12)2L1+(c2h26ch+12)2L1\kappa_{IUHOC}=2\frac{\left(12+c^{2}h^{2}\right)\prod_{i=0}^{L-2}\left[\left(c^{2}h^{2}+6ch+12\right)^{2^{i}}+\left(c^{2}h^{2}-6ch+12\right)^{2^{i}}\right]}{\left(c^{2}h^{2}+6ch+12\right)^{2^{L-1}}+\left(c^{2}h^{2}-6ch+12\right)^{2^{L-1}}}

be the condition number of HOC scheme preconditioned by IU method. Let us denote

A=c2h2+6ch+12,B=c2h26ch+12.A=c^{2}h^{2}+6ch+12,\quad B=c^{2}h^{2}-6ch+12.

Then A+B=2(12+c2h2)A+B=2\left(12+c^{2}h^{2}\right) and with this notations we shall have
kIUHOC=(A+B)i=0L2(A2i+B2i)A2L1+B2L1=(A+B)2(A2+B2)(A4+B4)(A2L2+B2L2)A2L1+B2L1k_{IUHOC}=\frac{(A+B)\prod_{i=0}^{L-2}\left(A^{2^{i}}+B^{2^{i}}\right)}{A^{2^{L-1}}+B^{2^{L-1}}}=\frac{(A+B)^{2}\left(A^{2}+B^{2}\right)\left(A^{4}+B^{4}\right)\ldots\left(A^{2^{L-2}}+B^{2^{L-2}}\right)}{A^{2^{L-1}}+B^{2^{L-1}}}.
Now, using the inequality (x+y)22(x2+y2)(x+y)^{2}\leqslant 2\left(x^{2}+y^{2}\right), for real numbers sequentially, we get

kIUHOC\displaystyle k_{IUHOC} 2(A2+B2)2(A4+B4)(A2L2+B2L2)A2L1+B2L122(A4+B4)2(A2L2+B2L2)A2L1+B2L1\displaystyle\leqslant\frac{2\left(A^{2}+B^{2}\right)^{2}\left(A^{4}+B^{4}\right)\ldots\left(A^{2^{L-2}}+B^{2^{L-2}}\right)}{A^{2^{L-1}}+B^{2^{L-1}}}\leqslant\frac{2^{2}\left(A^{4}+B^{4}\right)^{2}\ldots\left(A^{2^{L-2}}+B^{2^{L-2}}\right)}{A^{2^{L-1}}+B^{2^{L-1}}}\leqslant\ldots
2L2(A2L2+B2L2)2A2L1+B2L12L1(A2L1+B2L1)A2L1+B2L1=2L1\displaystyle\ldots\leqslant\frac{2^{L-2}\left(A^{2^{L-2}}+B^{2^{L-2}}\right)^{2}}{A^{2^{L-1}}+B^{2^{L-1}}}\leqslant\frac{2^{L-1}\left(A^{2^{L-1}}+B^{2^{L-1}}\right)}{A^{2^{L-1}}+B^{2^{L-1}}}=2^{L-1}

Remark 4 We notice that the upper bound of the condition number does not depent of the equation’s coefficient cc. From this point of view we can say that we have a uniform condition number estimate.

References

[1] Chebab, J.P., Miranville, A., Induced Hierarchical Preconditioners: The finite difference Case, Publication ANO-371 (1997).
[2] Chen, M., Temam, R., Incremental Unknowns for solving Partial Differential Equations, Numer. Math. 59 (1991), 255-271.
[3] Chen, M., Temam, R., Incremental Unknowns in Finite Differences: Condition Number of the Matrix, Siam J. on Matrix Analysis and Applications (SIMAX) 14 (1993) 2, 432-455.
[4] Chen, M., Temam, R., Incremental Unknowns for Solving Convection Diffusion Equations, (1993).
[5] Elman, H.C., Golub, G.H., Iterative Methods for Cyclically-Reduced Non-Self-Adjoint Linear Systems, Math. Comp. 54 (1990), 671-700.
[6] Kreiss,H., Oliger, J., Methods for the Approximate Solution of Time Dependent Problems, GARP Report No.10, 1973.
[7] Lele, S., Compact Finite Difference Schemes with Spectral-like Resolution, J. Comp. Phys. 103 (1992) 1, 16-42.
[8] MacKinnon,R.J., Carey,G.F., Superconvergent Derivatives: A Taylor Series Analysis, International Journal for Numerical Methods in Engineering 28 (1989), 489-509.
[9] MacKinnon,R.J., Carey,G.F., Nodal Superconvergence and Solution Enhancement for a class of finite element and finite difference methods, SIAM Journal on Scientific and Statistical Computing 11 (1990) 2, 343-353.
[10] MacKinnon,R.J., Johnson, R.W., Differential equation based representation of truncation errors for accurate numerical simulation, International Journal for Numerical Methods in Fluids 13 (1991), 739-757.
[11] MacKinnon,R.J., Langerman,M.A., A compact high-order finite-element method for elliptic transport problems with variable coefficients, Numerical Methods for Partial Diff. Equations 10 (1994), 1-19.
[12] Spotz,W.F., High-order compact finite difference schemes for computational mechanics, Ph.D. Dissertation, the University of Texas at Austin, December, 1995.
[13] Temam, R., Inertial Manifolds and Multigrid Methods, SIAM J. Math. Anal. 21 (1990) 1, 154-178.
"T.Popoviciu" Institute of Numerical Analysis P.O. Box 68

3400 Cluj-Napoca 1
e-mail: amuresan@ictp-acad.math.ubbcluj.ro
acurseu@ictp-acad.math.ubbcluj.ro

2001

Related Posts