Convergence of equilibria for numerical approximations of a suspension model

Contents lists available atScienceDirect

Computers and Mathematics with Applications

journal homepage:www.elsevier.com/locate/camwa

Convergence of equilibria for numerical approximations of a suspension model

J. Valero

^a,∗

, A. Giménez

^a

, O.V. Kapustyan

^b

, P.O. Kasyanov

^c

, J.M. Amigó

^a

aUniversidad Miguel Hernández de Elche, Centro de Investigación Operativa, 03202 Elche (Alicante), Spain

bTaras Shevchenko National University of Kyiv, Kyiv, Ukraine

cInstitute for Applied System Analysis, National Technical University of Ukraine ‘‘Kyiv Polytechnic Institute’’, Kyiv, Ukraine

a r t i c l e i n f o

Article history:

Received 23 December 2015 Received in revised form 11 April 2016 Accepted 23 May 2016

Available online xxxx

Keywords:

Non-Newtonian fluids Suspensions

Numerical approximations Finite-difference schemes Partial differential equations

a b s t r a c t

In this paper we study the numerical approximations of a non-Newtonian model for concentrated suspensions.

First, we prove that the approximative models possess a unique fixed point and study their convergence to a stationary point of the original equation.

Second, we implement an implicit Euler scheme, proving the convergence of these approximations as well.

Finally, numerical simulations are provided.

©2016 Elsevier Ltd. All rights reserved.

1. Introduction

Non-Newtonian (or complex) fluids often appear in nature and industry. Good examples of such fluids are toothpaste, ketchup, magma, blood, mucus or emulsions such as mayonnaise among many others. A special type of complex fluids are concentrated suspensions, which can be found, for example, in medicine (blood) or in building industry (cement). The dynamical behaviour of suspensions is still far from being well understood as developing a faithful mathematical model of such processes is not an easy task.

We are interested in an equation modelling suspensions which was proposed in [1]. In the last years, several authors have studied for this equation the existence and uniqueness of solutions [2,3], the asymptotic behaviour [4,5] and numerical approximations [6–8].

In our previous paper [8] we studied a sequence of approximative problems for this model, in which finite-difference schemes were used to deal with the partial derivative with respect to the spatial variable. The problem was split in three steps: a partial differential equation with a large diffusion, an infinite system of ordinary differential equations and finally a finite system of ordinary differential equations. For initial data satisfying suitable assumptions it was proved that the iterate limit of the solutions of the approximative problems in the spaceC

( [

0

,

^T

] ,

^L2

(

R

))

is equal to the solution of the original equation.

In this paper we extend the results from [8] in two ways.

∗Corresponding author.

E-mail addresses:[email protected](J. Valero),[email protected](A. Giménez),[email protected](O.V. Kapustyan),[email protected] (P.O. Kasyanov),[email protected](J.M. Amigó).

http://dx.doi.org/10.1016/j.camwa.2016.05.034 0898-1221/©2016 Elsevier Ltd. All rights reserved.

First, we study the convergence of the fixed points of the approximative problems. It is well-known [2] that for certain values of the parameters of the equation there exists a unique fixed point of the problem with support not included in the interval

[−

1

,

¹

]

. This equilibrium is asymptotically stable [4] and the numerical simulations in [7] suggest that every solution with initial data with support not included in

[−

1

,

¹

]

converges to this fixed point as time goes to

+∞

. We prove that each of the approximative problems possesses a unique fixed point and also that the iterate limit of the equilibria of the approximative problems in the spaceL²

(

R

)

is equal to the equilibrium of the original equation with support not included in the interval

[−

1

,

¹

]

.

Second, we complete the sequence of approximations of the problem by implementing an implicit Euler scheme for the discretization of the time derivative. We prove that the solution of the resulting system converges in the space C

( [

0

,

^T

] ,

^L2

(

R

))

to the solution of the finite system of ordinary differential equations approximating the original equation.

Finally, some numerical simulations are provided in the last section.

2. Previous results

In the previous paper [8] the authors considered the convergence of finite-difference approximations of the problem

∂

^p

∂

^t

−

D

(

^p

(

^t

)) ∂

^2p

∂σ

²

+

¹

T₀

χ

R\[−1,1]

(σ)

^p

=

^D

(

^p

(

^t

))

α δ

0

(σ) ,

⁽¹⁾

p

≥

0

,

^p

(

⁰

, σ) =

p⁰

(σ) ,

⁽²⁾

wherep

=

p

(

^t

, σ) ,

^t

∈ [

0

,

^T

] , σ ∈

_R,T₀and

α

are positive constants.

Here,

δ

0is the Dirac

δ

-function with support in the origin, D

(

^p

(

^t

)) = α

T₀



|σ|>1

p

(

^t

, σ)

^d

σ

and

χ

Iis the indicator function in the intervalI.

The functionp

(

^t

, σ )

is a probability density at timet, so for anyt

∈ [

0

,

^T

]

,



R

p

(

^t

, σ )

^d

σ =

1

,

⁽³⁾

p

(

^t

, σ ) ≥

0

,

^{for a.a.}

σ ∈

_R

.

It is well-known [2] that for anyp⁰

∈

L¹

(

R

) ∩

L^∞

(

R

)

^{such thatp0}

≥

0 a.e.,



Rp⁰

(σ)

^d

σ =

1,



R

| σ |

p⁰

(σ)

^d

σ < ∞

and D

(

^p0

) >

0 there exists a unique solutionp

=

p

(

^t

, σ)

^{of problem}(1)–(2), which satisfies(3).

We consider as a first step the approximative problem

∂

tp^c

−



D



p^c

(

^t

)  +

¹

c



∂

_{σ σ}^{2 p}

+

¹

T₀

χ

R\[−1,1]

(σ )

^pc

=

^D

(

^pc

(

^t

))

α δ

c

(σ ) ,

⁽⁴⁾

p^c

≥

0

,

^pc

(

⁰

, σ) =

p⁰_c

(σ ) ,

⁽⁵⁾

wherep^c

=

p^c

(

^t

, σ )

^,c

>

0 is a large parameter and the

δ

^-function

δ

0is replaced by the step continuous from the right function

δ

c

(σ) =



 

 

 

 

 

 

0

,

^if

σ < −

¹ 2c

,

c

,

^if

−

¹

2c

≤ σ <

¹ 2c

,

0

,

^if

σ ≥

¹

2c

.

We would like to highlight the fact that the new term¹_c

∂

_{σ σ}^{2 p}is an artificial diffusion which helps us to prove the convergence of the approximative solutions. Such a trick is very common in the numerical approximations of problems in Physics. Also,

[−

¹

2c

,

_2c¹

]

is the support of the map

δ

c, which approximates the

δ

^-function

δ

0. Therefore, whenc

→ +∞

, the artificial diffusion and the support of

δ

cconverge to 0 in unison.

Letp⁰_cbe such that

p⁰_c

∈

C₀^∞

(

R

),

^p0c

≥

0 a.e.,



R

p⁰_c

(σ )

^d

σ =

1

,

⁽⁶⁾

p⁰_c

→

p⁰ inL²

(

R

), σ

^p0c

→ σ

^{p0 inL1}

(

R

),

^asc

→ +∞ .

⁽⁷⁾

It is proved in [8, Theorem 3] that such approximation exists and that the unique solutionp^cto problem(4)–(5)converges to the unique solutionpto problem(1)–(2)in the spaceC

( [

0

,

^T

] ,

^X

)

^{, where}

X

=



p

∈

L²

(

R

) :



R

| σ | |

p

|

d

σ < +∞

 ,

endowed with the norm

∥

p

∥

_X

= ∥

p

∥

_L2(R)

+ 

R

| σ | |

p

|

d

σ

. It is important to remark that for allt

∈ [

0

,

^T

]

the solution satisfies



Rp^c

(

^t

, σ)

^d

σ =

1 andp^c

(

^t

, σ) ≥

0, for a.a.

σ ∈

_R, since it is a probability density.

Further, we shall consider the following approximating infinite system of ordinary differential equations:

dp^c_i^,h dt

−



D_h

(

^pch

(

^t

)) +

¹ c



_pc,h

i+1

(

^t

) −

2p^c_i^,h

(

^t

) +

p^c_i−^,h₁

(

^t

)

h²

+

¹

T₀

χ

Z\[−2n1,2n1]

(

ⁱ

)

^pci^,h

(

^t

) =

^Dh



p^c_h

(

^t

) 

α δ

cⁱ

,

⁽⁸⁾

p^c_i^,h

(

⁰

) =

p⁰_c,h,i

,

ⁱ

∈

_Z

,

⁽⁹⁾

whereh

>

⁰

σ

i

=

ih,i

∈

_Z,¹_h

=

2n₁, withn₁

∈

_N,

χ

Z\[−2n1,2n1]

(

ⁱ

) =



1

,

^ifi

̸∈

[

−

2n₁

,

²ⁿ1]

,

0

,

^otherwise

.

Here,p^c_h

(

^t

) =



p^c_i^,h

(

^t

) 

i∈_Zdenotes a sequence satisfying p^c_i^,h

(

^t

) ≃

p^c

(

^t

, σ

i

) ,

and we made the following approximations D



p^c

(

^t

) 

= α

T₀



|σ|>1

p^c

(

^t

, σ)

^d

σ ≃

D_h



p^c_h

(

^t

) 

= α

^h T₀



|_i|>2n1

p^c_i^,h

(

^t

),

⁽¹⁰⁾

δ

ⁱc

= δ

c

(σ

i

) =



₀

,

^ifi

< −

n_h,c

,

c

,

^if

−

n_h,c

≤

i

<

ⁿh,c

,

0

,

^ifi

≥

n_h,c

.

⁽¹¹⁾

Also,cis taken such that_2ch¹

=

n_h,c

∈

_Nandn_h,c

<

²ⁿ1.

We observe that the parameterhis the length of the intervals in the finite-difference approximation of the second derivative

∂

_{σ σ}² p, which is a diffusion term. The approximation is getting better ashgoes to 0. Also,n_h,c is the number of subintervals of lengthhof the interval

[

0

,

_2c¹

]

, the support of the function

δ

cin the positive semi-axis. The condition n_h,c

<

²ⁿ1is equivalent to say that_2c¹

<

1, that is, the support of

δ

cis strictly included in the interval

[−

1

,

¹

]

.

We define the following partition of the real line:

Ωh

= { σ

i

=

ih

}

_i∈Z

,

^Ii^h

= [ σ

i

, σ

i

+

h

).

Forp⁰_cfrom(6)we define the step function



^p

0

c,^h

(σ) = 

i∈_Z

p⁰_c

(

^ih

) χ

I_i^h

(σ )

and normalize it by setting

p⁰_c,h

(σ) = 

^p

0 c,^h

(σ )





^p0_c,h





L¹(R)

.

⁽¹²⁾

It holds thatp⁰_c,h

→

p⁰_cinL²

(

R

) ∩

L¹

(

R

)

^ash

→

0.

Further, we fixc

∈

_Nand take a sequenceh_n

→

0 such that _2ch¹

n

=

n_hn,c

∈

_N. Then conditions _h¹

n

=

2n₁

,

ⁿ1

∈

_N, n_hn,c

<

²ⁿ1are satisfied. We takep⁰_c,h

n,i

=

p⁰_c,h

n

(

^ihn

)

as the initial data in problem(8)and define the step functions p^c_hn

(

^t

, σ ) = 

i∈Z

p^c_i^,hn

(

^t

) χ

I_i^h

(σ ),

⁽¹³⁾

wherep^c_hn

(

^t

) = {

p^c_i^,hn

(

^t

) }

_i∈Zis the unique solution to problem(8)–(9). It is proved in [8, Theorem 2] that

p^c_hn

→

p^c strongly inC

( [

0

,

^T

];

L²

(

R

)),

⁽¹⁴⁾

wherep^cis the solution to problem(4)–(5)with initial datap⁰_c. As before, the solutionsp^c_hnsatisfy that



R

p^c_hn

(

^t

, σ)

^d

σ = 

i∈_Z

p^c_i^,hn

(

^t

)

^hn

=

1 and p^c_i^,hn

(

^t

) ≥

0

,

^{for anyt}

∈ [

0

,

^T

] ,

ⁱ

∈

_Z

.

Let us consider now finite-dimensional approximations. We define the operatorA^N_h

:

_R^2N+1

→

_R^2N+1by

A^N_h

:=

¹ h²

























1

−

1 0

· · ·

0 0 0

−

1 2

−

1 0

· · ·

0 0 0

−

1 2

−

1

... · · ·

0

...

⁰

... ... ...

⁰

...

0

· · · ... −

1 2

−

1 0 0 0

· · ·

0

−

1 2

−

1 0 0 0

· · ·

0

−

1 1

























(2N+1)×(2N+1)

.

⁽¹⁵⁾

Then we consider the finite-dimensional system



 



 



dp^c_h^,_,^N_i

dt

= −



D^N_h



p^c_h^,N

(

^t

) 

+

¹ c

 

A^N_hp^c_h^,N



i

−

¹

T₀

χ

Z\[−2n1,2n1]

(

ⁱ

)

^pch^,,^Ni

+

D^N_h



p^c_h^,N

(

^t

) 

α δ

ⁱc

,

p^c_h^,_,^N_i

(

⁰

) =

p^N_c,^,_h⁰_,i

, −

N

≤

i

≤

N

,

(16)

whereN

>

²ⁿ1,¹_h

=

2n₁

,

ⁿ1

∈

_N,_2ch¹

=

n_h,c

∈

_N

,

ⁿh,c

<

²ⁿ1and D^N_h



p^c_h^,N



= α

^h T₀



2n1<|i|≤N

p^c_h^,_,^N_i

.

What we have done is to cut the tails of the system(8)off in order to work with a finite number of equations. The condition N

>

²ⁿ1implies that we solve the problem for

σ

in an interval containing

[−

1

,

¹

]

. It is obvious thatNhas to be large to get good approximations.

For the initial datap⁰_c,hfrom(12)we consider the approximationsp^N_c,^,_h⁰given by

p^N_c,^,_h⁰_,i

=

^p

0 c,h,i



|_i|≤_N

hp⁰_c,h,i

,

^for

|

i

| ≤

N

,

⁽¹⁷⁾

wherep⁰_c,h,i

=

p⁰_c,h

(

^ih

)

. Then we define the step functions p^c_h^,N

(

^t

, σ ) = 

i∈Z

p^c_h^,_,^N_i

(

^t

) χ

Ii

(σ ),

⁽¹⁸⁾

where



p^c_h^,_,^N_i

( · ) 

|_i|≤_N is the unique solution to problem(16)with initial data(17)andp^c_h^,_,^N_i

(

^t

) =

0 if

|

i

| >

N. It is proved in [8, Section 5] that

p^c_h^,N

→

p^c_h inC

( [

0

,

^T

] ,

^L2

(

R

))

^asN

→ ∞ ,

wherep^c_his the function defined in(13)with initial data(12). Again, the property of being a probability density is satisfied:



R

p^c_h^,N

(

^t

, σ)

^d

σ = 

|i|≤N

p^c_h^,_,^N_i

(

^t

)

^h

=

1 and p^c_h^,_,^N_i

(

^t

) ≥

0

,

^{for anyt}

≥

0

, |

i

| ≤

N

.

3. Fixed points of approximations

Our aim in this section is to study the fixed points of the approximative problems and their convergence to the fixed points of the original problem(1). For simplicity, we shall consider the particular case whereT₀

=

1.

First, recall that for Eq.(1)withT₀

=

1 the fixed points, given by the solutions of D

(

^p

) ∂

^2p

∂σ

²

− χ

R\[−1,1]

(σ)

^p

+

^D

(

^p

)

α δ

0

(σ) =

0

,

⁽¹⁹⁾

are the following [2]:

•

Any probability densityp

( · )

with support in

[−

1

,

¹

]

solves(19). We note that all these solutions satisfyD

(

^p

) =

0;

•

If

α ≤

¹

2, there are no more equilibria. If

α >

¹₂, then there exists a unique fixed pointpwith positive value ofD

(

^p

)

^, which is given by

p

(σ ) =



 

 

 

 

 



 

 

 

 

 



√

D^∗ 2

α

^e(1

+σ )/√ D∗

,

^if

σ ≤ −

1

,

√

D^∗

+

1

2

α +

¹

2

α σ,

^if

−

1

≤ σ ≤

0

,

√

D^∗

+

1

2

α −

¹

2

α σ,

^{if 0}

≤ σ ≤

1

,

√

D^∗ 2

α

^e(1

−σ )/√ D∗

,

^if

σ ≥

1

,

(20)

where

D^∗

=



−

¹ 2

+

√

4

α −

1

2



2

(21)

andz

=

√

D^∗is the unique positive solution of the equation h

(

^z

) =

z²

+

z

− α +

¹

2

=

0

.

We observe that when

α >

¹₂the stationary pointpis asymptotically stable [4]. Moreover, the numerical simulations in [7] suggest that every solution with initial data satisfyingD



p⁰



>

0 converges to this fixed point as time goes to

+∞

_. We shall prove that for

α >

¹₂the approximative problems possess a unique fixed point converging to(20).

3.1. Equation with large diffusion

Let us consider now the fixed points of problem(4). In order to find them we fix firstD

>

0 and solve first the following ordinary differential equation:



D

+

¹

c



d²p^c

d

σ

²

− χ

R\[−1,^1]

(σ)

^pc

+

^D

α δ

c

(σ) =

0

.

We note that in this case, unlike problem(1), there is no stationary solutions withD

(

^p

) =

0.

Taking into account the conditionp^c

(σ) →

0, as

σ → ±∞

, it is not difficult to check that this equation possesses a unique solution defined by

p^c_D

(σ ) =



 

 

 

 

 

 

 

 

 

 

 



 

 

 

 

 

 

 

 

 

 

 



D 2

α 

D

+

¹

c

e^{(1+σ )/}

 D+¹_c

,

^if

σ ≤ −

1

,

D 2

α

1

+



D

+

¹

c

D

+

¹

c

+

^D

2

α 

D

+

¹

c

 σ,

^if

−

1

≤ σ ≤ −

¹ 2c

,

D

2

α

1

+



D

+

¹

c

D

+

¹

c

−

¹ 8

D

α 

D

+

¹

c



c

−

^Dc 2

α 

D

+

¹

c

σ

²

,

^if

−

¹

2c

≤ σ ≤

¹ 2c

,

D

2

α

1

+



D

+

¹

c

D

+

¹

c

−

^D

2

α 

D

+

¹

c

 σ,

^{if 1}

2c

≤ σ ≤

1

,

D

2

α 

D

+

¹

c

e^{(1−σ )/}

 D+¹_c

,

^if

σ ≥

1

.

Since D



p^c_D



= α



|σ|>1

p^c_D

(σ)

^d

σ =

D

,

^{for anyD}

>

⁰

,

in order to obtain a fixed point it remains to find a positive value ofDsuch that



Rp^c_D

(σ )

^d

σ =

1. Calculating the integral we obtain

1 48c²

α

2D D

+

¹

c



24c

+

24c²D

+

24c²



D

+

¹

C

+

12c²

−

1



=

1

,

and after the change of variablez

=



D

+

¹

c we finally have the equation g^c

(

^z

) =

z⁴

+

z³

−



1 24c²

+

¹

c

−

¹ 2

+ α



z²

−

¹

cz

−

¹ 2c

+

¹

24c³

=

0

.

⁽²²⁾

It follows from the Descarte’s rule of signs that for

α >

⁰

.

^{5 andc}large enough this polynomial possesses a unique positive rootz^c. More precisely,chas to satisfyc

>

^√1₁₂. We will takec

≥

1. Such condition is compatible with the meaning of the term¹_c

∂

_{σ σ}^{2 pin}(4), as this is an artificial diffusion that has to be small in order to approximate the original system properly, which means that we needcto be large.

If we pass to the limit asc

→ ∞

the polynomialg^c

(

^z

)

^{tends to} g

(

^z

) =

z²



z²

+

z

+

¹ 2

− α

 .

By continuity, the rootz^cconverges to the unique positive root ofh

(

^z

)

, which is equal toz^∗

=

√

D^∗. Therefore, D^c

= 

z^c



2

−

¹

c

→

D^∗

>

⁰

,

whereD^∗is given in(21). Hence,D^c

>

^{0, forc}large enough, and thus there is a unique stationary pointp^c

(σ) =

p^c_Dc

(σ )

^. Moreover, it is easy to see usingD^c

→

D^∗that

p^c

→

p inX

.

Therefore, we have proved the following result.

Theorem 1. Let

α >

⁰

.

^{5. Problem(4)}possesses a unique fixed point p^cfor c

≥

1and p^c

→

p in X

,

^{as c}

→ ∞ ,

where p is the unique fixed point of problem(1)such that D

(

^p

) >

^{0defined in(20).}

3.2. Lattice dynamical system

Further, we will study the fixed points of Eq.(8)with_2ch¹

=

n_h,c

∈

_Nandn_h,c

<

²ⁿ1

=

¹

h. As before, we fix firstD

>

⁰ and solve the following equation in differences



 

 

 

 

 



 

 

 

 

 



−

^D

+

¹

c

h²

(

^pi+1

−

2p_i

+

p_i−1

) +

p_i

=

0

,

^ifi

< −

2n₁

,

p_i+1

−

2p_i

+

p_i−1

=

0

,

^if

−

2n₁

≤

i

< −

n_h,c

,

−

^D

+

¹

c

h²

(

^pi+1

−

2p_i

+

p_i−1

) =

^D

α

^c

,

^if

−

n_h,c

≤

i

<

ⁿh,c

,

p_i+1

−

2p_i

+

p_i−1

=

0

,

^ifnh,^c

≤

i

≤

2n₁

,

−

^D

+

¹

c

h²

(

^pi+1

−

2p_i

+

p_i−1

) +

p_i

=

0

,

^ifi

>

²ⁿ1

,

(23)

whose solution, taking into account thatp_i

→

i→±∞0, is given by

p^c_i,^,_D^h

=



 

 

 



 

 

 



C₁

λ

ⁱ1

,

^ifi

< −

2n₁

,

A

+

Bi

,

^if

−

2n₁

≤

i

< −

n_h,c

,

E

+

Fi

−

^ch

2D 2



D

+

¹

c

 α

ⁱ

2

,

^if

−

n_h,c

≤

i

<

ⁿh,c

,

A

+

Bi

,

^ifnh,c

≤

i

≤

2n₁

,

C₂

λ

⁻1ⁱ

,

^ifi

>

²ⁿ1

,

(24)

provided that

λ

1

= λ

^c1^,,^Dh

=

1

+

^h

2

2



D

+

¹

c

 +

^h



D

+

¹

c



1

+

^h

2

4



D

+

¹

c

, λ

2

= λ

^c2^,,^Dh

=

¹

λ

^c1^,,^Dh

,

⁽²⁵⁾

and the constantsC₁

,

^A

,

^B

,

^E

,

^F

,

^A

,

^B

,

^C2satisfy the compatibility conditions



 

 

 

 

 

 

 

 

 

 

 

 

 

 



 

 

 

 

 

 

 

 

 

 

 

 

 

 



A

−

¹

hB

− λ

⁻1^1hC₁

=

0

,

A

−

^h

+

1

h B

− λ

⁻1^h+h1C₁

=

0

,

A

−

¹

2chB

−

E

+

¹

2chF

= −

^D

8c



D

+

¹

c

 α ,

A

−

^2ch

+

1

2ch B

−

E

+

^2ch

+

1

2ch F

= −

^D

(

¹

+

2ch

)

² 8c



D

+

¹

c

 α ,

A

+

¹

−

2ch

2ch B

−

E

−

¹

−

2ch

2ch F

= −

^D

(

¹

−

2ch

)

² 8c



D

+

¹

c

 α ,

A

+

¹

2chB

−

E

−

¹

2chF

= −

^D

8c



D

+

¹

c

 α ,

A

+

^h

+

1

h B

− λ

⁻1^h+h1C₂

=

0

,

A

+

¹

hB

− λ

⁻1^1hC₂

=

0

.

(26)

Solving this system we obtain:

C₁

= λ

1^1hD



D

+

¹

c

 α

 −

1

+ λ

⁻1¹



h

(

²

+

h

) −

2h² 4



−

1

+ λ

⁻1¹

 

1

+

h

− λ

⁻1¹

 ,

^C2

= λ

1^1hD



D

+

¹

c

 α

 −

1

+ λ

⁻1¹



h

(

²

−

h

) −

2h² 4



−

1

+ λ

⁻1¹

 

1

+

h

− λ

⁻1¹

,

⁽²⁷⁾

A

= 

^D D

+

¹

c

 α

 −

1

+ λ

⁻1¹

 (

²

+

h

) −

2h 4



−

1

+ λ

⁻1¹

 ,

^A

= 

^D D

+

¹

c

 α

 −

1

+ λ

⁻1¹

 (

²

−

h

) −

2h 4



−

1

+ λ

⁻1¹

 ,

B

= −

D



D

+

¹

c

 α

 −

1

+ λ

⁻1¹



h

(

²

+

h

) −

2h² 4



1

+

h

− λ

⁻1¹

 ,

^B

= 

^D D

+

¹

c

 α

 −

1

+ λ

⁻1¹



h

(

²

−

h

) −

2h² 4



1

+

h

− λ

⁻1¹

 ,

F

= −

^Dh 2



D

+

¹

c

 α (

¹

+

ch

) −

^D



D

+

¹

c

 α

 −

1

+ λ

⁻1¹



h

(

²

+

h

) −

2h² 4



1

+

h

− λ

⁻1¹

 ,

E

= −

^D

8c



D

+

¹

c

 α (

¹

+

2ch

) +

^D



D

+

¹

c

 α

 −

1

+ λ

⁻1¹

 (

²

+

h

) −

2h 4



−

1

+ λ

⁻1¹

 .

For simplicity of notation here and throughout the paper, if no confusion is possible, sometimes we omit the indexes c

,

^h

,

^Dand write just

λ

1.

We need to check first that

α

^h



|i|>2n1p^c_i,^,_D^h

=

D. Indeed, we can easily compute that

α

^h



|i|>2n1

p^c_i,^,_D^h

=

^Dh

2



D

+

¹

c

 λ

1

(λ

1

−

1

)

²

=

D 4



D

+

¹

c

 +

2h²

+

2h



4



D

+

¹

c

 +

h²



h

+



4



D

+

¹

c

 +

h²



2

=

D for anyD

>

⁰

.

We need to findD^c_h

>

0 such thatS_h^c

= 

i∈_Zp^c_i,^,_D^hc h

h

=

1. Using mathematical software we obtain

S^c_h

(

^D

) = 

i∈_Z

p^c_i,^,_D^hh

=

^b

c h

(

^D

)

w

h^c

(

^D

) ,

⁽²⁸⁾