Adapted Lie splitting method for convection-diffusion problems with singular convective term
Thi Tam Dang Trung Hau Hoang Giandomenico Orlandi


  



theoremTheorem[section]
lemmaLemma[section]
propositionProposition[section]
assumptionAssumption[section]

examExample[section]


example[1][Example]
		
  #1
remark[1][Remark]
		
  #1
corollary[theorem]Corollary

definition
definitionDefinition[section]




equationsection
arabic


plain
























	
1]Thi Tam DangThi-Tam.Dang@uibk.ac.at
1]Trung Hau HoangTrung-Hau.Hoang@uibk.ac.at
2]Giandomenico Orlandigiandomenico.orlandi@univr.it
[1]Department of Mathematics, University of Innsbruck, Austria
[2]Dipartimento di Informatica, Università di Verona, Italy



	Adapted Lie splitting method for convection-diffusion problems with singular convective term
    [
    
============================================================================================




Splitting methods are a widely used numerical scheme for solving convection-diffusion problems. However, they may lose stability in some situations, particularly when applied to convection-diffusion problems in the presence of an unbounded convective term. In this paper, we propose a new splitting method, called the “Adapted Lie splitting method”, which successfully overcomes the observed instability in certain cases. Assuming that the unbounded coefficient belongs to a suitable Lorentz space, we show that the adapted Lie splitting converges to first-order under the analytic semigroup framework. Furthermore, we provide numerical experiments to illustrate our newly proposed splitting approach.




§ INTRODUCTION

Let Ω be an open bounded subset of ℝ^N with N≥ 2. We consider convection-diffusion problems of the form

    {   ∂_tu(t,x) -D u(t,x) - ∇· ( c(x) u(t,x))  = f(t,x),   (t,x) ∈ [0,T]×Ω,
        u(t,x)|_∂Ω = b(t,x),
        u(0,x) = u_0(x),	
    	.

where D denotes a second-order strongly elliptic differential operator, the convective coefficient c(x) is possibly unbounded and belongs to the weak-Marcinkiewicz space L^N, ∞(Ω). We impose a possibly time-dependent Dirichlet boundary condition b: [0,T] ×∂Ω→ℝ, which is assumed to be continuously differentiable. In the context of convection-diffusion problems, the periodic boundary conditions are of significance within theoretical analyses, whereas the imposition of Dirichlet boundary conditions provides a practical approach to enforcing desired behaviors or constraints on the system, consistent with the realities of a wide range of physical and engineering problems.



In the context of linear homogeneous systems, equation (<ref>) is commonly referred to as the Fokker–Planck equation, which characterizes the evolution of certain Brownian motion processes (see <cit.>). When the convective coefficient c(x) remains bounded, the investigation of problem (<ref>) has garnered attention from numerous researchers. However, in various applications, particularly within the realm of semiconductor device diffusion models, the boundedness of the convective field may not be assured. The theoretical analysis of (<ref>) and its generalisations has led to the introduction of several notions of solution, notably including the concept of renormalized solution (see, e.g., <cit.>) and entropy solution (see <cit.>), that are designed to address scenarios where the convective coefficient lacks boundedness, and thus provide a particularly well-suited framework for analyzing the problem under consideration also from a numerical point of view. 





In the field of computational mathematics and scientific computing it is crucial to develop effective methods for solving convection-diffusion problems like (<ref>). Among the various methods being studied, splitting methods have received significant attention as a particularly promising approach (see, e.g., <cit.>). These methods decompose the problem into more manageable subproblems and offer a way to improve computational efficiency and effectiveness. Furthermore, splitting methods have an advantage in preserving some relevant features of the solution such as positivity, as long as the solutions of the subproblems also enjoy that feature (see <cit.>). This makes splitting methods appealing in practical applications where for instance maintaining positivity is crucial, such as modeling physical quantities like concentrations or densities. A standard splitting approach involves splitting the convection-diffusion problem (<ref>) into the following two subproblems:

    {   ∂_t v(t,x) = D v(t,x) + f(t,x), 
       	v(t,x)|_∂Ω= b(t,x),
    .
	
and

    {   ∂_t w(t,x)=  ∇· (c(x) w(t,x)),
       	w(t,x)|_Γ = b_1(t,x),
    	.
	
where Γ denotes the  inflow boundary and b(t,x)=b_1(t,x) for x ∈Γ.	

The diffusion problem (<ref>) can be efficiently treated by a fast Poisson solver (see <cit.>) (with complexity 𝒪(n log n)) where n is the number of degrees of freedom. The pure convection problem (<ref>) can be solved by appropriate methods. Our study focuses on the Lie splitting method, which is a well-known first-order method (see, e.g., <cit.>), for solving the convection-diffusion problem (<ref>).  Given an initial data u_n, fixing a time step size τ>0, one step of the classical Lie splitting iteration in the time interval [ t_n, t_n+1] is given by

    u_n+1 = ℒ_τu_n=  φ_τ^D,f∘φ_τ^C(u_n),

where φ_t^D,f and φ_t^C are the solution operators of problems (<ref>) and (<ref>) respectively. The operator ℒ_τ defined in (<ref>) is referred to as Lie splitting with time step size τ.



Splitting methods are shown to be effective when applied to convection-diffusion problems with a bounded convective term, i.e. when the coefficients are assumed to be sufficiently smooth (see <cit.>). However, in scenarios featuring unbounded convective coefficients many additional challenges are introduced, both for the theoretical analysis and from the numerical perspective. Our numerical experiments (see Section 5 below) have shown a notable instability of splitting methods when confronted with the presence of unbounded convective terms, due to singularity effects. This observation has prompted us to investigate potential modifications of splitting methods in order to minimize the instability encountered. Efforts to develop reliable numerical techniques that have the ability to handle the complexities of unbounded convective terms are motivated by the need to improve the accuracy and dependability of computational simulations in this field.



In this paper we present the “Adapted Lie splitting method”, an innovative approach developed to tackle the instability issue commonly found in the standard Lie splitting approach (<ref>). The proposed method aims to reduces the impact of singularity effects by decomposing the unbounded convective coefficient c(x) into a combination of a bounded function and a L^N,∞(Ω) function for which there is a control on its norm. Our newly approach thus centers around reformulating problem (<ref>) by strategically dividing the convection term into a controlled component and another that exhibits bounded convective behavior. This reformulation, which is detailed in Section 3, effectively avoids the inherent instability conditions encountered. Furthermore, our novel splitting approach offers a promising technique for achieving numerical stability and greatly enhancing accuracy in broader scenarios characterized by unbounded convective terms. 




The convergence results for the newly proposed "Adapted Lie splitting" method are obtained in our main result, Theorem <ref>, under suitable assumptions on the convective term (see Assumption <ref>). These assumptions ensure stability properties and uniform boundedness of the solution. Additionally, we conduct a thorough convergence analysis to demonstrate that the adapted Lie splitting method is a first-order convergent method.


The paper is structured as follows: Section 2 provides a review of the definition and properties of Lorentz spaces. In Section 3, we introduce an adapted Lie splitting method that effectively overcomes instability and significantly improves the accuracy of the numerical method.  In Section 4  our main theoretical convergence analysis result is presented in Theorem <ref>, which establishes the first-order convergence of the adapted Lie splitting under suitable assumptions on the convective coefficients, in an analytic semigroup framework. A rigorous error analysis is carried out to explain the behavior observed in the numerical simulations. To demonstrate the effectiveness of the proposed splitting method, numerical experiments have been conducted in both two- and three-dimensional spaces, as outlined in Sections 5 and 6. 






§ PRELIMINARY DEFINITIONS AND RESULTS

Let Ω be a bounded domain of ℝ^N. For any 1<p,q<∞ and h>0, we define the distribution function of f given by

    λ_g(h) = |{ x ∈Ω: | g(x) | >h }|.

The Lorentz space L^p,q(Ω) consists of all measurable functions g which satisfies the following quantity

    g _p,q^q := p ∫_0^∞[ λ_g(h)]^q/p h^q-1dh ≤∞.

The Lorentz space L^p,q(Ω) becomes a Banach space when endowed with a norm ·_p,q. In the case where p=q, the Lorentz space L^p,p(Ω) is idenitified as the Lebesgue space L^p(Ω). For q = ∞, the class L^p, ∞(Ω) coincides with the Marcinkiewicz class weak-L^p(Ω), and thus

    g_p, ∞^p:= sup_h>0 h^pλ_g(h) < ∞.

A typical example of an element of L^N, ∞(Ω) is the function g(x) = M / | x |, where M>0 is a constant. Further background on the Lorentz spaces can be found in  <cit.>.

For 1 ≤ q <p<r ≤∞, in the Lorentz spaces we have the following continuous inclusions

    L^r(Ω) ⊂ L^p,q(Ω) ⊂ L^p,r(Ω) ⊂ L^p, ∞(Ω) ⊂ L^q(Ω).

Let p, q ∈ [1, ∞] with 1/p+ 1/p^'= 1, 1/q + 1/q^'=1. For f ∈ L^p,q(Ω) and g ∈ L^p^',q^'(Ω) we have the Hölder-type inequality 

    fg_1 ≤ f_p,qg _p^', q^'.


Functions in L^∞(Ω) are not dense in L^p,∞(Ω). For any g ∈ L^N, ∞(Ω) we define

    dist_L^N, ∞(Ω)(g, L^∞(Ω))= inf_ϕ∈ L^∞(Ω) g-ϕ_L^N, ∞(Ω).

We have the  characterization 

    dist_L^N, ∞(Ω)(g, L^∞(Ω)) = lim_ K →∞ g -T_Kg_N, ∞,

where T_K denotes the truncation operator at level ± K, i.e. T_K(s)= s/| s |min{| s |, K }, for all s ∈ℝ. For example, the distance between g(x) = M /| x | in L^N, ∞(Ω) and L^∞(Ω) is given by the following formula:

    dist_L^N, ∞(Ω)(g, L^∞(Ω)) = M ω_N^1/N

where ω_N denotes the Lebesgue measure of the unit ball of ℝ^N. For further details regarding the distance to L^∞(Ω), we refer the reader to <cit.>.

For any given δ≥ 0, we define the subset X_δ(Ω) of L^N, ∞(Ω) consists of all functions g ∈ L^N, ∞ such that

    g -T_Kg_L^N, ∞(Ω)≤δ.

The existence  and uniqueness of a solution to the problem (<ref>) is given in the following theorem
(Existence and uniqueness,<cit.>)
	There exists a proper convex set X_δ(Ω) of the space L^N, ∞(Ω) such that, if 
	
    c ∈ X_δ(Ω)

	then the problem (<ref>) admits a unique solution u ∈ L^2(0,T; W^1,2(Ω)) ∩ C([0,T]; L^2(Ω)).

This section concludes with the Sobolev embedding theorem in Lorentz spaces, which is employed in the proof of Lemma <ref> below.
(Sobolev embedding <cit.>) 
	Assume that 1<p<N, q ≤ q ≤ p, then every function u ∈ W_0^1,1(Ω) verifying |∇ u |∈ L^p,q(Ω) actually belongs to L^p^*,q(Ω), where p^*:= Np/N-p is the Sobolev exponent of p and 
	
    u _p^*,q≤ S_N,p∇ u _p,q

where S_N,p is the Sobolev constant given by S_N,p= ω_N^-1/Np/N-p.









§ DESCRIPTION OF THE METHOD

In this section we propose a new adaptive splitting approach for problem (<ref>) that avoids the instability caused by the unbounded covective term. Let us provide first an auxiliary model problem.


 §.§ Model problem

For any fixed K >0 and  c(x) ∈ L^N, ∞(Ω), define

    θ(x) = T_K(c(x))/c(x),

where T_K is the usual truncation operator at level ± K. Consider the following auxiliary parabolic problem:

    {   ∂_tu(t) - Du(t) -  ∇·( (1- θ(·)) u(t) )  = f(t)+ ∇· ( θ(·)  u(t)),
        u(t)|_∂Ω = b(t),
        u(0)= u_0,
    	.

where D denotes a second-order strongly elliptic differential operator and θ(·) denotes the space-dependent adaptive coefficient. Henceforth, we adopt the notation u(t) instead of u(t,x) and θ(·) instead of θ(x) for the sake of notational brevity.

The second-order linear elliptic differential operator D is defined as follows:

    Du = ∑_i,j=1^N- ∂_i(a_ij(x) ∂_j u)  + ∑_i=1^Nα_i(x)∂_iu+ β(x) u,

where the matrix-value function (a_i,j(x)) ∈ℝ^N × N is assumed to be symmetric, i.e. a_i,j(x)= a_j,i(x), for 1 ≤ i,j ≤ N. Further, we assume that the coefficient a_ij, α_i, β are sufficiently smooth. 

For any x ∈Ω, σ∈ℝ^N, the operator D of the form (<ref>) is strongly elliptic if there exists a constant η such that

    ∑_i,j=1^N a_ij(x) σ_iσ_j≥η|σ|^2.


The operator (<ref>) on L^2(Ω) corresponds to elliptic operator in divergence form D written formally as

    Du = - div (a_ij∇ u ) + α·∇ u +  β u.

 This operator is associated with the bilinear form 

    a : H^1(Ω) × H^1(Ω) →ℤ

given by

    a(u,v) = ∫_Ω( ∑_i,j=1^N a_ij∂_iu ∂_jv̅ + α_j∂_ju v̅) + β u v̅ dx.

Let us moreover define the operator 

    𝒜u := - div(a_ij∇ u )+ α·∇ u -  div ((1- θ(·)) ∇ u) + β u.

The operator 𝒜 will be associated with the bilinear form 

    â : H^1(Ω) × H^1(Ω) →ℤ

given by

    â(u,v) = ∫_Ω( ∑_i,j=1^N a_ij∂_iu ∂_jv̅ + ∑_j=1^N( α_j∂_ju v̅ + (1- θ(·))u ∂_jv̅) + β u v̅)dx.



	The form â defined by (<ref>) is L^2(Ω)-elliptic.


Let us define â_1: H^1(Ω)× H^1(Ω) →ℤ is given by 

    â_1(u,v) = ã_1(u,v) + 	∫_Ω( ∑_j=1^N (1- θ(·))u ∂_jv̅) dx,
	
where ã_1(u,v): H^1(Ω)× H^1(Ω) →ℤ is defined as follows

    ã_1(u,v) = 	∫_Ω(  ∑_i,j=1^N a_ij∂_iu ∂_jv̅)dx.

The ellipticity condition (<ref>), the assumption (<ref>) with δ = η/2 S_N,2, and the Hölder inequality (<ref>), along with the Sobolev embedding theorem <ref>, are employed to obtain

    |â_1(u)|   = Reã_1(u)  +  c- T_Kc_L^N, ∞(Ω) u _L^2(Ω)∇ u _L^2(Ω)
       ≥η∇ u_L^2(Ω)^2- S_N,2 c-T_Kc_L^N, ∞(Ω)∇ u_L^2(Ω)^2
       ≥η/2∇ u _L^2(Ω)^2.

Thus, we have

    |â_1(u)|   + η/2 u_L^2(Ω)^2 ≥η/2 u_H^1(Ω)^2,

which implies that â_1 is L^2(Ω)-elliptic.
 
Let us define next â_2: H^1(Ω)× H^1(Ω) →ℤ by 

    â_2(u,v) = ∫_Ω(  ∑_j=1^Nα_j∂_ju v̅ + β u v̅)dx.

It is obvious that  â= â_1 + â_2.

The boundedness of the coefficients α_j, β and the Young's inequality imply that there exists a constant C such that	

    |â_2(u)|   ≤ C ∫_Ω|∇ u || u | dx + β_L^∞(Ω) u_L^2(Ω)^2
       ≤ C  u_H^1(Ω) u_L^2(Ω)
       ≤1/2( η u_H^1(Ω)^2 + 1/ηC^2  u_L^2(Ω)^2).

This leads to

    Re(â_1(u) + â_2(u)) + (  η/2 + C^2/2 η)  u_H^1(Ω)^2 ≥η/2 u _H^1(Ω)^2.

This concludes the proof. 










 §.§ Adapted Lie splitting method

In order to describe the construction of our proposed splitting method, let us first define 

    D̂u(t) = - D u(t) - ∇·( (1- θ(·)) u(t) ) ,

where θ(·) is given by (<ref>). As a consequence of Lemma <ref>, the operator D̂ remains a second-order strongly elliptic differential operator. The problem (<ref>) is then reformulated in terms of D̂ and θ as follows:

    {   ∂_t u(t) = D̂ u(t) + ∇· ( θ(·)  u(t))  +f(t),
        u(t)|_∂Ω =b(t),
        u(0) = u_0.	
     	.

We now proceed to the splitting of (<ref>) into the following two subproblems:

    {   ∂_t v(t) = D̂ v(t) + f(t),
       	v(t)|_∂Ω = b(t),
    	.

and 

    {   ∂_tw(t) = ∇·( θ(·) w(t) ),
         w(t)|_Γ = b_1(t).
    	.

We remark that in the splitting procedure, we impose a Dirichlet boundary condition for the diffusion equation (<ref>) while an inflow boundary condition is imposed for the pure convection equation (<ref>).


We define the iterative step of the adapted Lie splitting scheme with time step size τ>0 in the time interval [ t_n, t_n+1]  as follows:

    u_n+1 = ℒ_τ^Adu_n= φ_τ^D̂,f∘φ_τ^Ĉ(u_n),

where φ_t^D̂,h and φ_t^Ĉ are the solution operators of the problems (<ref>) and (<ref>), respectively. The operator ℒ_τ^Ad in (<ref>) is referred to as adapted Lie splitting scheme with time step size τ. 

Each step of the adapted Lie splitting scheme is explicitly carried out in the following way:






§ CONVERGENCE ANALYSIS

In this section we derive a thorough error analysis for adapted Lie splitting applied to (<ref>). More precisely, we prove that the adapted Lie splitting method is first-order convergent in the analytic semigroup setting and under appropriate assumptions.


 §.§ Analytical framework

The convergence analysis will be conducted within the context of analytic semigroups. To this end, we will rewrite the problem (<ref>) in a way that allows for the imposition of homogeneous Dirichlet boundary conditions. This transformation is only employed for the purposes of error analysis and does not affect our numerical scheme in general.
To achieve this, we introduce a smooth function, z(t), which is a solution to the following elliptic problem:

    {   D̂z(t) = 0,
       z(t)|_∂Ω=b(t),
    	.

where D̂ denotes a second-order strongly elliptic differential operator is given by (<ref>).  

Let us now define ũ(t)= u(t)-z(t) which satisfies the following abstract evolution equation

    {   ∂_tũ(t) +A ũ(t) = ∇· (θ(·) ũ(t)) + f(t),
        u(0)=u_0,
    	.

where the operator A is defined as: Aψ= -D̂ψ, for all ψ∈ D(A). For example, in the case of the operator D̂ on L^2(Ω), it can be shown that D(A)=H^2(Ω) ∩ H_0^1(Ω) ⊂ L^2(Ω). It is observed that the boundary condition for the equation <ref> is included in the domain of the operator A.

We make use of the following assumption for the operater A.

	Let X be a Banach space with norm ‖·‖. We assume that A is a linear operator on X and -A is the infinitesimal generator of an analytic semigroup e^-tA on X.

In the following we recall some properties of analytic semigroups which will be used extensively in the upcoming convergence analysis. For more background on the analytic framework, we refer to the books of Henry <cit.> or Pazy <cit.>. In particular, for any t ∈ [0,T] there exists a constant C such that

    ‖ e^-tA‖≤ C.

Moreover, for any α >0, the operator A enjoys the parabolic smoothing property

    A^α e^-tA≤C/t^α,

uniformly for t ∈ (0, T].



We now proceed by splitting the problem (<ref>) into two partial flows:

    ∂_tṽ(t)+ A ṽ(t)= f(t),

and 

    {   ∂_tw̃(t)= ∇·( θ(·) w̃(t)),
       w̃(t)|_Γ =0.
    	.

We note that the problem (<ref>) with the homeogeneous Dirichlet boundary conditions included in the domain of operator A.

The exact solution of (<ref>) at the time t_n+1 = t_n+τ is given by the variation-of-constants formula

    ũ(t_n+1)= e^-τ Aũ(t_n)+ ∫_0^τ e^-(τ-s)A( ∇· (θ(·)ũ(t_n+s) )  + f(t_n+s) ) ds.

This gives the exact solution of (<ref>) below

    u(t_n+1)   = z(t_n+1) +  e^-τ Aũ(t_n) + ∫_0^τ e^-(τ-s)A f(t_n+s)ds 
            + 	∫_0^τ e^-(τ-s)A∇·  (θ(·) ũ(t_n+s))  ds.

The upcoming subsections will investigate the local error and the global error of the adapted Lie splitting method applied to (<ref>).




 §.§ Local error

We first derive a bound for the local error for the adapted Lie splitting applied to (<ref>). For this purpose, let us consider one step of the numerical solution, starting at time t_n with the initial value w̃(t_n)= ũ(t_n) on the exact solution. The solution of (<ref>) with full step size τ can be expressed in a Taylor series with integral remainder as follows:

    w̃_n(τ)= w̃(t_n+τ) = ũ(t_n) + τ∇· (θ(·) ũ(t_n)) + ∫_0^τ(τ-s) w̃^''(t_n+ s) ds.

Next we integrate (<ref>) with initial data ṽ(t_n)= w̃_n(τ) to obtain

    ℒ_τ^Adũ(t_n)=	ṽ(t_n+τ)    = e^-τ Aũ(t_n)+ ∫_0^τ e^-(τ-s)Af(t_n+s) ds
            + τ e^-τ A∇· (θ(·) ũ(t_n))+ ∫_0^τ e^-τ A (τ-s) w̃^''(t_n+s)ds.

Therefore, the Lie splitting scheme applied to (<ref>) in the following gives the numerical solution

    u_n+1= ℒ_τ^Ad u(t_n)    =  z(t_n+1)+  e^-τ Aũ(t_n)+ ∫_0^τ e^-(τ-s)Af(t_n+s) ds
            +  τ e^-τ A∇· (θ(·) ũ(t_n))+ ∫_0^τ e^-τ A (τ-s) w̃^''(t_n+s)ds.

The local error d_n+1 = ℒ_τ^Ad(t_n)- u(t_n+1) is obtained by subtracting the exact solution (<ref>) from the numerical solution (<ref>). This gives the following representation of the local error:

    d_n+1    = τ e^-τ A∇· (θ(·) ũ(t_n)) - ∫_0^τ e^-(τ-s)A∇· (θ(·) ũ(t_n+s)) ds
            + ∫_0^τe^-τ A (τ -s) w̃^''(t_n+s) ds.

We set

    l_n(s) = e^-(τ-s)Al̃_n(s) ,   l̃_n(s)= ∇· (θ(·) ũ(t_n+s)).

The Taylor expansion of l_n(s) gives

    ∫_0^τ l_n(s)ds = τ l_n(0) + ∫_0^τ∫_0^σ l_n^'(σ) dσ ds.

This implies that

    d_n+1 = ∫_0^τ e^-τ A (τ -s) w̃^''(t_n+s) ds - ∫_0^τ∫_0^σ l_n^'(σ) dσ ds.

As a consequence of Theorem <ref>, the exact solution u(t_n) is continuously differentiable. Furthermore, by the boundedness property (<ref>) of an analytic semigroup e^-τ A, it can be shown that 

    ∫_0^τ e^-τ A (τ-s) w̃(t_n+s) ds ≤ C τ^2.

In order to bound d_n+1, it is sufficient to bound l_n^'(σ). The first derivative with respect to σ is given by

    l_n^'(σ) = A e^-(τ -σ)Al̃_n(σ) + e^-(τ-s)Al̃_n^'(σ).

By the definition of θ given by (<ref>), it holds that θ≤ K, for all fixed integers K >0. Using the boundedness property of the analytic semigroup (<ref>), the second term of the right hand side of (<ref>) is bounded. We observe that an additional A on the first term on the right hand side of (<ref>) can be compensated by applying the parabolic smoothing property (<ref>) in the error recursion for the global error (see the proof of Theorem <ref> below). Thus, the adapted Lie splitting has a consistency error proportional to τ. That is

    d_n+1 = A d̃_n+1+ 𝒪(τ^2),   d̃_n+1 = 𝒪(τ^2),

which completes the local error estimate for the adapted Lie splitting method applied to (<ref>).







 §.§ Global error

The first-order convergence of the adapted Lie splitting applied to (<ref>) is given in the next theorem. Before giving our main convergence analysis, we will make the following assumptions on the data of the problem (<ref>).

	Let D be a second-order strongly elliptic differential operator with smooth coefficients, f and b are twice continuously differentiable, and the convective coefficient c(x) satisfies the condition (<ref>). Moreover, we also assume that the initial data u_0 are spatially smooth and satisfy the boundary conditions.

As a consequence, under the aforementioned assumptions, Theorem <ref> implies that the solution u(t) of (<ref>) is continuously differentiable. This allows us to state our main convergence theorem.
 (Convergence of the adapted Lie splitting).
	Let the Assumptions <ref>-<ref> be in charge. Then there exists a constant τ_0>0 such that for all step sizes 0<τ≤τ_0 and t_n=nτ we have that the adapted Lie splitting (<ref>) applied to (<ref>) satisfies the global error bound
	
    ‖ u_n-u(t_n) ‖≤ Cτ(1 + | logτ| ),   0≤ nτ≤ T,

	where the constant C depends on T but is independent of τ and n.


Let us denote by e_n=u_n-u(t_n) the global error. Thus, we have

    e_n+1=ℒ_τ^Adu_n-ℒ_τ^Ad u(t_n)+d_n+1,
	
where d_n+1= ℒ_τ^Ad u(t_n)- u(t_n+1) denotes the local error. 

By (<ref>) and (<ref>), we have

    ℒ_τ^Adu_n-ℒ_τ^Ad u(t_n)
    	= e^-τ A e_n + 	τ R(u_n,u(t_n))+ 𝒪(τ^2),

where 

    R(u_n,u(t_n))=  e^-τ A( ∇· (θ(·) (u_n-u(t_n)) ) ).

Inserting (<ref>) and (<ref>) into (<ref>) yields

    e_n+1=d_n+1+τ R(u_n,u(t_n)) + 𝒪(τ^2),

As e_0 =0, solving the error recursion (<ref>) we get

    e_n   = ∑_k=0^n-1 e^-(n-k-1)τ A(d_k+1 + 𝒪(τ^2)) +  τ∑_k=0^n-1 e^-( n-k-1) τ A  R(u_k,u(t_k)).

The estimation of the local error ‖ d_k+1‖, given by (<ref>) allows us to obtain

    ‖ e_n‖   ≤∑_k=0^n-1 e^-(n-k-1)τ A (d_k+1 + 𝒪(τ^2))  +  τ∑_k=0^n-1 e^-( n-k-1) τ A  R(u_k,u(t_k)) 
       ≤∑_k=0^n-2 e^-(n-k-1)τ A (A 𝒪(τ^2)+ 𝒪(τ^2))   +  d_n +C τ∑_k=1^n-1  e^-(n-k)τ A∇· (θ(·) e_k).

Since d_n≤ Cτ, using the parabolic smoothing property (<ref>), we arrive at

    e_n   ≤ Cτ^2 ∑_k=0^n-2 e^-(n-k-1)τ A A + nC τ^2 + Cτ
            +C	τ∑_k=1^n-1 e^-(n-k)τ A A^1/2 A^-1/2∇ e_k
       ≤   Cτ^2 ∑_k=1^n-11/kτ + Cτ+ Cτ∑_k=1^n-11/( (n-k) τ)^1/2‖ e_k‖
       ≤ Cτ(1+ log| τ| )  + Cτ∑_k=1^n-1  t_n-k^-1/2‖ e_k‖.

This yields the desired bound by employing a discrete Gronwall lemma. This completes the proof.




	


§ NUMERICAL EXPERIMENTS IN TWO SPACE DIMENSION
 	
In this section, we present a series of numerical results for the convection-diffusion problem (<ref>) on the domain Ω = [-0.5,1]^2 , where the diffusion term is modeled by the Laplacian with the diffusion coefficient set to 0.01 and the source term f ≡ 0. The unbounded convective term is set to c(x, y) = 1/ | x | + 1/ | y | which belongs to the Lorentz spaces. The Laplacian is discretized by standard centered second-order finite differences, and the convection term is discretized by a first-order upwind scheme with 99 grid points with an equidistant mesh. The reference solution is computed using the ODE45 with absolute and relative tolerances are set to 10^-12. In all our simulations, the scheme (<ref>) is referred to as the classical Lie splitting, while the scheme proposed in Section 3.2 is referred to as the adapted Lie splitting.


 (time-independent boundary conditions) This problem considers the convection-diffusion equation with constant Dirichlet boundary conditions, which is given by u(t,x,-1/2)=u(t,x,1) =u(t,-1/2,y) = u(t,1,y)=1. We have chosen the initial data u(0,x, y)= 1+ sin(   2/3π( x+1/2) ) sin(   2/3π( y+1/2) ) that satisfies the presribed boundary conditions. Our numerical results are shown in Figure <ref>.
	

(time-dependent boundary conditions) In this example, we prescribe the convection-diffusion equation with time-dependent boundary conditions as follows: u(t,-1/2,y) = u(t,1,y) = 1 + sin(5t) and u(t,x,-1/2) = u(t,x,1) = 1 + sin(10t). As the initial data we have chosen u(0,x, y)= 1+ sin(   2/3π( x+1/2) ) sin(   2/3π( y+1/2) ). The numerical results are displayed in the Figure <ref>.
	

 	 				
As expected, Figure <ref> illustrates the instability effect observed in the classical Lie splitting (see the red line), whereas the adapted Lie splitting (see the magnetic line) appears to be perfectly stable in both cases with time-independent and time-dependent boundary conditions. Furthermore, our numerical simulations have demonstrated that for the time-independent boundary conditions, the adapted Lie splitting method achieves an accuracy factor of approximately 2 smaller than the classical Lie splitting method with the same time step size. In the case of time-dependent boundary conditions, the improved accuracy obtained by the adapted Lie splitting method is approximately a factor of 4 compared to the classical Lie splitting method. The lower accuracy is maintained when a smaller time step size is employed, as the two methods have an order of convergence of one. It is evident that the adapted Lie splitting method offers advantages over the traditional Lie splitting method, with greatly improved accuracy. 


	(homogeneous boundary conditions) In this example, we consider the convection-diffusion problem on the domain Ω = [-1,1]^2. The diffusion is modeled by the Laplacian with the confusion coefficient set to 0.1 and the convection is also set to c(x,y)=c(x,y) = 1/ | x | + 1/ | y |. We impose the equation with homogeneous Dirichlet boundary conditions. The initial data u_0(x,y) = sin( π x ) sin( π y ) is chosen to satisfy the boundary conditions. For the implementation, the equidistant mesh is used to discretize in space. Consequently, it is necessary to remove the element 0 in both the x-axis and the y-axis to avoid singularity (see Figure <ref>). Our numerical results are shown in Figure <ref>.
	
 	 	
It can be observed that the classical Lie splitting (the red line) is unstable as expected while the adapted Lis splitting (the magnetic line) appears to be perfectly stable (see Figure <ref>c). Furthermore, it can be seen that the improved accuracy obtained by the adapted Lie splitting method is approximately a factor of 18 over the classical Lie splitting method. 


§ NUMERICAL EXPERIMENTS IN THREE SPACE DIMENSION

In this section, we present numerical results for the convection-diffusion problem (<ref>) for the unknown u(t,x) in [0,T]× [-0.5,1]^3 of the form:

    {   ∂_tu-  u-  div{[ γx/| x |^2 + 1/2-N+ γ( 1- | x |^N-2-γ)x ] u}  = -  div(  x/| x |^N -γ) , 
        u|_∂Ω=0,
        u (0)=0,
    	.

where N/2 -1 < γ≤ N-1 with γ N-2. In the simulation, we have selected N=3 and γ = 2. The Laplacian is approximated by second-order finite differences methods, and the convection is discretized by a first-order upwind scheme with 19 grid points with an equidistant mesh. The numerical results are presented in Figure <ref>, which confirms that the classical Lie splitting (the red line) exhibits significant instability effects, in contrast to the remarkably improved stable behavior of the adapted Lie splitting (the magnetic line), which allows for large time step sizes. Furthermore, the adapted Lie splitting method shows superior accuracy compared to the classical Lie splitting.



 		








§ CONCLUSION AND OUTLOOK

We have proposed a newly adapted splitting approach that has proven successful in overcoming the instability encountered in convection-diffusion problems with unbounded convective coefficients. Our numerical experiments have demonstrated that the adapted Lie splitting is more accurate compared to the standard Lie splitting method. In addition, our proposed approach is potentially applicable to large systems of convection-diffusion, with a wider range of applicability in real-world phenomena. Extending our proposed approach to second-order convergent Strang splitting, or even higher-order splitting methods, has the potential to be beneficial, but it introduces additional complexities due to the unbounded property of the convective term. Such extensions may necessitate higher regularity assumptions on the solution. Furthermore, incompatibility effects may arise within the internal step of the Strang splitting scheme, resulting in order reduction (see, e.g., <cit.>). Nevertheless, despite the aforementioned complexities, the potential for enhanced convergence properties and computational efficiency of adapted splitting methods remains a compelling topic for future research.




*
siam



	
