Construction of consistent population functions

In this section we show by an explicit construction, that functions being consistent to the lattice-Boltzmann equation of arbitrarily high order do exist if we assume sufficient regularity of the solution to the target problem. In order to avoid the impression, that these consistent functions appear out of hot air we point out that the construction is strongly motivated by the results of the formal analysis. In fact, these approximate population functions are obtained by truncating the asymptotic expansion.

Before we formulate the smoothness requirements and describe the construction process, we state the following preparatory lemma.

Lemma 4.4. Let f ∈Cⁿ⁺¹([0, T]×R/LR,R) and 0≤ |h| ≤1, then:

The truncated Taylor series is split into two sums, where the indices β, γ are sup-posed to be non-negative. The second sum is abbreviated by S_n^∗.

S_n = X addends. However, we estimateS_n^∗simply by dropping the restriction 2β+γ ≥n+1 and adding up all derivatives of order ≤n.

|S_n^∗| ≤ |h|ⁿ⁺¹ X

4.3. Construction of consistent population functions 179

Estimating R_n yields:

|Rn| ≤ |h|ⁿ⁺¹ X

β+γ=n+1

|∂_t^β∂_x^γf(θ, ξ)|

≤ |h|ⁿ⁺¹ X

β+γ=n+1

k∂_t^β∂_x^γfkC⁰ = |h|ⁿ⁺¹ |f|Cⁿ⁺¹

Since

f(t+h², x+h) − Xn k=0

h^kT_k(∂_t, ∂_x)f(t, x) = S_n^∗ + R_n

the assertion follows by taking the modulus, applying the triangle inequality on the right hand side, using the estimates for R_n and S_n^∗ and observing, that kfkCⁿ⁺¹ =

kfkCⁿ + |f|Cⁿ⁺¹.

Existence and Smoothness Hypothesis: Let α ∈ N₀ be given. Consider the target problem with the initial and source functions

w∈C^2α+4(R/LR,R) g∈C^2α+2([0, T]×R/LR,R). (4.31) We assume, that there exist functionsu⁽⁰⁾, u⁽²⁾, ..., u^(2α) with

u^(2β) ∈C^2α−2β+4([0, T]×R/LR,R) β∈ {0,1, ..., α} (4.32) being solutions to the following coupled initial value problems (cf. eqn. (4.19)).

u^(2β)(0,·) = w δ_0β (4.33)

∂_tu^(2β)+ a∂_x−¹_θ(τ −¹₂)

| {z }

=ν

∂_x²+c

u^(2β) = EXT^(2β)+ LOW^(2β) (4.34)

Here δ_0β denotes the Kronecker-delta. The two source terms are given by:

LOW^(2α) :=

α−1

X

β=0

h

ωhP_2(α−β+1)i + ωaθhsP_2(α−β)+1i − chP_2(α−β)ii u^(2β) EXT^(2α) := hP2αig

We suppose, that the algorithmic parametersθ, τ are related to the diffusivity ν of the target problem according to (4.7). Then u⁽⁰⁾ is just the solution of the target problem (4.1). Furthermore let us define

u^(2α+2) ≡0, u^(2β+1) ≡0 β ∈ {0,1, ..., α+ 1}. For ease of notation we set u⁽⁻¹⁾≡0 andu⁽⁻²⁾ ≡0.

Remark 4.4. The smoothness hypothesis is motivated by the pragmatic reason, to obtain the main result (theorem 4.2) in classical sense. Since the polynomials defined by (4.15) mix time and space derivatives, we do not want to distinguish between temporal and spatial regularity. Apart from this, however, we do not

assume more regularity than necessary. In particular, the smoothness hypothesis guarantees, that the terms on the right hand side of (4.34) keep their classical meaning. It has to be checked separately, whether the required regularity of the data really implies the assumed regularity of the u^(2β)’s by successively solving the initial value problem (4.33),(4.34). Using the regularizing effect of the diffusive term (second derivative), the assumption may be relaxed; but we do not enter this issue here.

Due to the regularity assumption concerninggand theu^(2β)’s we can construct the following functions f⁽ⁿ⁾ forn∈ {0,1, ...,2α+ 3}

f⁽ⁿ⁾ :=

Xn k=0

P_n−k(τ, ∂t, S∂x)E^(k) + τ Pn−2gw, (4.35) whereE^(k) abbreviates D(u^(k)) +A(u^(k−1)) +R(u^(k−2)). It should be noticed, that the degree of the polynomialP_n−kisn−k. So, the highest (possibly mixed) deriva-tive contained in P_n−k(τ, ∂t, S∂x) is of order n−k≤2α+ 3−k. Since it is applied to u^(k)∈C^2α−k+4,u^(k−1) ∈C^2α−k+5 andu^(k−2) ∈C^2α−k+6 the sum is well defined.

Likewise, as we have assumed g ∈ C^2α+2, the highest derivative acting on g and being of order 2α+ 1, also exists in the classical sense. Actually, we have always presupposed one order of regularity more than necessary for the definition of the f⁽ⁿ⁾’s. This additional order is required to derive the estimate of lemma 4.6.

The approximate population function f_ǫ to be compared with the output of the lattice-Boltzmann algorithm is assembled by the following truncated asymptotic expansion

f_ǫ :=

2α+3X

n=0

f⁽ⁿ⁾ǫⁿ, (4.36)

where the f⁽ⁿ⁾’s appear as asymptotic orders. By virtue of their construction (cf.

4.16 resp. (4.35)) it is no surprise, that they exhibit the following algebraic property becoming crucial further below.

Lemma 4.5. u⁽ⁿ⁾ can be recovered from f⁽ⁿ⁾ (n = 0, ...,2α + 3) by computing its zeroth moment.

X

s∈S

f_s⁽ⁿ⁾=u⁽ⁿ⁾





 P

s∈Sfs^(2β) = u^(2β) for β ∈ {0,1, ..., α} P

s∈Sf_s^(2β+1) = 0 for β ∈ {0,1, ..., α+ 1} P

s∈Sf_s^(2α+2) = 0

Proof: For n = 0,1 the assertion can be directly verified using (4.13) and (4.14);

forn≥2 we take recourse to theorem 4.2. From the representation formula (4.35), we infer for even n = 2β + 2 and β ∈ {0,1, ..., α} by repeating the arithmetic manipulations performed already in the proof of theorem 4.2:

X

s∈S

f^(2β+2)=u^(2β+2)− ∂_t+a∂_x−¹_θ(τ −¹₂)∂_x²+c

u^(2β)+ EXT^(2β)+ LOW^(2β)

| {z }

=0

4.3. Construction of consistent population functions 181

Since u^(2β) satisfies (4.34), the underbraced term vanishes and the assertion is ob-tained. In the same way we get for odd n = 2β + 3 and β ∈ {0,1, ..., α} the expression¹³,

X

s∈S

f^(2β+3)=u^(2β+3)− ∂_t+a∂_x−¹_θ(τ −¹₂)∂_x²+c

u^(2β+1)+ LOW^(2β+1)

| {z }

=0

which vanishes sinceu⁽¹⁾=u⁽³⁾=...=u^(2α+3)= 0, thus concluding the proof.

Therefore we can write the zeroth moment off_ǫ u_ǫ := X

s∈S

f_ǫ,s = Xα β=0

u^(2β)ǫ^2β

as a truncated asymptotic expansion with theu^(2β)’s as coefficient functions. Finally let us define the following remainder terms,

r⁽ⁿ⁾_ǫ,s(t, x) := f_s⁽ⁿ⁾(t+ǫ², x+sǫ)−

2α+3X−n k=0

ǫ^k T_k(∂_t, s∂_x)f_s⁽ⁿ⁾(t, x) (4.37) satisfying the following estimate.

Lemma 4.6. There exists a constant K =K(a, c, τ, θ) depending only on the pa-rameters of the target problem¹⁴, such that

2α+3X

n=0

kr_ǫ,s⁽ⁿ⁾kC⁰ ǫⁿ ≤ KX^α

γ=0

ku^(2γ)kC^2(α−γ)+4 + kgkC^2α+2

ǫ^2α+4

for alls∈ S and 0≤ǫ≤1.

Proof: In order to estimate rǫ,s^(2β) we apply lemma 4.4 tofs^(2β). kr^(2β)_ǫ,s kC⁰ = sup

x∈R/LR t∈[0,T−ǫ²]

f_s⁽ⁿ⁾(t+ǫ², x+sǫ)−

2α+3X−n k=0

ǫ^kT_k(∂_t, s∂_x)f_s⁽ⁿ⁾(t, x)

≤ kf_s⁽ⁿ⁾kC^2α+4−2β ǫ^2α+4−2β

Using the notation of (4.8), equation (4.35) is written forn= 2β withβ∈ {0, ..., α+

1} more explicitly in the following form:

f^(2β) = Xβ γ=0

h

P_2β−2γD^(2γ) + P_{2β−2γ−1}A^(2γ) + P_{2β−2γ−2}R^(2γ)i

+ τ P_2β−2gw

13Note, that LOW^(2β+1) is a linear combination of derivatives ofu⁽¹⁾, ..., u^(2β−1). For the exact definition see equation (4.20). In particular, LOW^(2β+1) vanishes ifu⁽¹⁾, ..., u^(2β−1)vanish.

14Recall: aadvection speed,creactivity,τ collision time,θweight parameter.

Using the triangle inequality, this leads to:

kf_s⁽ⁿ⁾kC^2α+4−2β ≤nX^β

γ=0

P_2β−2γDs^(2γ)+P_{2β−2γ−1}A^(2γ)s +P_{2β−2γ−2}R^(2γ)s

C^2α+4−2β

+ τ P_2β−2gw_s

C^2α+4−2β

o

(4.38) Let us now simplify the right hand side by estimating it further. Since P_2β−2g contains (mixed) derivatives ofg up to order 2β−2, we see, thatkP_2β−2gkC^2α+4−2β

involves derivatives up to order (2β−2) + (2α+ 4−2β) = 2α+ 2. Therefore we get kτ P_2β−2gw_skC^2α+4−2β ≤ G^(2β)kgkC^2α+2

with G^(2β):=τ sum of absolute values of polynomial coefficients of P_2β−2(τ,·,·) . Similarly, we proceed with the first addend. Introducing the polynomials p^(2β)_s,2γ, whose coefficients depend only on the parameters a, c, τ and θ via Ds(1), As(1), Rs(1),

p^(2β)_s,2γ(ϑ, ς) :=

D^s(1)P_2(β−γ)+A^s(1)P_{2(β−γ)−1}+R^s(1)P_{2(β−γ)−2}

(τ, ϑ, ς) we define the constants K_2γ^(2β) analogously to the G^(2β)’s by:

K_2γ^(2β) := sum of absolute values of polynomial coefficients of p^(2β)_⊕,2γ

So, the K_2γ^(2β)’s depend also exclusively on a, c, τ and θ. In the representation of f^(2β) derivatives of u^(2γ) occur up to order 2β −2γ (for γ ≤ β) contained in P_2β−2γ =P_2β−2γ(τ, ∂_t, s∂_x). That is why the highest derivative of u^(2γ) occurring on the right hand side of (4.38) is of order (2α+ 4−2β) + (2β−2γ) = 2α+ 4−2γ.

Hence we obtain:

P_2β−2γDs^(2γ)+P_{2β−2γ−1}A^(2γ)s +P_{2β−2γ−2}R^(2γ)s

C^2α+4−2β ≤ K_2γ^(2β)ku^(2γ)kC^2α+4−2γ

So we conclude:

kr^(2β)_ǫ,s kC⁰ ≤ X^β

γ=0

K_2γ^(2β)ku^(2γ)kC^2α+4−2γ + G^(2β) kgkC^2α+2

ǫ^2α+4−2β

The casen= 2β+ 1 with β ∈ {0, ..., α+ 1} is treated in the same way, resulting in an analogous inequality:

kr^(2β+1)_ǫ,s kC⁰ ≤ Xβ

γ=0

K_2γ^(2β+1)ku^(2γ)kC^2α+4−2γ + G^(2β+1)kgkC^2α+2

ǫ^2α+3−2β

Multiplying these equations byǫ^2β resp. ǫ^2β+1, adding them together and summing over β ∈ {0, ..., α+ 1} yields, after changing the order of the double summation on

4.3. Construction of consistent population functions 183 by definition. So, we obtain the assertion by choosing

K := max n^2α+3X

Recalling the definition of the constants G⁽ⁿ⁾ and K_2γ⁽ⁿ⁾, we see that K does nei-ther depend on ǫ nor on the u^(2γ)’s and g. In fact, it is only determined by the polynomials P_k (k∈ {0,1, ...,2α+ 3}) and the indicated parameters.

We are now well prepared to prove that ˆf_h solves the discrete lattice-Boltzmann equation on the gridGh approximately, up to a remainder-term of order O(h^2α+4) = O(N_h^−2α−4).

Theorem 4.2. Suppose that the existence and smoothness hypothesis is satisfied.

May f_ǫ ∈ F(S ×[0, T]×R/LR,R) be constructed as described above. Then the subsequent estimate holds true for any admissible mesh-size h ∈ HL and all time steps n∈T_h\ {N_T(h)} with 1≤p≤ ∞.

ˆf_h(n+ 1,·)−B_h(n+ 1)ˆf_h(n,·)

L^p ≤ C_p h^2α+4 (4.39) The constant C_p is independent of h and n

C_p := √^p

where K is the constant from lemma 4.6. Thus f_ǫ is consistent of order 2α+ 4 to the lattice-Boltzmann equation with respect to the supremum-norm and a fortiori with the weaker L^p-norms (1≤p <∞).

For the proof of the theorem, the recursive representation formula (4.12) for the asymptotic order functions will be of great help:

f⁽⁰⁾ = E⁽⁰⁾, f⁽ⁿ⁾ = E⁽ⁿ⁾ − Xn k=1

T_kf^n−k + δ_2nτ gw n≥1. (4.41)

Proof: Let us first consider the case of the supremum norm (p=∞). The asserted inequality is just a short notation of

|D_h(n, i, s)| :=

...

Figure 4.3: Change of indexation in a double sum. See equation (4.43).

that is asserted to hold for any s∈ S,(t_n+h², x_i)∈ Gh. Since we want to prove the estimate independently ofxi ∈ Lh, i.e. in the supremum norm, we have to verify the following inequality, obtained by replacingx_i withx_i+sh, for any (t_n+h², x_i)∈ Gh possible without misinterpretation. Recall in particular T_k=T_k(∂_t, s∂_x).

(∗) :=

Solving (4.37) for the shifted order function, multiplying by hⁿ and summing over n∈ {0,1, ..., α+ 3} yields By summing also the representation formula (4.41) over n we obtain

f_h,s = Es⁽⁰⁾ + kas second index and rename, finally, ν inn), we will change the indexation of the double sum such that it can be put together with the double sum in (4.42).

f_h,s = Eh,s(u_h) − τ

4.3. Construction of consistent population functions 185

Let us now plug (4.42) and (4.43) into (∗). By virtue of lemma (4.5) we have Eh P

˜ s∈Sf_h,˜s

= Eh(u_h). This enables us to put together the equilibrium terms coming directly from the algorithm with those from (4.43). Thus we get:

(∗) = The double sums can be simplified essentially:

2α+3X

Using again (4.43) for the resubstitution (replace third, fourth and fifth term to-gether by −f_h,s) we arrive at Taking into account the estimate of the remainder terms in lemma (4.6), we find the following inequality implying the assertion.

(∗) ≤ K X^α

Initialization and Convergence. Next, let us consider the case α= 0 in detail.

This case is, of course, the simplest as well as the most important one, since it is sufficient to prove convergence of the lattice-Boltzmann scheme. The consideration

of higher orders leads to a truncated asymptotic expansion of the zeroth moment, that additionally describes the leading error terms being not of interest here.

We have already remarked, that u⁽⁰⁾ =u, whereu is just the solution of the target problem (4.1). As no other asymptotic orders are involved, we will always write u instead of u⁽⁰⁾. Let us assume, that the existence and smoothness hypothesis is fulfilled for u and g, such that we have u ∈ C⁴([0, T]× R/LR,R) and g ∈ C²([0, T]×R/LR,R). Then it is possible to construct recursively the following functions by means of the polynomialsT_k(see equation (4.10) and table 4.2), where T_k abbreviates T_k(∂_t, s∂_x).

f⁽⁰⁾:= D(u)

f⁽¹⁾:= A(u) − τ T₁f⁽⁰⁾ f⁽²⁾:= R(u) − τn

T₁f⁽¹⁾ + T₂f⁽⁰⁾o

+ τ gw f⁽³⁾:= −τn

T₁f⁽²⁾ + T₂f⁽¹⁾ + T₃f⁽⁰⁾ o

(4.44)

Note, that this construction procedure is equivalent to (4.35). The approximate population function f_ǫ is then defined by

f_ǫ := f⁽⁰⁾ + ǫf⁽¹⁾ + ǫ²f⁽²⁾ + ǫ³f⁽³⁾ . (4.45) Explicitly we obtain for the asymptotic orders:

f_s⁽⁰⁾(t, x) =w_su(t, x) f_s⁽¹⁾(t, x) =w_sn

θ a s − τ s ∂_xo u(t, x) f_s⁽²⁾(t, x) =−τw_sn

c + θas²∂_x − (τ− ¹₂)s²∂_x² + ∂_to

u(t, x) + τw_sg(t, x) f_s⁽³⁾(t, x) =τ²ws

n

sc∂x+θas³∂_x²−(τ −¹2)s³∂_x³+s∂t∂x

o

u(t, x) + τwss∂xg(t, x)

− τw_s ∂_t+ ¹₂s²∂_x²

θas−τ s∂_x

u(t, x) − τw_sn

s∂_t∂_x+¹₆s³∂_x³o u(t, x)

Considering these functions as polynomials in s∈ S, we observe, thatfs⁽¹⁾ and fs⁽³⁾

have odd parity, while f_s⁽⁰⁾ and f_s⁽²⁾ are of even parity. Hence the zeroth moments of f_s⁽¹⁾ andf_s⁽³⁾ vanish. Since u is the solution of the target equation we find

X

s∈S

f_s⁽²⁾ = −τn

cu+a∂_xu−¹_θ(τ− ¹₂)∂²_xu+∂_tu−go

= 0.

So we obtain, without explicitly referring to lemma 4.5, X

s∈S

fǫ,s = X

s∈S

f_s⁽⁰⁾ + hf_s⁽¹⁾ + h²f_s⁽²⁾ + h³f_s⁽³⁾

= u .

4.3. Construction of consistent population functions 187

Let us now estimate the remainder terms:

|r⁽⁰⁾_h,s(t, x)| ≤ kukC⁴ h⁴

|r⁽¹⁾_h,s(t, x)| ≤ θ|a| kukC³ h³ + τk∂_xukC³ h³ ≤ (θ|a|+τ)

| {z }

=:K1

kukC⁴ h³

|r⁽²⁾_h,s(t, x)| ≤ τ ckukC² +τ|a|θk∂_xukC²+τ(τ −¹₂)k∂_x²ukC² +τk∂_tukC² +τkgkC²

≤ τ c+θτ|a|+τ(τ −¹₂) +τ

| {z }

=:K2

kukC⁴h² + τkgkC²h²

|r⁽²⁾_h,s(t, x)| ≤

τ² c+|a|θ+τ +¹₂

+³₂τ(|a|θ+τ +⁷₆τ

| {z }

=:K3

kukC⁴h+τkgkC²h

As already mentioned, the lattice-Boltzmann scheme is classically initialized by the equilibrium associated with the initial condition for the target problem. The pre-ceding analysis of the lattice-Boltzmann scheme suggests, however, to initialize the population function by the initial value of the approximated population function.

Thus there are four further possibilities of initialization depending on how many orders should be taken into account.

equilibrium: F_h(0, i) = Eh w(xi)

needs:w 0^th order : F_h(0, i) = f_h⁽⁰⁾(0, x_i) ” w 1^storder : F_h(0, i) = P¹

n=0

f_h⁽ⁿ⁾(0, x_i)hⁿ ” w, ∂_xw 2^nd order : F_h(0, i) = P²

n=0

f_h⁽ⁿ⁾(0, x_i)hⁿ ” w, ∂_xw, ∂_x²w, g 3^rdorder : F_h(0, i) = P³

n=0

f_h⁽ⁿ⁾(0, x_i)hⁿ ” w, ∂_xw, ∂_x²w, ∂_x³w, g, ∂_xg

The right column indicates the quantities being necessary for each specific initial-ization. Usually, the initial valuewand the sourcegare provided in the grid nodes at least. For initializations of first and higher order, additionally spatial derivatives of w and g are needed, which must be supplied numerically (e.g. numeric differ-entiation) in the worst case, if w and g are not prescribed analytically. Moreover, the expressions forf⁽²⁾ and f⁽³⁾ contain temporal derivatives of the solution. These can not be computed directly, but they can be procured by taking recourse to the target equation. Solving it for∂_tuyields with u(0,·) =w:

∂_tu(0,·) =ν∂_x²w−a∂_xw−cu+g(0,·)

∂_t∂_xu(0,·) =ν∂_x³w−a∂_x²w−cu+∂_xg(0,·)

Thus, the knowledge of w, g is sufficient to compute f⁽²⁾(0,·) and f⁽³⁾(0,·).

Obviously, the initialization via equilibrium or zeroth order is consistent of order 1; generally the initialization up to order n is consistent of order n+ 1. After

theorem 4.2 our approximated population function is consistent of order 4 to the discrete lattice-Boltzmann equation. In virtue of theorem 4.1, we conclude under the stability assumption, that the inherent deviation is for any fixed, finite time of magnitude h². Thus it is reasonable to choose the initialization so that the inherited initial deviation is of the same size. Any higher order initialization does not considerably decrease the error, whereas a lower order initialization might spoil the result since the inherited initial deviation is then expected to dominate the inherent deviation. Therefore, the initialization should be preferably done up to first order:

F_h(0, i) =w(x_i)w + h θa−τ ∂_x

w(x_i)Sw. (4.46) Theorem 4.3. Consider the target problem (4.1) with w∈C⁴(R/LR,R) and g ∈ C²([0, T]×R/LR,R) and assume u∈ C⁴([0, T]×R/LR,R). If assumption (4.27) is satisfied, then the zeroth moment U_h of the population function converges in the grid nodes to the solution of the target equation:

kU_h(N_t,·)−u(Nˆ _t,·)k^p_L≤3^p−1p E_pJ_p+E_pC_pT h².

Im Dokument Analysis of Lattice-Boltzmann Methods : Asymptotic and Numeric Investigation of a Singularly Perturbed System (Seite 190-200)