A positivity-preserving, energy stable and convergent numerical scheme for the Poisson-Nernst-Planck system
Chun Liu, Cheng Wang, Steven M. Wise, Xingye Yue, Shenggao Zhou
AA positivity-preserving, energy stable and convergent numericalscheme for the Poisson-Nernst-Planck system
Chun Liu ∗ Cheng Wang † Steven M. Wise ‡ Xingye Yue § Shenggao Zhou ¶ September 18, 2020
Abstract
In this paper we propose and analyze a finite difference numerical scheme for the Poisson-Nernst-Planck equation (PNP) system. To understand the energy structure of the PNP model,we make use of the Energetic Variational Approach (EnVarA), so that the PNP system couldbe reformulated as a non-constant mobility H − gradient flow, with singular logarithmic energypotentials involved. To ensure the unique solvability and energy stability, the mobility functionis explicitly treated, while both the logarithmic and the electric potential diffusion terms aretreated implicitly, due to the convex nature of these two energy functional parts. The positivity-preserving property for both concentrations, n and p , is established at a theoretical level. Thisis based on the subtle fact that the singular nature of the logarithmic term around the value of0 prevents the numerical solution reaching the singular value, so that the numerical scheme isalways well-defined. In addition, an optimal rate convergence analysis is provided in this work,in which many highly non-standard estimates have to be involved, due to the nonlinear paraboliccoefficients. The higher order asymptotic expansion (up to third order temporal accuracy andfourth order spatial accuracy), the rough error estimate (to establish the (cid:96) ∞ bound for n and p ),and the refined error estimate have to be carried out to accomplish such a convergence result. Inour knowledge, this work will be the first to combine the following three theoretical propertiesfor a numerical scheme for the PNP system: (i) unique solvability and positivity, (ii) energystability, and (iii) optimal rate convergence. A few numerical results are also presented in thisarticle, which demonstrates the robustness of the proposed numerical scheme. Key words and phrases : Poisson-Nernst-Planck (PNP) system, logarithmic energy potential,positivity preserving, energy stability, optimal rate convergence analysis, higher order asymp-totic expansion
AMS subject classification : 35K35, 35K55, 65M12, 65M06, 82C70 ∗ Department of Applied Mathematics, Illinois Institute of Technology, IL 60616, USA ([email protected]) † Department of Mathematics, University of Massachusetts, North Dartmouth, MA 02747, USA([email protected]) ‡ Department of Mathematics, The University of Tennessee, Knoxville, TN 37996, USA ([email protected]) § Department of Mathematics, Soochow University, Suzhou 215006, P.R. China ([email protected]) ¶ Department of Mathematics, Soochow University, Suzhou 215006, P.R. China (Corresponding Author:[email protected]) a r X i v : . [ m a t h . NA ] S e p Introduction
We consider the two-particle Poisson-Nernst-Planck (PNP) system of equations ∂ t n = D n ∆ n − z e k B θ ∇ · ( D n n ∇ φ ) , (1.1) ∂ t p = D p ∆ p + z e k B θ ∇ · ( D p p ∇ φ ) , (1.2) − ε ∆ φ = z e ( p − n ) + ρ f , (1.3)where k B is the Boltzmann constant; θ is the absolute temperature; n and p are the concentrationsof negatively and positively charged ions, respectively; ε is the dielectric coefficient of the solution; z is valence of ions; e is the charge an electron; φ is the electric potential; and D n and D p arediffusion/mobility coefficients. Boundary conditions are very important for PNP systems and mustbe handled carefully [12]. However, we will assume periodic boundary conditions in this work forsimplicity of presentation. The analysis could be extended to more complicated, more physicalboundary conditions. In addition, for simplicity of presentation in the theoretical analysis, we as-sume that source term, ρ f , associated to the background fixed charge density, vanishes everywhere.The extension to a non-zero source term is straightforward.The PNP system is one of the most extensively studied models for the transport of chargedparticles in many physical and biological problems, including free electrons in semiconductors [19,25, 26]; fuel cells [29, 32]; ionic particles in electrokinetic fluids [2, 18, 24]; phase separation andpolarization for ionic liquids [14]; and ion channels in cell membranes [1, 10, 30]. The EnergeticVariational Approach (EnVarA) [9] shows that the PNP system is the gradient flow with respectto a particular free energy. In more detail, the free energy functional of a two-particle mixture maybe formulated as E ( n, p ) = (cid:90) Ω (cid:26) k B θ (cid:18) n ln nn + p ln pp (cid:19)(cid:27) d x + z e ε (cid:107) n − p (cid:107) H − , (1.4)under the assumption that n − p is of mean zero, where n and p are reference concentrations.The H − norm is defined via (cid:107) f (cid:107) H − := (cid:113) ( f, f ) H − , where ( f, g ) H − := ( ∇ ψ f , ∇ ψ g ) L , and ψ f ∈ ˚ H (Ω) := H (Ω) ∩ ˚ L (Ω) is the solution to − ∆ ψ f = f ∈ ˚ L (Ω) := (cid:8) f ∈ L (Ω) (cid:12)(cid:12) ( f, L = 0 (cid:9) . Formally, then (cid:107) f (cid:107) H − = (cid:0) f, ( − ∆) − f (cid:1) L . The PNP system (1.1) – (1.3) is the following H − -like gradient flow: ∂ t n = ∇ · (cid:18) D n k B θ n ∇ µ n (cid:19) , ∂ t p = ∇ · (cid:18) D p k B θ p ∇ µ p (cid:19) , (1.5)where µ n and µ p are chemical potentials given by µ n := δ n E = k B θ (ln nn + 1) + z e ε ( − ∆) − ( n − p ) = k B θ (ln nn + 1) − z e φ, (1.6) µ p := δ p E = k B θ (ln pp + 1) + z e ε ( − ∆) − ( p − n ) = k B θ (ln pp + 1) + z e φ, (1.7)2nd φ is the periodic and mean-zero solution to − ε ∆ φ = z e ( p − n ) . Of course, for the system to make sense, we require that the initial data satisfy1 | Ω | (cid:90) Ω n ( x , d x = 1 | Ω | (cid:90) Ω p ( x , d x > . Notice that non-constant coefficient mobility functions are involved in the formulated gradient flow.It is clear that the PDE solutions are conserved, positive (in the sense that n, p >
0, point-wise) and energy dissipative. There are a number of papers describing numerical methods for thePNP system. However, the theoretical analysis for numerical approximations turns out to be verychallenging, in particular for those based on the EnVarA formulation. First, the positivity of n and p have to be enforced to make the numerical scheme well-defined in the EnVarA formulation. Someexisting works have reported a positivity-preserving analysis [5, 6, 12, 17, 22, 20, 21, 34], whilemany of these analyses come from the maximum principle argument, instead of in the variationalframework. Second, the energy stability has also played a central role in the study of gradient flows.Such a stability analysis has appeared in a few existing numerical works [11, 23, 27], while the uniquesolvability and positivity-preserving analysis have been missing. Furthermore, there have been afew existing works for the convergence analysis [4, 31, 35], while these convergence estimates havebeen based on the perfect Laplacian operator structure for n and p , instead of the H − gradient flowstructure, so that the energy estimate is not available. Many other numerical schemes have beenreported [13, 16, 27, 28, 33, 36, 40]. However, no existing work has combined the following threetheoretical features in the numerical analyses: (i) unique solvability/positivity preserving property,(ii) energy stability in the variational framework, and (iii) optimal rate convergence analysis.In this paper we construct and analyze a finite difference numerical scheme, which preserves allthree important theoretical features. For the energy stability property, the numerical scheme has tobe based on the variational structure of the original PNP system. The mobility function is explicitlyupdated in the scheme to enforce the strictly elliptic nature of the operator associated with thetemporal derivative part in the H − gradient flow. For the chemical potential part, all the terms aretreated implicitly, because of the convex nature of both the logarithmic and the electric potentialdiffusion energy parts (in terms of n and p ). Moreover, the positivity-preserving property, for both n and p , will be theoretically established. Such an analysis is based on the fact that the numericalsolution is equivalent to the minimization of the numerical energy functional, and the singularnature of the logarithmic term around the value of 0 prevents the numerical solution reaching asingular value. As a result, the numerical scheme is always well-defined, and the unique solvabilityanalysis results from the convex nature of the implicit parts in the scheme. Such a techniquehas been successfully applied to the Cahn-Hilliard model [3, 7, 8], while its application to thePNP system will be involved more subtle details, due to the non-constant mobility. Furthermore,the energy stability comes directly from the corresponding convexity analysis, combined with thepositivity of the mobility functions.We provide an optimal rate convergence analysis for the proposed numerical scheme. The vari-ational structure and the non-constant mobility make this analysis highly challenging, especiallywhen compared with existing convergence estimates in [4, 31, 35], wherein a perfect Laplacian oper-ator is kept in tact. To overcome such a well-known difficulty, several highly non-standard estimateshave to be introduced, due to the nonlinear parabolic coefficients. The higher order asymptoticexpansion, up to the third order temporal accuracy and fourth order spatial accuracy, has to beperformed with a careful linearization technique. Such a higher order asymptotic expansion enables3ne to obtain a rough error estimate, so that to the (cid:96) ∞ bound for n and p could be derived. This (cid:96) ∞ estimate yields the upper and lower bounds of the two variables, and these bounds play a crucialrole in the subsequent analysis. Finally, the refined error estimate is carried out to accomplish thedesired convergence result. To our knowledge, it will be the first work to combine three theoreticalproperties for any numerical scheme for the PNP system: unique solvability/positivity-preserving,energy stability, and optimal rate convergence analysis.The rest of the article is organized as follows. In Section 2 we propose the fully discrete numericalscheme. The detailed proof for the positivity-preserving property of the numerical solution isprovided in Section 3, and the energy stability analysis is established in Section 4. The optimalrate convergence analysis is presented in Section 5. Some numerical results are provided in Section 6.Finally, the concluding remarks are given in Section 7. We introduce the dimensionless dependent variables ˆ n := n/n , ˆ p := p/p , with c = n = p , andˆ φ := φ/φ , with φ = k B θ z e . We use the dimensionless independent variables ˆ x := x/L and ˆ t := t/T , with L = (cid:115) εk B θ ( z e ) c and T = L D n . Define ˆ D := D p /D n . Then the dimensionless dynamical equations may be written (after droppingthe hats on the parameters and variables) as ∂ t n = ∇ · ( ∇ n − n ∇ φ ) , (2.1) ∂ t p = D ∇ · ( ∇ p + p ∇ φ ) , (2.2) − ∆ φ = p − n. (2.3)This system dissipates the dimensionless energy E ( n, p ) = (cid:90) Ω (cid:26) n ln n + p ln p + 12 ( n − p )( − ∆) − ( n − p ) (cid:27) d x , (2.4)and may be viewed as the following conserved gradient flow: ∂ t n = ∇ · ( n ∇ µ n ) , ∂ t p = D ∇ · ( p ∇ µ p ) , (2.5)where µ n and µ p are the dimensionless chemical potentials given by µ n := δ n E = ln n + 1 + ( − ∆) − ( n − p ) = ln n + 1 − φ, (2.6) µ p := δ p E = ln p + 1 + ( − ∆) − ( p − n ) = ln p + 1 + φ, (2.7)and φ is the periodic solution to − ∆ φ = p − n. Consequently, the energy is dissipated at the rate d t E = − (cid:90) Ω (cid:110) n |∇ µ n | + D p |∇ µ p | (cid:111) d x ≤ . .2 The finite difference spatial discretization We use the notation and results for some discrete functions and operators from [15, 38, 39]. LetΩ = ( − L x , L x ) × ( − L y , L y ) × ( − L z , L z ), where for simplicity, we assume L x = L y = L z =: L > N ∈ N be given, and define the grid spacing h := LN , i.e., a uniform spatial mesh size is takenfor simplicity of presentation. We define the following two uniform, infinite grids with grid spacing h > E := { p i + / | i ∈ Z } , C := { p i | i ∈ Z } , where p i = p ( i ) := ( i − / ) · h . Consider thefollowing 3-D discrete N -periodic function spaces: C per := { ν : C × C × C → R | ν i,j,k = ν i + αN,j + βN,k + γN , ∀ i, j, k, α, β, γ ∈ Z } , E xper := (cid:110) ν : E × C × C → R (cid:12)(cid:12)(cid:12) ν i + ,j,k = ν i + + αN,j + βN,k + γN , ∀ i, j, k, α, β, γ ∈ Z (cid:111) , in which identification ν i,j,k = ν ( p i , p j , p k ) is taken. The spaces E yper and E zper are analogouslydefined. The functions of C per are called cell centered functions . The functions of E xper , E yper , and E zper , are called east-west , north-south , and up-down face-centered functions , respectively. We alsodefine the mean zero space˚ C per := ν ∈ C per (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ν := h | Ω | N (cid:88) i,j,k =1 ν i,j,k , and denote (cid:126) E per := E xper ×E yper ×E zper . In addition, we introduce the important difference and averageoperators on the spaces: A x ν i + / ,j,k := 12 ( ν i +1 ,j,k + ν i,j,k ) , D x ν i + / ,j,k := 1 h ( ν i +1 ,j,k − ν i,j,k ) ,A y ν i,j + / ,k := 12 ( ν i,j +1 ,k + ν i,j,k ) , D y ν i,j + / ,k := 1 h ( ν i,j +1 ,k − ν i,j,k ) ,A z ν i,j,k + / := 12 ( ν i,j,k +1 + ν i,j,k ) , D z ν i,j,k + / := 1 h ( ν i,j,k +1 − ν i,j,k ) , with A x , D x : C per → E xper , A y , D y : C per → E yper , A z , D z : C per → E zper . Likewise, a x ν i,j,k := 12 (cid:0) ν i + / ,j,k + ν i − / ,j,k (cid:1) , d x ν i,j,k := 1 h (cid:0) ν i + / ,j,k − ν i − / ,j,k (cid:1) ,a y ν i,j,k := 12 (cid:0) ν i,j + / ,k + ν i,j − / ,k (cid:1) , d y ν i,j,k := 1 h (cid:0) ν i,j + / ,k − ν i,j − / ,k (cid:1) ,a z ν i,j,k := 12 (cid:0) ν i,j,k + / + ν i,j,k − / (cid:1) , d z ν i,j,k := 1 h (cid:0) ν i,j,k + / − ν i,j,k − / (cid:1) , with a x , d x : E xper → C per , a y , d y : E yper → C per , and a z , d z : E zper → C per . The discrete gradient ∇ h : C per → (cid:126) E per and the discrete divergence ∇ h · : (cid:126) E per → C per are given by ∇ h ν i,j,k = (cid:0) D x ν i + / ,j,k , D y ν i,j + / ,k , D z ν i,j,k + / (cid:1) , ∇ h · (cid:126)f i,j,k = d x f xi,j,k + d y f yi,j,k + d z f zi,j,k , where (cid:126)f = ( f x , f y , f z ) ∈ (cid:126) E per . The standard 3-D discrete Laplacian, ∆ h : C per → C per , becomes∆ h ν i,j,k := ∇ h · ( ∇ h ν ) i,j,k = d x ( D x ν ) i,j,k + d y ( D y ν ) i,j,k + d z ( D z ν ) i,j,k = 1 h ( ν i +1 ,j,k + ν i − ,j,k + ν i,j +1 ,k + ν i,j − ,k + ν i,j,k +1 + ν i,j,k − − ν i,j,k ) . D is a periodic scalar function that is defined at all of the face center points and (cid:126)f ∈ (cid:126) E per , then D (cid:126)f ∈ (cid:126) E per , assuming point-wise multiplication, and we may define ∇ h · (cid:0) D (cid:126)f (cid:1) i,j,k = d x ( D f x ) i,j,k + d y ( D f y ) i,j,k + d z ( D f z ) i,j,k . Specifically, if ν ∈ C per , then ∇ h · ( D∇ h · ) : C per → C per is defined point-wise via ∇ h · (cid:0) D∇ h ν (cid:1) i,j,k = d x ( D D x ν ) i,j,k + d y ( D D y ν ) i,j,k + d z ( D D z ν ) i,j,k . In addition, the following grid inner products are defined: (cid:104) ν, ξ (cid:105) := h N (cid:88) i,j,k =1 ν i,j,k ξ i,j,k , ν, ξ ∈ C per , [ ν, ξ ] x := (cid:104) a x ( νξ ) , (cid:105) , ν, ξ ∈ E xper , [ ν, ξ ] y := (cid:104) a y ( νξ ) , (cid:105) , ν, ξ ∈ E yper , [ ν, ξ ] z := (cid:104) a z ( νξ ) , (cid:105) , ν, ξ ∈ E zper . [ (cid:126)f , (cid:126)f ] := [ f x , f x ] x + [ f y , f y ] y + [ f z , f z ] z , (cid:126)f i = ( f xi , f yi , f zi ) ∈ (cid:126) E per , i = 1 , . Subsequently, we define the following norms for cell-centered functions. If ν ∈ C per , then (cid:107) ν (cid:107) := (cid:104) ν, ν (cid:105) ; (cid:107) ν (cid:107) pp := (cid:104)| ν | p , (cid:105) , for 1 ≤ p < ∞ , and (cid:107) ν (cid:107) ∞ := max ≤ i,j,k ≤ N | ν i,j,k | . The gradient norms areintroduced as follows: (cid:107)∇ h ν (cid:107) := [ ∇ h ν, ∇ h ν ] = [ D x ν, D x ν ] x + [ D y ν, D y ν ] y + [ D z ν, D z ν ] z , ∀ ν ∈ C per , (cid:107)∇ h ν (cid:107) pp := [ | D x ν | p , x + [ | D y ν | p , y + [ | D z ν | p , z , ∀ ν ∈ C per , ≤ p < ∞ . Higher order norms can be defined. For example, (cid:107) ν (cid:107) H h := (cid:107) ν (cid:107) + (cid:107)∇ h ν (cid:107) , (cid:107) ν (cid:107) H h := (cid:107) ν (cid:107) H h + (cid:107) ∆ h ν (cid:107) , ∀ ν ∈ C per . Lemma 2.1 ([37, 39]) . Let D be an arbitrary periodic, scalar function defined on all of the facecenter points. For any ψ, ν ∈ C per and any (cid:126)f ∈ (cid:126) E per , the following summation-by-parts formulasare valid: (cid:104) ψ, ∇ h · (cid:126)f (cid:105) = − [ ∇ h ψ, (cid:126)f ] , (cid:104) ψ, ∇ h · ( D∇ h ν ) (cid:105) = − [ ∇ h ψ, D∇ h ν ] . (2.8) For simplicity, we denote ( M mn ) i,j,k = n mi,j,k , ( M mp ) i,j,k = Dp mi,j,k , and introduce the followingmobility function at the face-centered mesh points:( ˘ M mn ) i + / ,j,k := A x ( M mn ) i + / ,j,k , ( ˘ M mn ) i,j + / ,k := A y ( M mn ) i,j + / ,k , (2.9)( ˘ M mn ) i,j,k + / := A z ( M mn ) i,j,k + / , with similar definitions for ˘ M mp . We use the following semi-implicit scheme: given n m , p m ∈ C per ,find n m +1 , p m +1 ∈ C per such that n m +1 − n m ∆ t = ∇ h · (cid:16) ˘ M mn ∇ h µ m +1 n (cid:17) , (2.10) p m +1 − p m ∆ t = ∇ h · (cid:16) ˘ M mp ∇ h µ m +1 p (cid:17) , (2.11) µ m +1 n = ln n m +1 + ( − ∆ h ) − ( n m +1 − p m +1 ) , (2.12) µ m +1 p = ln p m +1 + ( − ∆ h ) − ( p m +1 − n m +1 ) . (2.13)6 Positivity-preserving and unique solvability analyses
Recall the average operator: f = | Ω | (cid:104) f, (cid:105) . It is obvious that the numerical scheme (2.10) – (2.13)is mass conservative, so that n m = n := β , p m = p := β , with 0 < β , ∀ m ≥ . The following preliminary estimates, which are proved in the recent paper [3], are recalled. For any ϕ ∈ ˚ C per , there exists a unique ψ ∈ ˚ C per that solves L ˘ M ( ψ ) := −∇ h · ( ˇ M∇ h ψ ) = ϕ. (3.1)The following discrete norm may be defined: (cid:107) ϕ (cid:107) L − M = (cid:113) (cid:104) ϕ, L − M ( ϕ ) (cid:105) . (3.2)If ˇ M ≡
1, we have L ˘ M ( ψ ) = − ∆ h ψ and define (cid:107) ϕ (cid:107) − ,h = (cid:112) (cid:104) ϕ, ( − ∆ h ) − ( ϕ ) (cid:105) . (3.3) Lemma 3.1 ([3]) . Suppose that ϕ (cid:63) , ˆ ϕ ∈ C per , with ˆ ϕ − ϕ (cid:63) ∈ ˚ C per . Assume that < ˆ ϕ i,j,k , ϕ (cid:63)i,j,k ≤ M h , for all ≤ i, j, k ≤ N , where M h > may depend on h . The following estimate is valid: (cid:107) ( − ∆ h ) − ( ˆ ϕ − ϕ (cid:63) ) (cid:107) ∞ ≤ ˜ C M h , (3.4) where ˜ C > only depends on Ω . Lemma 3.2 ([3]) . Suppose that ϕ , ϕ ∈ C per , with ϕ − ϕ ∈ ˚ C per . Assume that (cid:107) ϕ (cid:107) ∞ , (cid:107) ϕ (cid:107) ∞ ≤ M h , and ˘ M ≥ M at a point-wise level, for some constant M > that is independent of h . Thenwe have the following estimate: (cid:13)(cid:13)(cid:13) L − M ( ϕ − ϕ ) (cid:13)(cid:13)(cid:13) ∞ ≤ C := ˜ C M − h − / , (3.5) where ˜ C > depends only upon M h and Ω . The positivity-preserving and unique solvability properties are established in the following the-orem.
Theorem 3.1.
Given n m , p m ∈ C per , with < n mi,j,k , p mi,j,k , ≤ i, j, k ≤ N , and n m − p m ∈ ˚ C per ,there exists a unique solution ( n m +1 , p m +1 ) ∈ [ C per ] to the numerical scheme (2.10) – (2.13) , with < n m +1 i,j,k , p m +1 i,j,k , ≤ i, j, k ≤ N and n m +1 − p m +1 ∈ ˚ C per .Proof. Suppose, as before, that n m = p m = β >
0. Define ν m := n m − β and ρ m := p m − β .The numerical solution of (2.10) – (2.13) is equivalent to the minimization of the following discreteenergy functional: J mh ( ν, ρ ) = 12∆ t (cid:32) (cid:107) ν − ν m (cid:107) L − M mn + (cid:107) ρ − ρ m (cid:107) L − M mp (cid:33) + (cid:104) ( ν + β ) ln( ν + β ) + ( ρ + β ) ln( ρ + β ) , (cid:105) + 12 (cid:107) ν − ρ (cid:107) − ,h , (3.6)7ver the admissible set˚ A h := (cid:26) ( ν, ρ ) ∈ (cid:104) ˚ C per (cid:105) (cid:12)(cid:12)(cid:12)(cid:12) < ν i,j,k + β , ρ i,j,k + β < M h , ≤ i, j, k ≤ N (cid:27) , (3.7)where M h := ( β | Ω | ) / h . We observe that J mh ( n, p ) is a strictly convex function over this domain.Next, we prove that there exists a minimizer of J mh ( n, p ) over the domain ˚ A h .Consider the following closed domain: for δ > A h,δ := (cid:26) ( ν, ρ ) ∈ (cid:104) ˚ C per (cid:105) (cid:12)(cid:12)(cid:12)(cid:12) δ ≤ ν i,j,k + β , ρ i,j,k + β ≤ M h − δ, ≤ i, j, k ≤ N (cid:27) . (3.8)Since ˚ A h,δ is a compact set in the hyperplane H := { ( ν, ρ ) | ν = ρ = 0 } , there exists a (not necessarilyunique) minimizer of J mh ( ν, ρ ) over ˚ A h,δ . The key point of the positivity analysis is that, such aminimizer could not occur at one of the boundary points (in H ) if δ is sufficiently small.Let us suppose that the minimizer of J mh ( ν, ρ ) occurs at a boundary point of ˚ A h,δ . Without lossof generality, we assume the minimizer is ( ν (cid:63)i,j,k , ρ (cid:63)i,j,k ), with ν (cid:63)i ,j ,k + β = δ , at some grid point( i , j , k ). Suppose that ν (cid:63) attains its maximum value at the point ( i , j , k ). By the fact that ν (cid:63) = 0, it is obvious that ν (cid:63)i ,j ,k ≥ ψ ∈ ˚ C per , d s J mh ( ν (cid:63) + sψ, ρ (cid:63) ) | s =0 = 1∆ t (cid:68) L − M mn ( ν (cid:63) − ν m ) , ψ (cid:69) + (cid:104) ln ( ν (cid:63) + β ) + 1 , ψ (cid:105) + (cid:10) ( − ∆ h ) − ( ν (cid:63) − ρ (cid:63) ) , ψ (cid:11) . Let us pick the direction ψ ∈ ˚ C per , such that ψ i,j,k = δ i,i δ j,j δ k,k − δ i,i δ j,j δ k,k , where δ k,(cid:96) is the Kronecker delta function. Then,1 h d s J mh ( ν (cid:63) + sψ, ρ (cid:63) ) | s =0 = ln (cid:32) ν (cid:63)i ,j ,k + β ν (cid:63)i ,j ,k + β (cid:33) + ( − ∆ h ) − ( ν (cid:63) − ρ (cid:63) ) i ,j ,k − ( − ∆ h ) − ( ν (cid:63) − ρ (cid:63) ) i ,j ,k + 1∆ t (cid:16) L − M mn ( ν (cid:63) − ν m ) i ,j ,k − L − M mn ( ν (cid:63) − ν m ) i ,j ,k (cid:17) . (3.9)Because n (cid:63)i ,j ,k = ν (cid:63)i ,j ,k + β = δ and n (cid:63)i ,j ,k = ν (cid:63)i ,j ,k + β ≥ β , we have ln (cid:32) ν (cid:63)i ,j ,k + β ν (cid:63)i ,j ,k + β (cid:33) ≤ ln δβ . (3.10)For the third and fourth terms appearing in (3.9), we apply Lemma 3.1 and obtain − C M h ≤ ( − ∆ h ) − ( ν (cid:63) − ρ (cid:63) ) i ,j ,k − ( − ∆ h ) − ( ν (cid:63) − ρ (cid:63) ) i ,j ,k ≤ C M h . (3.11)Similarly, for the last two terms appearing in (3.9), an application of Lemma 3.2 indicates that − C M − h − / ≤ L − M mn ( ν (cid:63) − ν m ) i ,j ,k − L − M mn ( ν (cid:63) − ν m ) i ,j ,k ≤ C M − h − / . (3.12)8onsequently, a substitution of (3.10) – (3.12) into (3.9) yields1 h d s J mh ( ν (cid:63) + sψ, ρ (cid:63) ) | s =0 ≤ ln δβ + 2 ˜ C M h + 2 ˜ C M − ∆ t − h − / . (3.13)Define D := 2 ˜ C M h + 2 ˜ C M − ∆ t − h − / , and note that D is a constant for fixed ∆ t and h , though it is singular, as ∆ t, h →
0. For anyfixed ∆ t and h , we may choose δ > δβ + D < . (3.14)This in turn gaurantees that d s J mh ( ν (cid:63) + sψ, ρ (cid:63) ) | s =0 < . (3.15)This contradicts the assumption that J mh has a minimum at ( ν (cid:63) , ρ (cid:63) ), since the directional derivativeis negative in a direction pointing into the interior of ˚ A h,δ .Using similar arguments, we can also prove that, the global minimum of J mh ( ν, ρ ) over ˚ A h,δ could not possibly occur at a boundary point satisfying ρ (cid:63)i ,j ,k + β = δ , if δ is small enough. Thedetails are left to interested readers.Therefore, the global minimum of J mh ( ν, ρ ) over ˚ A h,δ could only possibly occur at an interiorpoint, for δ > J mh ( ν, ρ ) is a smooth function, we conclude that there mustbe a solution ( ν i,j,k , ρ i,j,k ) ∈ ˚ A h , so that d s J mh ( ν + sψ, ρ + sφ ) | s =0 = 0 , ∀ ( ψ, φ ) ∈ (cid:104) ˚ C per (cid:105) , (3.16)which is equivalent to the numerical solution of (2.10) – (2.13). The existence of a positive numericalsolution is established.Finally, since J mh ( ν, ρ ) is a strictly convex function over ˚ A h , the uniqueness analysis for thisnumerical solution is straightforward. The proof of Theorem 3.1 is complete. With the positivity-preserving and unique solvability properties for the numerical scheme (2.10) –(2.13) established, we now prove energy stability. We introduce the following discrete energy: E h ( n, p ) := (cid:104) n ln n + p ln p, (cid:105) + 12 (cid:107) n − p (cid:107) − ,h . (4.1) Theorem 4.1.
For the numerical solution (2.10) – (2.13) , we have E h ( n m +1 , p m +1 ) + ∆ t (cid:0)(cid:2) ˇ M mn ∇ h µ m +1 n , ∇ h µ m +1 n (cid:3) + (cid:2) ˇ M mp ∇ h µ m +1 p , ∇ h µ m +1 p (cid:3)(cid:1) ≤ E h ( n m , p m ) , (4.2) so that E h ( n m , p m ) ≤ E h ( n , p ) ≤ C , for all m ∈ N , where C > is a constant independent of h .Proof. Taking discrete inner products of (2.10) with µ m +1 n and of (2.11) with µ m +1 p , we obtain (cid:10) n m +1 − n m , µ m +1 n (cid:11) + (cid:10) p m +1 − p m , µ m +1 p (cid:11) + ∆ t (cid:0)(cid:2) ˇ M mn ∇ h µ m +1 n , ∇ h µ m +1 n (cid:3) + (cid:2) ˇ M mp ∇ h µ m +1 p , ∇ h µ m +1 p (cid:3)(cid:1) = 0 . (4.3)9n the other hand, the convexity of the energy terms (cid:104) n ln n, (cid:105) , (cid:104) p ln p, (cid:105) and (cid:107) n − p (cid:107) − ,h implythat (cid:10) n m +1 − n m , ln n m +1 (cid:11) ≥ (cid:10) n m +1 ln n m +1 , (cid:11) − (cid:104) n m ln n m , (cid:105) , (4.4) (cid:10) p m +1 − p m , ln p m +1 (cid:11) ≥ (cid:10) p m +1 ln p m +1 , (cid:11) − (cid:104) p m ln p m , (cid:105) , (4.5) (cid:10) n m +1 − n m , ( − ∆ h ) − ( n m +1 − p m +1 ) (cid:11) + (cid:10) p m +1 − p m , ( − ∆ h ) − ( p m +1 − n m +1 ) (cid:11) ≥ (cid:16)(cid:13)(cid:13) n m +1 − p m +1 (cid:13)(cid:13) − ,h − (cid:107) n m − p m (cid:107) − ,h (cid:17) . (4.6)Substitution of (4.4) – (4.6) into (4.3) leads to (4.2), so that the unconditional energy stability isproved.Finally, that there is a constant C > h , such that E h ( n , p ) ≤ C ,follows from a consistency argument. The details are left to the interested reader. (cid:96) ∞ (0 , T ; (cid:96) ) ∩ (cid:96) (0 , T ; H h ) Now we proceed into the convergence analysis. Let ( N , P , Φ) be the exact PDE solution for thenon-dimensional PNP system (2.1) – (2.3). With sufficiently regular initial data, it is reasonableto assume that the exact solution has regularity of class R , where N , P ∈ R := H (0 , T ; C per (Ω)) ∩ H (cid:0) , T ; C (Ω) (cid:1) ∩ L ∞ (cid:0) , T ; C (Ω) (cid:1) . (5.1)In addition, we assume that the following separation property is valid for the exact solution: N ≥ (cid:15) , P ≥ (cid:15) , for some (cid:15) > , (5.2)which we assume holds at a point-wise level. Define N N ( · , t ) := P N N ( · , t ), P N ( · , t ) := P N P ( · , t ),the (spatial) Fourier projection of the exact solution into B K , the space of trigonometric polynomialsof degree to and including K (with N = 2 K + 1). The following projection approximation isstandard: if ( N , P ) ∈ L ∞ (0 , T ; H (cid:96) per (Ω)), for any (cid:96) ∈ N with 0 ≤ k ≤ (cid:96) , (cid:107) N N − N (cid:107) L ∞ (0 ,T ; H k ) ≤ Ch (cid:96) − k (cid:107) N (cid:107) L ∞ (0 ,T ; H (cid:96) ) , (cid:107) P N − P (cid:107) L ∞ (0 ,T ; H k ) ≤ Ch (cid:96) − k (cid:107) P (cid:107) L ∞ (0 ,T ; H (cid:96) ) . (5.3)Notice that the Fourier projection estimate (5.3) does not preserve the positivity of the variables,while we could take h sufficiently small (corresponding to a large N ) so that N N ≥ (cid:15) , P N ≥ (cid:15) .By N mN , P mN we denote N N ( · , t m ) and P N ( · , t m ), respectively, with t m = m · ∆ t . Since( N N , P N ) ∈ B K , the mass conservative property is available at the discrete level: N mN = 1 | Ω | (cid:90) Ω N N ( · , t m ) d x = 1 | Ω | (cid:90) Ω N N ( · , t m − ) d x = N m − N , (5.4) P mN = 1 | Ω | (cid:90) Ω P N ( · , t m ) d x = 1 | Ω | (cid:90) Ω P N ( · , t m − ) d x = P m − N , (5.5)for any m ∈ N . On the other hand, the solution of (2.10) – (2.11) is also mass conservative at thediscrete level: n m = n m − , p m = p m − , ∀ m ∈ N . (5.6)As indicated before, we use the mass conservative projection for the initial data: n = P h N N ( · , t =0), p = P h P N ( · , t = 0), that is( n ) i,j,k := N N ( p i , p j , p k , t = 0) , ( p ) i,j,k := P N ( p i , p j , p k , t = 0) . (5.7)10or the exact electric potential Φ, we denote its Fourier projection as Φ N . The error grid functionis defined as e mn := P h N mN − n m , e mp := P h P mN − p m , e mφ := P h Φ mN − φ m , ∀ m ∈ N . (5.8)Therefore, it follows that e mn = e mp = 0, for any m ∈ N , so that the discrete norm (cid:107) · (cid:107) − ,h is welldefined for the error grid function.The following theorem is the main result of this section. Theorem 5.1.
Given initial data N ( · , t = 0) , P ( · , t = 0) ∈ C (Ω) , suppose the exact solutionfor the PNP system (2.1) – (2.2) is of regularity class R . Then, provided ∆ t and h are sufficientlysmall, and under the linear refinement requirement C h ≤ ∆ t ≤ C h , we have (cid:107) e mn (cid:107) + (cid:107) e mp (cid:107) + (cid:16) ∆ t m (cid:88) k =1 ( (cid:107)∇ h e kn (cid:107) + (cid:107)∇ h e kp (cid:107) ) (cid:17) / + (cid:107) e mφ (cid:107) H h ≤ C (∆ t + h ) , (5.9) for all positive integers m , such that t m = m ∆ t ≤ T , where C > is independent of ∆ t and h . (2.10) – (2.13) : asymptotic expansionof the numerical solution By consistency, the project solution N N , P N solves the discrete equations (2.10) – (2.13) with afirst order accuracy in time and second order accuracy in space. Meanwhile, it is observed that thisleading local truncation error will not be enough to recover an a-priori (cid:96) ∞ bound for the numericalsolution to recover the separation property. To remedy this, we use a higher order consistencyanalysis, via a perturbation argument, to recover such a bound in later analysis. In more detail,we need to construct supplementary fields, N ∆ t, , N ∆ t, , P ∆ t, , P ∆ t, and ˇ N , ˇ P satisfyingˇ N = N N + P N (∆ t N ∆ t, + ∆ t N ∆ t, + h N h, ) , ˇ P = P N + P N (∆ t P ∆ t, + ∆ t P ∆ t, + h P h, ) , (5.10)so that a higher O (∆ t + h ) consistency is satisfied with the given numerical scheme (2.10) – (2.13).The constructed fields N ∆ t,j , N h, , P ∆ t,j , P h, which will be found using a perturbation expansion,will depend solely on the exact solution ( N , P ).The following truncation error analysis for the temporal discretization can be obtained by usinga straightforward Taylor expansion, as well as the estimate (5.3) for the projection solution: N m +1 N − N mN ∆ t = ∇ · (cid:16) N mN ∇ (ln N m +1 N + ( − ∆) − ( N m +1 N − P m +1 N )) (cid:17) +∆ t ( G (0) n ) m + O (∆ t ) + O ( h m ) , (5.11) P m +1 N − P mN ∆ t = ∇ · (cid:16) D P mN ∇ (ln P m +1 N + ( − ∆) − ( P m +1 N − N m +1 N )) (cid:17) +∆ t ( G (0) p ) m + O (∆ t ) + O ( h m ) . (5.12)Here m ≥ G (0) n G (0) p are smooth enough in the sense that their derivativesare bounded. 11he leading order temporal correction function ( N ∆ t, , P ∆ t, ) is given by solving the followingequations: ∂ t N ∆ t, = ∇ · (cid:16) N ∆ t, ∇ (ln N N + ( − ∆) − ( N N − P N ))+ N N ∇ ( 1 N N N ∆ t, + ( − ∆) − ( N ∆ t, − P ∆ t, )) (cid:17) − G (0) n , (5.13) ∂ t P ∆ t, = ∇ · (cid:16) D P ∆ t, ∇ (ln P N + ( − ∆) − ( P N − N N ))+ D P N ∇ ( 1 P N P ∆ t, + ( − ∆) − ( P ∆ t, − N ∆ t, )) (cid:17) − G (0) p . (5.14)Existence of a solution of the above linear PDE system is straightforward. Note that the solutiondepends only on the projection solution ( N N , P N ). In addition, the derivatives of ( N ∆ t, , P ∆ t, ) invarious orders are bounded. Of course, an application of the semi-implicit discretization (as givenby (5.11) – (5.12)) to (5.13) – (5.14) implies that N m +1∆ t, − N m ∆ t, ∆ t = ∇ · (cid:16) N m ∆ t, ∇ (ln N m +1 N + ( − ∆) − ( N m +1 N − P m +1 N ))+ N mN ∇ ( 1 N m +1 N N m +1∆ t, + ( − ∆) − ( N m +1∆ t, − P m +1∆ t, )) (cid:17) − ( G (0) n ) m + ∆ t h n + O (∆ t ) , (5.15) P m +1∆ t, − P m ∆ t, ∆ t = ∇ · (cid:16) D P m ∆ t, ∇ (ln P m +1 N + ( − ∆) − ( P m +1 N − N m +1 N ))+ D ( P N ) m ∇ ( 1 P m +1 N P m +1∆ t, + ( − ∆) − ( P m +1∆ t, − N m +1∆ t, )) (cid:17) − ( G (0) p ) m + ∆ t h n + O (∆ t ) . (5.16)Therefore, a combination of (5.11) – (5.12) and (5.15) – (5.16) leads to the second order temporaltruncation error for ˇ N := N N + ∆ t P N N ∆ t, , ˇ P := P N + ∆ t P N P ∆ t, :ˇ N m +11 − ˇ N m ∆ t = ∇ · (cid:16) ˇ N m ∇ (ln ˇ N m +11 + ( − ∆) − ( ˇ N m +11 − ˇ P m +11 )) (cid:17) +∆ t ( G (1) n ) m + O (∆ t ) + O ( h m ) , (5.17)ˇ P m +11 − ˇ P m ∆ t = ∇ · (cid:16) D ˇ P m ∇ (ln ˇ P m +11 + ( − ∆) − (ˇ P m +11 − ˇ N m +11 )) (cid:17) +∆ t ( G (1) p ) m + O (∆ t ) + O ( h m ) . (5.18)In the derivation of (5.17) – (5.18), the following linearized expansions have been utilized:ln ˇ N = ln( N N + ∆ t P N N ∆ t, ) = ln N N + ∆ t P N N ∆ t, N N + O (∆ t ) , (5.19)ln ˇ P = ln( P N + ∆ t P N P ∆ t, ) = ln P N + ∆ t P N P ∆ t, P N + O (∆ t ) . (5.20)Similarly, the next order temporal correction function ( N ∆ t, , P ∆ t, ) are given by the following12inear equations: ∂ t N ∆ t, = ∇ · (cid:16) N ∆ t, ∇ (ln ˇ N + ( − ∆) − ( ˇ N − ˇ P ))+ ˇ N ∇ ( 1ˇ N N ∆ t, + ( − ∆) − ( N ∆ t, − P ∆ t, )) (cid:17) − G (1) n , (5.21) ∂ t P ∆ t, = ∇ · (cid:16) D P ∆ t, ∇ (ln ˇ P + ( − ∆) − (ˇ P − ˇ N ))+ D ˇ P ∇ ( 1ˇ P P ∆ t, + ( − ∆) − ( P ∆ t, − N ∆ t, )) (cid:17) − G (1) p , (5.22)and the solution depends only on the exact solution ( N , P ), with derivatives of various orders staybounded. In turn, an application of the semi-implicit discretization to (5.21) – (5.22) implies that N m +1∆ t, − N m ∆ t, ∆ t = ∇ · (cid:16) N m ∆ t, ∇ (ln ˇ N m +11 + ( − ∆) − ( ˇ N m +11 − ˇ P m +11 ))+ ˇ N m ∇ ( 1ˇ N m +11 N m +1∆ t, + ( − ∆) − ( N m +1∆ t, − P m +1∆ t, )) (cid:17) − ( G (1) n ) m , (5.23) P m +1∆ t, − P m ∆ t, ∆ t = ∇ · (cid:16) D P m ∆ t, ∇ (ln ˇ P m +11 + ( − ∆) − (ˇ P m +11 − ˇ N m +11 ))+ D ˇ P m ∇ ( 1ˇ P m +11 P m +1∆ t, + ( − ∆) − ( P m +1∆ t, − N m +1∆ t, )) (cid:17) − ( G (1) p ) m . (5.24)Subsequently, a combination of (5.21) – (5.22) and (5.23) – (5.24) yields the third order temporaltruncation error for ˇ N := ˇ N + ∆ t P N N ∆ t, , ˇ P := ˇ P + ∆ t P N P ∆ t, :ˇ N m +12 − ˇ N m ∆ t = ∇ · (cid:16) ˇ N m ∇ (ln ˇ N m +12 + ( − ∆) − ( ˇ N m +12 − ˇ P m +12 )) (cid:17) +∆ t ( G (2) n ) m + O (∆ t ) + O ( h m ) , (5.25)ˇ P m +12 − ˇ P m ∆ t = ∇ · (cid:16) D ˇ P m ∇ (ln ˇ P m +12 + ( − ∆) − (ˇ P m +12 − ˇ N m +12 )) (cid:17) +∆ t ( G (2) p ) m + O (∆ t ) + O ( h m ) . (5.26)In fact, similar linearized expansions (as in (5.19) – (5.20)) have been used in the derivation.Next, we construct the spatial correction term ( N h, , P h, ) to upgrade the spatial accuracy order.The following truncation error analysis for the spatial discretization can be obtained by using astraightforward Taylor expansion for the constructed profile ( ˇ N , ˇ P ):ˇ N m +12 − ˇ N m ∆ t = ∇ h · (cid:16) A ( ˇ N m ) ∇ h (ln ˇ N m +12 + ( − ∆ h ) − ( ˇ N m +12 − ˇ P m +12 )) (cid:17) +∆ t ( G (2) n ) m + h ( H (0) n ) m + O (∆ t + h ) , (5.27)ˇ P m +12 − ˇ P m ∆ t = ∇ h · (cid:16) D A (ˇ P m ) ∇ h (ln ˇ P m +12 + ( − ∆ h ) − (ˇ P m +12 − ˇ N m +12 )) (cid:17) +∆ t ( G (2) p ) m + h ( H (0) p ) m + O (∆ t + h ) , (5.28)in which the average operator is taken in a similar form as (2.9). The spatially discrete functions H (0) n , H (0) p are smooth enough in the sense that their discrete derivatives are bounded. We alsonotice that there is no O ( h ) truncation error term, due to the fact that the centered difference usedin the spatial discretization gives local truncation errors with only even order terms, O ( h ), O ( h ),13tc. Subsequently, the spatial correction function ( N h, , P h, ) is given by solving the following linearPDE system: given by the following linear equations: ∂ t N h, = ∇ · (cid:16) N h, ∇ (ln ˇ N + ( − ∆) − ( ˇ N − ˇ P ))+ ˇ N ∇ ( 1ˇ N N h, + ( − ∆) − ( N h, − P h, )) (cid:17) − H (0) n , (5.29) ∂ t P h, = ∇ · (cid:16) D P h, ∇ (ln ˇ P + ( − ∆) − (ˇ P − ˇ N ))+ D ˇ P ∇ ( 1ˇ P P h, + ( − ∆) − ( P h, − N h, )) (cid:17) − H (0) p , (5.30)and the solution depends only on the exact solution ( N , P ), with the divided differences of variousorders stay bounded. In turn, an application of a full discretization to (5.29) – (5.30) implies that N m +1 h, − N mh, ∆ t = ∇ h · (cid:16) A ( N mh, ) ∇ h (ln ˇ N m +12 + ( − ∆ h ) − ( ˇ N m +12 − ˇ P m +12 ))+ A ( ˇ N m ) ∇ ( 1ˇ N m +12 N m +1 h, + ( − ∆ h ) − ( N m +1 h, − P m +1 h, )) (cid:17) − ( H (0) n ) m + O (∆ t + h ) , (5.31) P m +1 h, − P mh, ∆ t = ∇ h · (cid:16) D A ( P mh, ) ∇ h (ln ˇ P m +12 + ( − ∆ h ) − (ˇ P m +12 − ˇ N m +12 ))+ D A (ˇ P m ) ∇ h ( 1ˇ P m +12 P m +1 h, + ( − ∆ h ) − ( P m +1 h, − N m +1 h, )) (cid:17) − ( H (0) p ) m + O (∆ t + h ) . (5.32)Finally, a combination of (5.29) – (5.30) and (5.31) – (5.32) yields the higher order temporaltruncation error for ( ˇ N , ˇ P ) (as given by (5.10)):ˇ N m +1 − ˇ N m ∆ t = ∇ h · (cid:16) A ( ˇ N m ) ∇ h (ln ˇ N m +1 + ( − ∆ h ) − ( ˇ N m +1 − ˇ P m +1 )) (cid:17) + τ m +1 n , (5.33)ˇ P m +1 − ˇ P m ∆ t = ∇ h · (cid:16) D A (ˇ P m ) ∇ h (ln ˇ P m +1 + ( − ∆ h ) − (ˇ P m +1 − ˇ N m +1 )) (cid:17) + τ m +1 p , (5.34)where (cid:107) τ m +1 n (cid:107) , (cid:107) τ m +1 p (cid:107) ≤ C (∆ t + h ) . Again, the linear expansions have been extensively utilized.
Remark 5.1.
Trivial initial data N ∆ t,j ( · , t = 0) , P ∆ t,j ( · , t = 0) ≡ are given ( j = 1 , ) as in (5.13) – (5.14) , (5.21) – (5.22) , respectively. Similar trivial initial data is also imposed to ( N h, , P h, ) as in (5.29) – (5.30) . Therefore, using similar arguments as in (5.4) – (5.6) , we conclude that n ≡ ˇ N , p ≡ ˇ P , n k = n , p k = p , ∀ k ≥ , (5.35) and ˇ N k = 1 | Ω | (cid:90) Ω ˇ N ( · , t k ) d x = 1 | Ω | (cid:90) Ω ˇ N d x = n , ∀ k ≥ , (5.36)ˇ P k = 1 | Ω | (cid:90) Ω ˇ P ( · , t k ) d x = 1 | Ω | (cid:90) Ω ˇ P d x = p , ∀ k ≥ , (5.37)14 here the first step of (5.36) is based on the fact that ˇ N ∈ B K , and the second step comes from themass conservative property of ˇ N at the continuous level. These two properties will be used in lateranalysis.In addition, since ( ˇ N , ˇ P ) is mass conservative at a discrete level, as given by (5.36) – (5.37) , weobserve that the local truncation error τ n , τ p has a similar property: τ m +1 n = τ m +1 p = 0 , ∀ m ≥ . (5.38) Remark 5.2.
Since the temporal and spatial correction functions ( N ∆ t,j , P ∆ t,j ) , ( N h, , P h, ) arebounded, we recall the separation property (5.2) for the exact solution, and obtain a similar propertyfor the constructed profile ( ˇ N , ˇ P ) : ˇ N ≥ (cid:15) (cid:63) , ˇ P ≥ (cid:15) (cid:63) , for (cid:15) (cid:63) > , (5.39) in which the projection estimate (5.3) has been repeatedly used. Notice that we could take ∆ t and h sufficiently small so that (5.39) is valid for a modified value (cid:15) (cid:63) , such as (cid:15) (cid:63) = (cid:15) . Such a uniformbound will be used in the convergence analysis.In addition, since the correction functions only depend on ( N N , P N ) and the exact solution,its W , ∞ norm will stay bounded. In turn, we are able to obtain a discrete W , ∞ bound for theconstructed profile ( ˇ N , ˇ P ) : (cid:107) ˇ N k (cid:107) ∞ ≤ C (cid:63) , (cid:107) ˇ P k (cid:107) ∞ ≤ C (cid:63) , (cid:107)∇ h ˇ N k (cid:107) ∞ ≤ C (cid:63) , (cid:107)∇ h ˇ P k (cid:107) ∞ ≤ C (cid:63) , ∀ k ≥ . (5.40) Remark 5.3.
The reason for such a higher order asymptotic expansion and truncation error es-timate is to justify an a-priori (cid:96) ∞ bound of the numerical solution, which is needed to obtain theseparation property, similarly formulated as (5.39) for the constructed approximate solution. Withsuch a property valid for both the constructed approximate solution and the numerical solution, thenonlinear error term could be appropriately analyzed in the (cid:96) ∞ (0 , T ; (cid:96) ) convergence estimate. Instead of a direct analysis for the error function defined in (5.8), we introduce alternate numericalerror functions:˜ n m := P h ˇ N m − n m , ˜ p m := P h ˇ P m − p m , ˜ φ m := ( − ∆ h ) − (˜ p m − ˜ n m ) , ∀ m ∈ N . (5.41)The advantage of such a numerical error function is associated with its higher order accuracy, whichcomes from the higher order consistency estimate (5.33) – (5.34). Again, since ˜ n m = ˜ p m = 0, whichcomes from the fact (5.35) – (5.37), for any m ≥
0, we conclude that the discrete norm (cid:107) · (cid:107) − ,h iswell defined for the error grid function (˜ n m , ˜ p m ).In turn, subtracting the numerical scheme (2.10) – (2.13) from the consistency estimate (5.33)– (5.34) yields˜ n m +1 − ˜ n m ∆ t = ∇ h · (cid:16) A ( n m ) ∇ h ˜ µ m +1 n + A (˜ n m ) ∇ h V m +1 n (cid:17) + τ m +1 n , (5.42)˜ p m +1 − ˜ p m ∆ t = ∇ h · (cid:16) D A ( p m ) ∇ h ˜ µ m +1 p + D A (˜ p m ) ∇ h V m +1 p (cid:17) + τ m +1 p , (5.43)where ˜ µ m +1 n = ln ˇ N m +1 − ln n m +1 + ( − ∆ h ) − (˜ n m +1 − ˜ p m +1 ) , (5.44) V m +1 n = ln ˇ N m +1 + ( − ∆ h ) − ( ˇ N m +1 − ˇ P m +1 ) , (5.45)˜ µ m +1 p = ln ˇ P m +1 − ln p m +1 + ( − ∆ h ) − (˜ p m +1 − ˜ n m +1 ) , (5.46) V m +1 p = ln ˇ P m +1 + ( − ∆ h ) − (ˇ P m +1 − ˇ N m +1 ) . (5.47)15ince V m +1 n and V m +1 p only depend on the exact solution and the constructed profiles, we assumea discrete W , ∞ bound: (cid:107)V m +1 n (cid:107) W , ∞ h , (cid:107)V m +1 p (cid:107) W , ∞ h ≤ C (cid:63) . (5.48)To proceed with the nonlinear analysis, we make the following a-priori assumption at the previoustime step: (cid:107) ˜ n m (cid:107) , (cid:107) ˜ p m (cid:107) ≤ ∆ t + h . (5.49)Such an a-priori assumption will be recovered by the optimal rate convergence analysis at thenext time step, as will be demonstrated later. In turn, a discrete W , ∞ bound is available for thenumerical error function at the previous time step, with the help of inverse inequality: (cid:107) ˜ n m (cid:107) ∞ ≤ C (cid:107) ˜ n m (cid:107) h ≤ C (∆ t + h ) h ≤ C (∆ t + h ) ≤ , (5.50) (cid:107)∇ h ˜ n m (cid:107) ∞ ≤ C (cid:107) ˜ n m (cid:107) ∞ h ≤ C (∆ t + h ) h ≤ C (∆ t + h ) ≤ , (5.51)where the linear refinement constraint C h ≤ ∆ t ≤ C h has been used. By similar arguments, (cid:107) ˜ p m (cid:107) ∞ ≤ C (∆ t + h ) ≤ (cid:107)∇ h ˜ p m (cid:107) ∞ ≤ C (∆ t + h ) ≤ . (5.52)Subsequently, the following W , ∞ h bound is available for the numerical solution at the previous timestep: (cid:107) n m (cid:107) ∞ ≤ (cid:107) ˇ N m (cid:107) ∞ + (cid:107) ˜ n m (cid:107) ∞ ≤ ˜ C := C (cid:63) + 1 , (5.53) (cid:107) p m (cid:107) ∞ ≤ (cid:107) ˇ P m (cid:107) ∞ + (cid:107) ˜ p m (cid:107) ∞ ≤ ˜ C , (5.54) (cid:107)∇ h n m (cid:107) ∞ ≤ (cid:107)∇ h ˇ N m (cid:107) ∞ + (cid:107)∇ h ˜ n m (cid:107) ∞ ≤ C (cid:63) + 1 = ˜ C , (5.55) (cid:107)∇ h p m (cid:107) ∞ ≤ (cid:107)∇ h ˇ P m (cid:107) ∞ + (cid:107)∇ h ˜ p m (cid:107) ∞ ≤ C (cid:63) + 1 = ˜ C , (5.56)with the regularity assumption (5.40) applied. In addition, because of the (cid:96) ∞ estimate (5.50),(5.52) for the numerical error function, we can bound it by (cid:15) (cid:63) : (cid:107) ˜ n m (cid:107) ∞ ≤ C (∆ t + h ) ≤ (cid:15) (cid:63) (cid:107) ˜ p m (cid:107) ∞ ≤ C (∆ t + h ) ≤ (cid:15) (cid:63) , (5.57)so that the separation property is also valid for the numerical solution at the previous time step: n m ≥ ˇ N m − (cid:107) ˜ n m (cid:107) ∞ ≥ (cid:15) (cid:63) p m ≥ ˇ P m − (cid:107) ˜ p m (cid:107) ∞ ≥ (cid:15) (cid:63) , (5.58)where the separation estimate (5.39) has been utilized.Taking a discrete inner product with (5.42), (5.43) by ˜ µ n +1 n , ˜ µ n +1 p , respectively, leads to (cid:104) ˜ n m +1 , ˜ µ m +1 n (cid:105) + (cid:104) ˜ p m +1 , ˜ µ m +1 p (cid:105) + ∆ t ( (cid:104)A ( n m ) ∇ h ˜ µ m +1 n , ∇ h ˜ µ m +1 n (cid:105) + D (cid:104)A ( p m ) ∇ h ˜ µ m +1 p , ∇ h ˜ µ m +1 p (cid:105) )= (cid:104) ˜ n m , ˜ µ m +1 n (cid:105) + (cid:104) ˜ p m , ˜ µ m +1 p (cid:105) + ∆ t ( (cid:104) τ m +1 n , ˜ µ m +1 n (cid:105) + (cid:104) τ m +1 p , ˜ µ m +1 p (cid:105) ) − ∆ t ( (cid:104)A (˜ n m ) ∇ h V m +1 n , ∇ h ˜ µ m +1 n (cid:105) + D (cid:104)A (˜ p m ) ∇ h V m +1 p , ∇ h ˜ µ m +1 p (cid:105) ) . (5.59)Because of the separation estimate (5.58), at a point-wise level, the following inequalities are avail-able: (cid:104)A ( n m ) ∇ h ˜ µ m +1 n , ∇ h ˜ µ m +1 n (cid:105) ≥ (cid:15) (cid:63) (cid:107)∇ h ˜ µ m +1 n (cid:107) , (5.60) (cid:104)A ( p m ) ∇ h ˜ µ m +1 p , ∇ h ˜ µ m +1 p (cid:105) ≥ (cid:15) (cid:63) (cid:107)∇ h ˜ µ m +1 p (cid:107) . (5.61)16y the mean-free property (5.38) for the local truncation error terms, the following estimate canbe derived: (cid:104) τ m +1 n , ˜ µ m +1 n (cid:105) ≤ (cid:107) τ m +1 n (cid:107) − ,h · (cid:107)∇ h ˜ µ m +1 n (cid:107) ≤ (cid:15) (cid:63) (cid:107) τ m +1 n (cid:107) − ,h + 18 (cid:15) (cid:63) (cid:107)∇ h ˜ µ m +1 n (cid:107) , (5.62) (cid:104) τ m +1 p , ˜ µ m +1 p (cid:105) ≤ (cid:107) τ m +1 p (cid:107) − ,h · (cid:107)∇ h ˜ µ m +1 p (cid:107) ≤ D(cid:15) (cid:63) (cid:107) τ m +1 p (cid:107) − ,h + 18 D(cid:15) (cid:63) (cid:107)∇ h ˜ µ m +1 p (cid:107) . (5.63)For the two terms (cid:104) ˜ n m , ˜ µ m +1 n (cid:105) and (cid:104) ˜ p m , ˜ µ m +1 p (cid:105) , an application of the Cauchy inequality reveals that (cid:104) ˜ n m , ˜ µ m +1 n (cid:105) ≤ (cid:107) ˜ n m (cid:107) − ,h · (cid:107)∇ h ˜ µ m +1 n (cid:107) ≤ (cid:15) (cid:63) ∆ t (cid:107) ˜ n m (cid:107) − ,h + 18 (cid:15) (cid:63) ∆ t (cid:107)∇ h ˜ µ m +1 n (cid:107) , (5.64) (cid:104) ˜ p m , ˜ µ m +1 p (cid:105) ≤ (cid:107) ˜ p m (cid:107) − ,h · (cid:107)∇ h ˜ µ m +1 p (cid:107) ≤ D(cid:15) (cid:63) ∆ t (cid:107) ˜ p m (cid:107) − ,h + 18 D(cid:15) (cid:63) ∆ t (cid:107)∇ h ˜ µ m +1 p (cid:107) . (5.65)For the last two terms on the right hand side of (5.59), we see that −(cid:104)A (˜ n m ) ∇ h V m +1 n , ∇ h ˜ µ m +1 n (cid:105) ≤ (cid:107)∇ h V m +1 n (cid:107) ∞ · (cid:107)A (˜ n m ) (cid:107) · (cid:107)∇ h ˜ µ m +1 n (cid:107) ≤ C (cid:63) (cid:107) ˜ n m (cid:107) · (cid:107)∇ h ˜ µ m +1 n (cid:107) ≤ C (cid:63) ) (cid:15) (cid:63) (cid:107) ˜ n m (cid:107) + 18 (cid:15) (cid:63) (cid:107)∇ h ˜ µ m +1 n (cid:107) , (5.66)and, similarly, − D (cid:104)A (˜ p m ) ∇ h V m +1 p , ∇ h ˜ µ m +1 p (cid:105) ≤ C (cid:63) ) D(cid:15) (cid:63) (cid:107) ˜ p m (cid:107) + 18 D(cid:15) (cid:63) (cid:107)∇ h ˜ µ m +1 p (cid:107) . (5.67)A substitution of (5.60) – (5.67) into (5.59) leads to (cid:104) ˜ n m +1 , ˜ µ m +1 n (cid:105) + (cid:104) ˜ p m +1 , ˜ µ m +1 p (cid:105) + (cid:15) (cid:63) t ( (cid:107)∇ h ˜ µ m +1 n (cid:107) + D (cid:107)∇ h ˜ µ m +1 p (cid:107) ) ≤ (cid:15) (cid:63) ∆ t (cid:107) ˜ n m (cid:107) − ,h + 2 D(cid:15) (cid:63) ∆ t (cid:107) ˜ p m (cid:107) − ,h + 2∆ t(cid:15) (cid:63) (cid:107) τ m +1 n (cid:107) − ,h + 2∆ tD(cid:15) (cid:63) (cid:107) τ m +1 p (cid:107) − ,h + 2( C (cid:63) ) ( (cid:15) (cid:63) ) − ∆ t ( (cid:107) ˜ n m (cid:107) + D − (cid:107) ˜ p m (cid:107) ) . (5.68)Moreover, the detailed expansions in (5.44) and (5.46) reveal the following identities: (cid:104) ln ˇ N m +1 − ln n m +1 , ˜ n m +1 (cid:105) = (cid:104) ln ˇ N m +1 − ln n m +1 , ˇ N m +1 − n m +1 (cid:105) ≥ , (5.69) (cid:104) ln ˇ P m +1 − ln p m +1 , ˜ p m +1 (cid:105) = (cid:104) ln ˇ P m +1 − ln p m +1 , ˇ P m +1 − p m +1 (cid:105) ≥ , (5.70)and (cid:104) ( − ∆ h ) − (˜ n m +1 − ˜ p m +1 ) , ˜ n m +1 (cid:105) + (cid:104) ( − ∆ h ) − (˜ p m +1 − ˜ n m +1 ) , ˜ p m +1 (cid:105) = (cid:107) ˜ n m +1 − ˜ p m +1 (cid:107) − ,h ≥ , (5.71)where the positivities of ( n m +1 , p m +1 ) and ( ˇ N m +1 , ˇ P m +1 ) have been applied in the derivationof (5.69) and (5.70). Then we conclude that (cid:104) ˜ n m +1 , ˜ µ m +1 n (cid:105) + (cid:104) ˜ p m +1 , ˜ µ m +1 p (cid:105) ≥ . (5.72)17or the right hand side of (5.68), the following estimates are available, which come from the a-prioriassumption (5.49): 2 (cid:15) (cid:63) ∆ t (cid:107) ˜ n m (cid:107) − ,h ≤ C(cid:15) (cid:63) ∆ t (cid:107) ˜ n m (cid:107) ≤ C (∆ t + h ) , (5.73)2 D(cid:15) (cid:63) ∆ t (cid:107) ˜ p m (cid:107) − ,h ≤ CD(cid:15) (cid:63) ∆ t (cid:107) ˜ p m (cid:107) ≤ C (∆ t + h ) , (5.74)2∆ t(cid:15) (cid:63) (cid:107) τ m +1 n (cid:107) − ,h ≤ C ∆ t (cid:107) τ m +1 n (cid:107) ≤ C (∆ t + ∆ th ) , (5.75)2∆ tD(cid:15) (cid:63) (cid:107) τ m +1 p (cid:107) − ,h ≤ C ∆ t (cid:107) τ m +1 p (cid:107) ≤ C (∆ t + ∆ th ) , (5.76)2( C (cid:63) ) ( (cid:15) (cid:63) ) − ∆ t (cid:107) ˜ n m (cid:107) ≤ C ∆ t (cid:107) ˜ n m (cid:107) ≤ C (∆ t + h ) , (5.77)2( C (cid:63) ) ( (cid:15) (cid:63) ) − D − ∆ t (cid:107) ˜ p m (cid:107) ≤ C ∆ t (cid:107) ˜ p m (cid:107) ≤ C (∆ t + h ) , (5.78)where the fact that (cid:107) f (cid:107) − ,h ≤ C (cid:107) f (cid:107) , as well as the linear refinement constraint C h ≤ ∆ t ≤ C h ,have been repeatedly applied. Going back (5.68), we obtain (cid:15) (cid:63) t ( (cid:107)∇ h ˜ µ m +1 n (cid:107) + D (cid:107)∇ h ˜ µ m +1 p (cid:107) ) ≤ C (∆ t + h ) , (5.79)so that (cid:107)∇ h ˜ µ m +1 n (cid:107) , (cid:107)∇ h ˜ µ m +1 p (cid:107) ≤ C (∆ t + h ) . (5.80)Meanwhile, the error evolutionary equation (5.42) implies that (cid:107) ˜ n m +1 − ˜ n m (cid:107) ≤ ∆ t ( (cid:107)∇ h · ( A ( n m ) ∇ h ˜ µ m +1 n ) (cid:107) + (cid:107)∇ h · ( A (˜ n m ) ∇ h V m +1 n ) (cid:107) ) + ∆ t (cid:107) τ m +1 n (cid:107) . (5.81)Furthermore, the following estimate is available for the first term, based on a detailed nonlinearexpansion in the finite difference space, as well as repeated applications of discrete H¨older inequality: (cid:107)∇ h · ( A ( n m ) ∇ h ˜ µ m +1 n ) (cid:107) ≤ C ( (cid:107) n m (cid:107) ∞ · (cid:107)∇ h ∇ h ˜ µ m +1 n (cid:107) + (cid:107)∇ h n m (cid:107) ∞ · (cid:107)∇ h ˜ µ m +1 n (cid:107) ) ≤ C ˜ C ( (cid:107)∇ h ∇ h ˜ µ m +1 n (cid:107) + (cid:107)∇ h ˜ µ m +1 n (cid:107) ) ≤ C (∆ t + h ) , (5.82)in which the a-priori estimates (5.53), (5.55) have been used in the second step, and the followinginverse inequality has been applied in the last step: (cid:107)∇ h ∇ h ˜ µ m +1 n (cid:107) ≤ C (cid:107)∇ h ˜ µ m +1 n (cid:107) h ≤ C (∆ t + h ) . (5.83)The second term on the right hand side of (5.81) could be similarly analyzed: (cid:107)∇ h · ( A (˜ n m ) ∇ h V m +1 n ) (cid:107) ≤ C ( (cid:107) ˜ n m (cid:107) · (cid:107)∇ h ∇ h V m +1 n (cid:107) ∞ + (cid:107)∇ h ˜ n m (cid:107) · (cid:107)∇ h V m +1 n (cid:107) ∞ ) ≤ CC (cid:63) ( (cid:107) ˜ n m (cid:107) + (cid:107)∇ h ˜ n m (cid:107) ) ≤ C (∆ t + h ) , (5.84)in which the regularity assumption (5.48) has been recalled in the second step, while an inverseinequality (cid:107)∇ h f (cid:107) ≤ C (cid:107) f (cid:107) h has been applied in the last step. Therefore, a combination of (5.82),(5.84) and (5.81) results in (cid:107) ˜ n m +1 − ˜ n m (cid:107) ≤ C (∆ t + h ) + C (∆ t + h ) + C (∆ t + ∆ th ) ≤ C (∆ t + h ) . (5.85)18 similar estimate could be derived for (cid:107) ˜ p m +1 − ˜ p m (cid:107) : (cid:107) ˜ p m +1 − ˜ p m (cid:107) ≤ C (∆ t + h ) . (5.86)As a consequence, a combination with the a-priori error bound (5.49) (at the previous time step)results in a rough error estimate for ˜ n m +1 , ˜ p m +1 : (cid:107) ˜ n m +1 (cid:107) + (cid:107) ˜ p m +1 (cid:107) ≤ (cid:107) ˜ n m (cid:107) + (cid:107) ˜ p m (cid:107) + (cid:107) ˜ n m +1 − ˜ n m (cid:107) + (cid:107) ˜ p m +1 − ˜ p m (cid:107) ≤ ˆ C (∆ t + h ) , (5.87)under the linear refinement requirement C h ≤ ∆ t ≤ C h , with ˆ C dependent on the physicalparameters. Subsequently, an application of 3-D inverse inequality implies that (cid:107) ˜ n m +1 (cid:107) ∞ + (cid:107) ˜ p m +1 (cid:107) ∞ ≤ C ( (cid:107) ˜ n m +1 (cid:107) + (cid:107) ˜ p m +1 (cid:107) ) h ≤ ˆ C (∆ t + h ) , (5.88)where ˆ C := C ˆ C , under the same linear refinement requirement. Because of the accuracy order,we could take ∆ t and h sufficient small so thatˆ C (∆ t + h ) ≤ (cid:15) (cid:63) , (5.89)so that (cid:107) ˜ n m +1 (cid:107) ∞ + (cid:107) ˜ p m +1 (cid:107) ∞ ≤ (cid:15) (cid:63) . (5.90)Its combination with (5.39), the separation property for the constructed approximate solution,leads to a similar property for the numerical solution at time step t m +1 : (cid:15) (cid:63) ≤ n m +1 ≤ C (cid:63) + (cid:15) (cid:63) ≤ ˜ C and (cid:15) (cid:63) ≤ p m +1 ≤ C (cid:63) + (cid:15) (cid:63) ≤ ˜ C . (5.91)Such a uniform (cid:107) · (cid:107) ∞ bound will play a very important role in the refined error estimate. Remark 5.4.
In the rough error estimate (5.87) , we see that the accuracy order is lower thanthe one given by the a-priori-assumption (5.49) . Therefore, such a rough estimate could not beused for a global induction analysis. Instead, the purpose of such an estimate is to establish auniform (cid:107) · (cid:107) ∞ bound, via the technique of inverse inequality, so that a discrete separation propertybecomes available for the numerical solution, as well as its maximum values. With such a propertyestablished for the numerical solution, the refined error analysis will yield much sharper estimates. Before proceeding into the refined error estimate, the following two preliminary results are needed.
Lemma 5.1.
Under the a-priori (cid:107) · (cid:107) ∞ estimate (5.53) , (5.58) for the numerical solution at theprevious time step and the rough (cid:107) · (cid:107) ∞ estimate (5.91) for the one at the next time step, we have (cid:104)A ( n m ) ∇ h (ln ˇ N m +1 − ln n m +1 ) , ∇ h ˜ n m +1 (cid:105) ≥ γ (0) n (cid:107)∇ h ˜ n m +1 (cid:107) − M (0) n (cid:107) ˜ n m +1 (cid:107) − M (1) n h , (5.92) D (cid:104)A ( p m ) ∇ h (ln ˇ P m +1 − ln p m +1 ) , ∇ h ˜ p m +1 (cid:105) ≥ γ (0) p (cid:107)∇ h ˜ p m +1 (cid:107) − M (0) p (cid:107) ˜ p m +1 (cid:107) − M (1) p h , (5.93) where the constants γ (0) n , γ (0) p , M (0) n , M (0) p , M (1) n , M (1) p only depend on (cid:15) (cid:63) , C (cid:63) , ˜ C , D and | Ω | . roof. Looking at a single mesh cell ( i, j, k ) → ( i + 1 , j, k ), we make the following observation D x (ln ˇ N m +1 − ln n m +1 ) i + / ,j,k = 1 h (ln N m +1 i +1 ,j,k − ln N m +1 i,j,k ) − h (ln n m +1 i +1 ,j,k − ln n m +1 i,j,k )= 1 ξ N D x N m +1 i + / ,j,k − ξ n D x n m +1 i + / ,j,k = (cid:18) ξ N − ξ n (cid:19) D x N m +1 i + / ,j,k + 1 ξ n D x ˜ n m +1 i + / ,j,k , (5.94)in which the mean value theorem has been repeatedly applied, where ξ N is between N m +1 i +1 ,j,k and N m +1 i,j,k and ξ n is between n m +1 i +1 ,j,k and n m +1 i,j,k . (5.95)In turn, its product with D x ˜ n i + / ,j,k leads to D x ˜ n i + / ,j,k · D x (ln ˇ N m +1 − ln n m +1 ) i + / ,j,k = (cid:18) ξ N − ξ n (cid:19) D x N m +1 i + / ,j,k · D x ˜ n i + / ,j,k + 1 ξ n | D x ˜ n m +1 i + / ,j,k | . (5.96)For the second part, the rough (cid:107) · (cid:107) ∞ estimate (5.91) for n m +1 implies that 0 < ξ n ≤ ˜ C , which inturn gives 1 ξ n ≥ C and 1 ξ n | D x ˜ n m +1 i + / ,j,k | ≥ C | D x ˜ n m +1 i + / ,j,k | . (5.97)For the first term on the right hand side of (5.96), we begin with the following identity:1 ξ N = ln N m +1 i +1 ,j,k − ln N m +1 i,j,k N m +1 i +1 ,j,k − N m +1 i,j,k = ln (cid:16) N m +1 i +1 ,j,k − N m +1 i,j,k N m +1 i,j,k (cid:17) N m +1 i +1 ,j,k − N m +1 i,j,k . (5.98)By setting t (0) N = N m +1 i +1 ,j,k − N m +1 i,j,k N m +1 i,j,k , the following Taylor expansion is available:ln(1 + t (0) N ) = t (0) N −
12 ( t (0) N ) + 13 ( t (0) N ) −
14 ( t (0) N ) + 15(1 + η N ) ( t (0) N ) , (5.99)with η N between 0 and t (0) N . Its substitution into (5.98) yields1 ξ N = 1 N m +1 i,j,k − N m +1 i +1 ,j,k − N m +1 i,j,k N m +1 i,j,k ) + ( N m +1 i +1 ,j,k − N m +1 i,j,k ) N m +1 i,j,k ) − ( N m +1 i +1 ,j,k − N m +1 i,j,k ) N m +1 i,j,k ) + 15(1 + η N ) ( N m +1 i +1 ,j,k − N m +1 i,j,k ) ( N m +1 i,j,k ) . (5.100)A similar equality could be derived for ξ n :1 ξ n = 1 n m +1 i,j,k − n m +1 i +1 ,j,k − n m +1 i,j,k n m +1 i,j,k ) + ( n m +1 i +1 ,j,k − n m +1 i,j,k ) n m +1 i,j,k ) − ( n m +1 i +1 ,j,k − n m +1 i,j,k ) n m +1 i,j,k ) + 15(1 + η n ) ( n m +1 i +1 ,j,k − n m +1 i,j,k ) ( n m +1 i,j,k ) , (5.101)20ith η n between 0 and t (0) n = n m +1 i +1 ,j,k − n m +1 i,j,k n m +1 i,j,k . In addition, the following estimates are derived: (cid:12)(cid:12)(cid:12) N m +1 i,j,k − n m +1 i,j,k (cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12) ˜ n m +1 i,j,k N m +1 i,j,k n m +1 i,j,k (cid:12)(cid:12)(cid:12) ≤ (cid:15) (cid:63) ) | ˜ n m +1 i,j,k | , (5.102)and (cid:12)(cid:12)(cid:12) N m +1 i +1 ,j,k − N m +1 i,j,k ( N m +1 i,j,k ) − n m +1 i +1 ,j,k − n m +1 i,j,k ( n m +1 i,j,k ) (cid:12)(cid:12)(cid:12) ≤ (cid:12)(cid:12)(cid:12) ˜ n m +1 i +1 ,j,k − ˜ n m +1 i,j,k ( n m +1 i,j,k ) (cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12) ( N m +1 i +1 ,j,k − N m +1 i,j,k )( n m +1 i +1 ,j,k − N m +1 i,j,k )˜ n m +1 i,j,k ( N m +1 i,j,k ) ( n m +1 i,j,k ) (cid:12)(cid:12)(cid:12) ≤ (cid:15) (cid:63) ) ( | ˜ n m +1 i,j,k | + | ˜ n m +1 i +1 ,j,k | )+ 2 C (cid:63) ( C (cid:63) + ˜ C ) ( (cid:15) (cid:63) ) | ˜ n m +1 i,j,k |≤ Q (2) ( | ˜ n m +1 i,j,k | + | ˜ n m +1 i +1 ,j,k | ) , (5.103)where Q (2) := 4( (cid:15) (cid:63) ) + 8 C (cid:63) ( C (cid:63) + ˜ C )( (cid:15) (cid:63) ) , and the rough (cid:107) · (cid:107) ∞ estimate (5.91), the regularity assumption (5.40), and the separation prop-erty (5.39) have been extensively applied. The two other difference terms could be similarly ana-lyzed: (cid:12)(cid:12)(cid:12) ( N m +1 i +1 ,j,k − N m +1 i,j,k ) ( N m +1 i,j,k ) − ( n m +1 i +1 ,j,k − n m +1 i,j,k ) ( n m +1 i,j,k ) (cid:12)(cid:12)(cid:12) ≤ Q (3) ( | ˜ n m +1 i,j,k | + | ˜ n m +1 i +1 ,j,k | ) , (5.104) (cid:12)(cid:12)(cid:12) ( N m +1 i +1 ,j,k − N m +1 i,j,k ) ( N m +1 i,j,k ) − ( n m +1 i +1 ,j,k − n m +1 i,j,k ) ( n m +1 i,j,k ) (cid:12)(cid:12)(cid:12) ≤ Q (4) ( | ˜ n m +1 i,j,k | + | ˜ n m +1 i +1 ,j,k | ) , where Q (3) , Q (4) only depend on (cid:15) (cid:63) , C (cid:63) and ˜ C . For the remainder terms, we observe that | N m +1 i +1 ,j,k − N m +1 i,j,k | = h | D x N m +1 i + / ,j,k | ≤ h (cid:107) D x N m +1 (cid:107) ∞ ≤ C (cid:63) h, (5.105) | t (0) N | = (cid:12)(cid:12)(cid:12) N m +1 i +1 ,j,k − N m +1 i,j,k N m +1 i,j,k (cid:12)(cid:12)(cid:12) ≤ C (cid:63) ( (cid:15) (cid:63) ) − h ≤ Q (5) h ≤ , (5.106)where Q (5) = C (cid:63) ( (cid:15) (cid:63) ) − and where we have used (cid:15) (cid:63) ≤ N m +1 i,j,k . Furthermore | η N | ≤ , so that | η N | ≥
12 and (cid:12)(cid:12)(cid:12) η N ) (cid:12)(cid:12)(cid:12) ≤ . (5.107)Finally, |R | = (cid:12)(cid:12)(cid:12) η N ) ( N m +1 i +1 ,j,k − N m +1 i,j,k ) ( N m +1 i,j,k ) (cid:12)(cid:12)(cid:12) ≤ · ( C (cid:63) h ) ( (cid:15) (cid:63) ) ≤ Q (6) h , (5.108)21ith Q (6) = C (cid:63) ) (cid:15) (cid:63) ) . The other remainder term has a similar bound |R | = (cid:12)(cid:12)(cid:12) η n ) ( n m +1 i +1 ,j,k − n m +1 i,j,k ) ( n m +1 i,j,k ) (cid:12)(cid:12)(cid:12) ≤ · ( ˜ C h ) ( (cid:15) (cid:63) ) ≤ Q (7) h , (5.109)with Q (7) = C (cid:15) (cid:63) ) . Consequently, a combination of (5.102) – (5.105), (5.108) and (5.109) indicatesthat (cid:12)(cid:12)(cid:12) ξ N − ξ n (cid:12)(cid:12)(cid:12) ≤ Q (0) ( | ˜ n m +1 i,j,k | + | ˜ n m +1 i +1 ,j,k | ) + Q (cid:63) h , (5.110)with Q (0) = 2( (cid:15) (cid:63) ) + 12 Q (2) + 13 Q (3) + 14 Q (4) and Q (cid:63) = Q (6) + Q (7) . Then we arrive at an estimate for the first part on the right hand side of (5.96): (cid:18) ξ N − ξ n (cid:19) D x N m +1 i + / ,j,k · D x ˜ n i + / ,j,k ≥ − ( Q (0) ( | ˜ n m +1 i,j,k | + | ˜ n m +1 i +1 ,j,k | ) + Q (cid:63) h ) · C (cid:63) · | D x ˜ n i + / ,j,k |≥ − ( Q (0) ( | ˜ n m +1 i,j,k | + | ˜ n m +1 i +1 ,j,k | ) + Q (cid:63) h ) ( C (cid:63) ) ˜ C − (4 ˜ C ) − | D x ˜ n m +1 i + / ,j,k | . (5.111)Subsequently, a combination of (5.96), (5.96) and (5.111) results in D x ˜ n i + / ,j,k · D x (ln ˇ N m +1 − ln n m +1 ) i + / ,j,k ≥ − ( Q (0) ( | ˜ n m +1 i,j,k | + | ˜ n m +1 i +1 ,j,k | ) + Q (cid:63) h ) ( C (cid:63) ) ˜ C + 34 ˜ C | D x ˜ n m +1 i + / ,j,k | ≥
34 ˜ C | D x ˜ n m +1 i + / ,j,k | − Q (0) C (cid:63) ) ˜ C ( | ˜ n m +1 i,j,k | + | ˜ n m +1 i +1 ,j,k | ) − Q (cid:63) C (cid:63) ) ˜ C h . (5.112)Notice that this inequality is valid at a point-wise level. With summation over space, and keepingin mind of the a-priori (cid:107) · (cid:107) ∞ estimate (5.53), (5.58) for n m , we obtain (cid:104)A ( n m ) ∇ h (ln ˇ N m +1 − ln n m +1 ) , ∇ h ˜ n m +1 (cid:105)≥ (cid:15) (cid:63) ·
34 ˜ C (cid:107)∇ h ˜ n m +1 | − Q (0) C (cid:63) ) ˜ C (cid:107) ˜ n m +1 (cid:107) − Q (cid:63) C (cid:63) ) ˜ C | Ω | h . (5.113)This proves the first nonlinear estimate (5.92), by setting γ (0) n = (cid:15) (cid:63) C , M (0) n = 8( Q (0) C (cid:63) ) ˜ C , and M (1) n = 2( Q (cid:63) C (cid:63) ) ˜ C | Ω | . The second nonlinear estimate (5.93) could be derived exactly in the samemanner. The details are skipped for the sake of brevity.The next preliminary estimate is more straightforward. Lemma 5.2.
For ˜ φ k (for any k ≥ ) defined in (5.41) , we have the estimate (cid:107)∇ h ˜ φ k (cid:107) ≤ ˜ C (cid:107) ˜ n k − ˜ p k (cid:107) , (5.114) for some constant ˜ C > that is independent of h .Proof. Inequality (5.114) is a direct consequence of the standard estimate: (cid:107) f (cid:107) − ,h ≤ C (cid:107) f (cid:107) , forany f with f = 0. 22ow we proceed with the refined error estimate. Taking a discrete inner product with (5.42),(5.43) by 2˜ n m +1 , 2˜ p m +1 , respectively, leads to1∆ t ( (cid:107) ˜ n m +1 (cid:107) − (cid:107) ˜ n m (cid:107) + (cid:107) ˜ n m +1 − ˜ n m (cid:107) + (cid:107) ˜ p m +1 (cid:107) − (cid:107) ˜ p m (cid:107) + (cid:107) ˜ p m +1 − ˜ p m (cid:107) )+2( (cid:104)A ( n m ) ∇ h ˜ µ m +1 n , ∇ h ˜ n m +1 (cid:105) + D (cid:104)A ( p m ) ∇ h ˜ µ m +1 p , ∇ h ˜ p m +1 (cid:105) )= 2( (cid:104) τ m +1 n , ˜ n m +1 (cid:105) + (cid:104) τ m +1 p , ˜ p m +1 (cid:105) ) − (cid:104)A (˜ n m ) ∇ h V m +1 n , ∇ h ˜ n m +1 (cid:105) + D (cid:104)A (˜ p m ) ∇ h V m +1 p , ∇ h ˜ p m +1 (cid:105) ) , (5.115)where summation-by-parts has been applied. For the local truncation error terms, similar estimatescould be derived:2 (cid:104) τ m +1 n , ˜ n m +1 (cid:105) ≤ (cid:107) τ m +1 n (cid:107) + (cid:107) ˜ n m +1 (cid:107) , (cid:104) τ m +1 p , ˜ p m +1 (cid:105) ≤ (cid:107) τ m +1 p (cid:107) + (cid:107) ˜ p m +1 (cid:107) . (5.116)For the nonlinear diffusion error inner product on the left hand side, we see that (cid:104)A ( n m ) ∇ h ˜ µ m +1 n , ∇ h ˜ n m +1 (cid:105) = (cid:104)A ( n m ) ∇ h (ln ˇ N m +1 − ln n m +1 ) , ∇ h ˜ n m +1 (cid:105) + (cid:104)A ( n m ) ∇ h ˜ φ m +1 , ∇ h ˜ n m +1 (cid:105) . (5.117)The second part has the following lower bound (cid:104)A ( n m ) ∇ h ˜ φ m +1 , ∇ h ˜ n m +1 (cid:105) ≥ − ˜ C (cid:107)∇ h ˜ φ m +1 (cid:107) · (cid:107)∇ h ˜ n m +1 (cid:107) ≥ − ˜ C ˜ C (cid:107) ˜ n m +1 − ˜ p m +1 (cid:107) · (cid:107)∇ h ˜ n m +1 (cid:107) ≥ − ( ˜ C ˜ C ) (cid:107) ˜ n m +1 − ˜ p m +1 (cid:107) − γ (0) (cid:107)∇ h ˜ n m +1 (cid:107) , in which the inequality (5.114) (in Lemma 5.2) has been applied in the second step. Its substitutioninto (5.117), combined with the preliminary estimate (5.92) (in Lemma 5.1), leads to (cid:104)A ( n m ) ∇ h ˜ µ m +1 n , ∇ h ˜ n m +1 (cid:105)≥ γ (0) n (cid:107)∇ h ˜ n m +1 (cid:107) − M (0) n (cid:107) ˜ n m +1 (cid:107) − M (1) n h −
2( ˜ C ˜ C ) ( (cid:107) ˜ n m +1 (cid:107) + (cid:107) ˜ p m +1 (cid:107) ) . (5.118)A similar lower bound could be derived for the other nonlinear error inner product on the left handside; the details are skipped for the sake of brevity: D (cid:104)A ( p m ) ∇ h ˜ µ m +1 p , ∇ h ˜ p m +1 (cid:105)≥ γ (0) p (cid:107)∇ h ˜ p m +1 (cid:107) − M (0) p (cid:107) ˜ p m +1 (cid:107) − M (1) p h − D ˜ C ˜ C ) ( (cid:107) ˜ n m +1 (cid:107) + (cid:107) ˜ p m +1 (cid:107) ) . (5.119)For the last two nonlinear error inner product terms on the right hand side, a direct application ofCauchy inequality is applied: − (cid:104)A (˜ n m ) ∇ h V m +1 n , ∇ h ˜ n m +1 (cid:105) ≤ (cid:107)∇ h V m +1 n (cid:107) ∞ · (cid:107)A (˜ n m ) (cid:107) · (cid:107)∇ h ˜ n m +1 (cid:107) ≤ C (cid:63) (cid:107) ˜ n m (cid:107) · (cid:107)∇ h ˜ n m +1 (cid:107) ≤ C (cid:63) ) ( γ (0) n ) − (cid:107) ˜ n m (cid:107) + 12 γ (0) n (cid:107)∇ h ˜ n m +1 (cid:107) , (5.120) − D (cid:104)A (˜ p m ) ∇ h V m +1 p , ∇ h ˜ p m +1 (cid:105) ≤ DC (cid:63) ) ( γ (0) p ) − (cid:107) ˜ p m (cid:107) + 12 γ (0) p (cid:107)∇ h ˜ p m +1 (cid:107) , (5.121)with the regularity assumption (5.48) recalled. 23inally, a substitution of (5.116), (5.118) – (5.119) and (5.120) – (5.121) into (5.115) results in1∆ t ( (cid:107) ˜ n m +1 (cid:107) − (cid:107) ˜ n m (cid:107) + (cid:107) ˜ p m +1 (cid:107) − (cid:107) ˜ p m (cid:107) ) + γ (0) n (cid:107)∇ h ˜ n m +1 (cid:107) + γ (0) p (cid:107)∇ h ˜ p m +1 (cid:107) ≤ M (2) ( (cid:107) ˜ n m +1 (cid:107) + (cid:107) ˜ p m +1 (cid:107) ) + M (3) ( (cid:107) ˜ n m (cid:107) + (cid:107) ˜ p m (cid:107) ) + M (4) h + (cid:107) τ m +1 n (cid:107) + (cid:107) τ m +1 p (cid:107) , (5.122)where M (2) = 4( ˜ C ˜ C ) (1 + D ) + 2( M (0) n + M (0) p ) + 1 , (5.123) M (3) = 2( C (cid:63) ) ( γ (0) n ) − + 2( DC (cid:63) ) ( γ (0) p ) − , (5.124) M (4) = 2( M (1) n + M (1) p ) . (5.125)Therefore, an application of discrete Gronwall inequality leads to the desired higher order conver-gence estimate (cid:107) ˜ n m +1 (cid:107) + (cid:107) ˜ p m +1 (cid:107) + (cid:16) ∆ t m +1 (cid:88) k =1 ( (cid:107)∇ h ˜ n k (cid:107) + (cid:107)∇ h ˜ p k (cid:107) ) (cid:17) / ≤ C (∆ t + h ) , (5.126)based on the higher order truncation error accuracy, (cid:107) τ m +1 n (cid:107) , (cid:107) τ m +1 p (cid:107) ≤ C (∆ t + h ). Thiscompletes the refined error estimate. Recovery of the a-priori assumption (5.49)With the higher order error estimate (5.126) at hand, we notice that the a-priori assumptionin (5.49) is satisfied at the next time step t m +1 : (cid:107) ˜ n m +1 (cid:107) , (cid:107) ˜ p m +1 (cid:107) ≤ ˆ C (∆ t + h ) ≤ ∆ t + h , (5.127)provided ∆ t and h are sufficiently small. Therefore, an induction analysis could be applied. Thisfinishes the higher order convergence analysis.As a result, the convergence estimate (5.9) for the variable ( n, p ) is a direct consequenceof (5.126), combined with the definition (5.10) of the constructed approximate solution ( ˇ N , ˇ P ),as well as the projection estimate (5.3).In terms of the convergence estimate for the electric potential variable φ , we recall the definitionfor ˜ φ k in (5.41) and observe that (cid:107) ˜ φ m (cid:107) H h ≤ C (cid:107) ∆ h ˜ φ m (cid:107) ≤ Cε (cid:107) ˜ n m − ˜ p m (cid:107) ≤ ˆ C (∆ t + h ) , (5.128)where ˆ C = C ˆ C . Then (cid:107) ˜ φ m − e mφ (cid:107) H h ≤ C (cid:107) ∆ h ( ˜ φ m − e mφ ) (cid:107) ≤ ˆ C (∆ t + h ) , (5.129)and ( − ∆ h )( ˜ φ m − e mφ ) = P N (∆ t P ∆ t, + ∆ t P ∆ t, + h P h, − ∆ t N ∆ t, − ∆ t N ∆ t, − h N h, ) + τ mφ , (5.130)where the discrete elliptic regularity has been applied in (5.128), (5.129), and the truncation errorfor φ is defined as τ mφ = ( − ∆ h )Φ N − (ˇ P m − ˇ N m ).Finally, we arrive at (cid:107) e mφ (cid:107) H h ≤ (cid:107) ˜ φ m (cid:107) H h + (cid:107) ˜ φ m − e mφ (cid:107) H h ≤ ˆ C (∆ t + h ) + ˆ C (∆ t + h ) ≤ ( ˆ C + 1)(∆ t + h ) . (5.131)This completes the proof of Theorem 5.1. 24 Numerical results
To get numerical solutions, we need to solve the fully nonlinear scheme (2.10) – (2.13) at each timestep. We propose an iterative method as follows. First, the initial value for the nonlinear iterationis taken as n m +1 , := n m , p m +1 , := p m , and φ m +1 , := φ m . Subsequently, given the k -th iteratenumerical solution n m +1 ,k , p m +1 ,k , φ m +1 ,k , we obtain the first stage of the ( k + 1)-th iterate bysolving n m +1 , ∗ − ∆ t ∇ h · (cid:18) ˘ M mn ∇ h (cid:18) n m +1 , ∗ n m +1 ,k (cid:19)(cid:19) = n m + ∆ t ∇ h · (cid:16) ˘ M mn ∇ h (cid:16) ln n m +1 ,k − φ m +1 ,k (cid:17)(cid:17) ,p m +1 , ∗ − ∆ t ∇ h · (cid:18) ˘ M mp ∇ h (cid:18) p m +1 , ∗ p m +1 ,k (cid:19)(cid:19) = p m + ∆ t ∇ h · (cid:16) ˘ M mp ∇ h (cid:16) ln p m +1 ,k + φ m +1 ,k (cid:17)(cid:17) , − ∆ h φ m +1 , ∗ = p m +1 , ∗ − n m +1 , ∗ . (6.1)In addition, to make the nonlinear iteration smoother, we then obtain n m +1 ,k +1 , p m +1 ,k +1 , and φ m +1 ,k +1 by (cid:16) n m +1 ,k +1 , p m +1 ,k +1 , φ m +1 ,k +1 (cid:17) = ω r (cid:16) n m +1 ,k , p m +1 ,k , φ m +1 ,k (cid:17) + (1 − ω r ) (cid:0) n m +1 , ∗ , p m +1 , ∗ , φ m +1 , ∗ (cid:1) , (6.2)where ω r ∈ (0 ,
1) is a relaxation parameter. We notice that, two linear systems for n and p ,associated with M -matrices, need to be solved in the ( k + 1)-th iteration algorithm (6.1). In fact,(6.1) could be viewed as a linearized Newton iteration for the proposed numerical scheme (2.10) –(2.13), at least in the ln n and ln p nonlinear parts. It is expected that, under a sufficient conditionon the time step size ∆ t , such a linearized iteration algorithm guarantees positive concentrationsat a discrete level in each iteration stage, and an iteration convergence to the proposed numericalscheme (2.10) – (2.13) is also available. The detailed analysis will be left in the future works.In the following, we demonstrate the performance of the proposed numerical scheme in a twodimensional setting. With rescaling, the computational domain becomes Ω = ( − , . Also, wetake the parameters z = 1, n = p = c = 0 . L = 13 . D n = D p = D with D beingthe diffusion constant of sodium ions in water. To test accuracy, we consider the following exact solution n = e − t sin(2 πx ) cos(2 πy ) + 2 ,p = e − t cos(2 πx ) sin(2 πy ) + 2 ,ψ = e − t sin(2 πx ) sin(2 πy ) , (6.3)to the PNP equations with source terms: ∂ t n = ∇ · ( ∇ n − n ∇ φ ) + f n , (6.4) ∂ t p = ∇ · ( ∇ p + p ∇ φ ) + f p , (6.5) − ∆ φ = p − n + ρ f . (6.6)Here the source terms f n , f p , and ρ f , and the initial conditions are obtained with the known exactsolution. 25 (cid:96) ∞ error in p Order (cid:96) ∞ error in n Order (cid:96) ∞ error in ψ Order0.1 1.898E-2 - 1.898E-2 - 1.200E-1 -0.05 4.864E-3 1.96 4.864E-3 1.96 3.001E-2 2.000.025 1.231E-3 1.98 1.231E-3 1.98 7.524E-3 2.000.0125 3.093E-4 1.99 3.093E-4 1.99 1.882E-3 2.00Table 1: The (cid:96) ∞ error and convergence order for the numerical solutions of p , n , and ψ with∆ t = h .To verify the accuracy of the proposed scheme (2.10) to (2.13), we perform numerical tests usingvarious mesh resolution with ∆ t = h . Note that such a mesh ratio is chosen for the purpose ofaccuracy tests rather than the stability concern. As shown in Table 1, the (cid:96) ∞ error for numericalsolutions of p, n , and φ at time T = 0 . We also conduct numerical simulations to test the performance of the proposed scheme in preservingphysical properties at discrete level. The numerical schemes are applied to solve the rescaled PNPequations (6.4) to (6.6) without sources terms in the Nernst–Planck equations, but with a fixedcharge distribution given by ρ f ( x, y ) = e − [ ( x + ) +( y + ) ] − e − [ ( x + ) +( y − ) ] − e − [ ( x − ) +( y + ) ] + e − [ ( x − ) +( y − ) ] . The initial data for concentrations are given by p (0 , x, y ) = 0 . n (0 , x, y ) = 0 . . -2024 0.350.40.45 E h Total Mass of pT C M i n C M i n × -3 Figure 1: Left: The snapshots of ψ , p , and n at time T = 0 . T = 0 .
1, and T = 5. Right: Theevolution of discrete energy E h , total mass of p , and the minimum concentration C Min .Figure 1 displays snapshots of the electrostatic potential and concentrations at time T = 0 . T = 0 .
1, and T = 5. One observes that the concentrations of cations and anions develop peaks26nd valleys due to electrostatic interactions, and that the electrostatic potential initially inducedby the fixed charges gets screened quickly by attracted mobile ions carrying opposite charges, astime evolves. At T = 5, the system nearly reaches equilibrium.By periodic boundary conditions, the total mass of concentrations is conserved in time evolution.This is verified in the right panel of the Fig. 1, in which the total mass of the cations conversesperfectly. In addition, the discrete energy E h decreases monotonically, being consistent with ouranalysis; cf. Theorem 4.1. Of interest is the evolution of the minimum concentration that is definedby C Min := Min { Min i,j,k n mi,j,k , Min i,j,k p mi,j,k } . The evolution of C Min , together with the inset plot,demonstrates that the numerical solution of concentration remain positive all the time. In summary,our numerical tests further confirm that the proposed numerical scheme respects mass conservation,energy dissipation, and positivity at discrete level.
A finite difference numerical scheme is proposed and analyzed for the Poisson-Nernst-Planck (PNP)system. The Energetic Variational Approach (EnVarA) is taken, so that the PNP system couldbe reformulated as a non-constant mobility H − gradient flow, with singular logarithmic energypotentials involved. In the proposed numerical algorithm, the mobility function is explicitly treatedto ensure the unique solvability, while both the logarithmic and the electric potential diffusion termsare treated implicitly, because of their convex natures. The positivity-preserving property for both n and p are theoretically established, which is based on the subtle fact that, the singular nature of thelogarithmic term around the value of 0 prevents the numerical solution reaching the singular value.As a result, the numerical scheme is always well-defined. The energy stability of the numericalscheme comes from the convex nature of the energy functional in terms of n and p , combined withtheir positivity property. In addition, an optimal rate convergence analysis is provided in this work.To overcome a well-known difficulty associated with the non-constant mobility, many highly non-standard estimates have to be involved, due to the nonlinear parabolic coefficients. The higher orderasymptotic expansion, up to third order temporal accuracy and fourth order spatial accuracy, hasto be performed with a careful linearization technique. Such a higher order asymptotic expansionenable one to obtain a rough error estimate, so that to the (cid:96) ∞ bound for n and p could be derived.This (cid:96) ∞ estimate yields the upper and lower bounds of the two variables, and these bounds haveplayed a crucial role in the subsequent analysis. Finally, the refined error estimate are carriedout to accomplish the desired convergence result. It the first work to combine three theoreticalproperties for any numerical scheme to the PNP system: unique solvability/positivity-preserving,energy stability and optimal rate convergence analysis. A few numerical results are also presentedin this article, which demonstrates the robustness of the proposed numerical scheme. Acknowledgements
This work is supported in part by the National Science Foundation (USA) grants NSF DMS-1759535, NSF DMS-1759536 (C. Liu), NSF DMS-2012669 (C. Wang), NSF DMS-1719854, DMS-2012634 (S. Wise), National Natural Science Foundation of China 11971342 (X. Yue), 21773165,Young Elite Scientist Sponsorship Program by Jiangsu Association for Science and Technology,Natural Science Foundation of Jiangsu Province, China, and National Key R&D Program of China2018YFB0204404 (S. Zhou). 27 eferences [1] M.Z. Bazant, K. Thornton, and A. Ajdari. Diffuse-charge dynamics in electrochemical systems.
Phys. Rev. E , 70(2):021506, 2004.[2] Y. Ben and H.C. Chang. Nonlinear Smoluchowski slip velocity and micro-vortex generation.
J. Fluid Mech. , 461:229–238, 2002.[3] W. Chen, C. Wang, X. Wang, and S.M. Wise. Positivity-preserving, energy stable numericalschemes for the Cahn-Hilliard equation with logarithmic potential.
J. Comput. Phys.: X ,3:100031, 2019.[4] J. Ding, C. Wang, and S. Zhou. Optimal rate convergence analysis of a second order numericalscheme for the Poisson-Nernst-Planck system.
Numer. Math. Theor. Meth. Appl. , 12:607–626,2019.[5] J. Ding, Z. Wang, and S. Zhou. Positivity preserving finite difference methods for Poisson-Nernst-Planck equations with steric interactions: Application to slit-shaped nanopore conduc-tance.
J. Comput. Phys. , 397:108864, 2019.[6] J. Ding, Z. Wang, and S. Zhou. Structure-preserving and efficient numerical methods for iontransport.
J. Comput. Phys. , 2020.[7] L. Dong, C. Wang, H. Zhang, and Z. Zhang. A positivity-preserving, energy stable andconvergent numerical scheme for the Cahn-Hilliard equation with a Flory-Huggins-deGennesenergy.
Commun. Math. Sci. , 17:921–939, 2019.[8] L. Dong, C. Wang, H. Zhang, and Z. Zhang. A positivity-preserving second-order BDF schemefor the Cahn-Hilliard equation with variable interfacial parameters.
Commun. Comput. Phys. ,28:967–998, 2020.[9] B. Eisenberg, Y. Hyon, and C. Liu. Energy variational analysis of ions in water and channels:Field theory for primitive models of complex ionic fluids.
J. Chem. Phys. , 133(10):104104,2010.[10] R.S. Eisenberg. Computing the field in proteins and channels.
J. Mem. Biol. , 150:1–25, 1996.[11] A. Flavell, J. Kabre, and X. Li. An energy-preserving discretization for the Poisson-Nernst-Planck equations.
J. Comput. Electron , 16:431–441, 2017.[12] A. Flavell, M. Machen, R. Eisenberg, J. Kabre, C. Liu, and X. Li. A conservative finitedifference scheme for Poisson-Nernst-Planck equations.
J. Comput. Electron. , 13:235–249,2014.[13] H. Gao and D. He. Linearized conservative finite element methods for the Nernst–Planck–Poisson equations.
J. Sci. Comput. , 72:1269–1289, 2017.[14] N. Gavish and A. Yochelis. Theory of phase separation and polarization for pure ionic liquids.
J. Phys. Chem. Lett. , 7:1121–1126, 2016.[15] J. Guo, C. Wang, S.M. Wise, and X. Yue. An H convergence of a second-order convex-splitting, finite difference scheme for the three-dimensional Cahn-Hilliard equation. Commun.Math. Sci. , 14:489–515, 2016. 2816] D. He and K. Pan. An energy preserving finite difference scheme for the Poisson-Nernst-Plancksystem.
Appl. Math. Comput. , 287:214–223, 2015.[17] J. Hu and X. Huang. A fully discrete positivity-preserving and energy-dissipative finite differ-ence scheme for Poisson–Nernst–Planck equations.
Numer. Math. , 145:77–115, 2020.[18] R.J. Hunter.
Foundations of Colloid Science . Oxford University Press, Oxford, UK, 2001.[19] J.W. Jerome.
Analysis of Charge Transport. Mathematical Theory and Approximation of Semi-conductor Models . Springer-Verlag, New York, 1995.[20] H. Liu and W. Maimaitiyiming. Efficient, positive, and energy stable schemes for multi-dpoisson-nernst-planck systems.
ArXiv:2001.08350 , 2020.[21] H. Liu and W. Maimaitiyiming. Second order unconditional positivity-preserving schemes forsimulation of ionic channels.
Commun. Comput. Phys. , 2020.[22] H. Liu and Z. Wang. A free energy satisfying finite difference method for Poisson-Nernst-Planckequations.
J. Comput. Phys. , 268:363–376, 2014.[23] H. Liu and Z. Wang. A free energy satisfying discontinuous Galerkin method for one-dimensional Poisson-Nernst-Planck systems.
J. Comput. Phys. , 328:413–437, 2017.[24] J. Lyklema.
Fundamentals of Interface and Colloid Science. Volume II: Solid-liquid Interfaces .Academic Press Limited, San Diego, CA, 1995.[25] P.A. Markowich.
The Stationary Seminconductor Device Equations . Springer-Verlag, Vienna,Austria, 1986.[26] P.A. Markowich, C.A. Ringhofer, and C. Schmeiser.
Seminconductor Equations . Springer-Verlag, New York, 1990.[27] M. Metti, J. Xu, and C. Liu. Energetically stable discretizations for charge transport andelectrokinetic models.
J. Comput. Phys. , 306:1–18, 2016.[28] M. Mirzadeh and F. Gibou. A conservative discretization of the Poisson-Nernst-Planck equa-tions on adaptive cartesian grids.
J. Comput. Phys. , 274:633–653, 2014.[29] I. Nazarov and K. Promislow. The impact of membrane constraint on PEM fuel cell watermanagement.
J. Electrochem. Soc. , 154(7):623–630, 2007.[30] W. Nonner, D.P. Chen, and B. Eisenberg. Progress and prospects in permeation.
J. Gen.Physiol. , 113:773–782, 1999.[31] A. Prohl and M. Schmuck. Convergent discretizations for the Nernst–Planck–Poisson system.
Numer. Math. , 111:591–630, 2009.[32] K. Promislow and J.M. Stockie. Adiabatic relaxation of convective-diffusive gas transport ina porous fuel cell electrode.
SIAM J. Appl. Math. , 62(1):180–205, 2001.[33] Y. Qian, C. Wang, and S. Zhou. A positive and energy stable numerical scheme for thePoisson-Nernst-Planck-Cahn-Hilliard equations with steric interactions.
J. Comput. Phys. ,2020. Submitted and in review: arXiv:2002.09690.2934] F. Siddiqua, Z. Wang, and S. Zhou. A modified Poisson-Nernst-Planck model with excludedvolume effect: Theory and numerical implementation.
Commun. Math. Sci. , 16:251–271, 2018.[35] Y. Sun, P. Sun, B. Zheng, and G. Lin. Error analysis of finite element method for Poisson-Nernst-Planck equations.
J. Comput. Appl. Math. , 301:28–43, 2016.[36] B. Tu, M. Chen, Y. Xie, L. Zhang, B. Eisenberg, and B. Lu. A parallel finite element simulatorfor ion transport through three-dimensional ion channel systems.
J. Comput. Chem. , 287:214–223, 2015.[37] C. Wang and S.M. Wise. An energy stable and convergent finite-difference scheme for themodified phase field crystal equation.
SIAM J. Numer. Anal. , 49:945–969, 2011.[38] S.M. Wise. Unconditionally stable finite difference, nonlinear multigrid simulation of theCahn-Hilliard-Hele-Shaw system of equations.
J. Sci. Comput. , 44:38–68, 2010.[39] S.M. Wise, C. Wang, and J. Lowengrub. An energy stable and convergent finite-differencescheme for the phase field crystal equation.
SIAM J. Numer. Anal. , 47:2269–2288, 2009.[40] S. Xu, M. Chen, S. Majd, X. Yue, and C. Liu. Modeling and simulating asymmetrical con-ductance changes in gramicidin pores.