The Finite Element Approximation of the Nonlinear Poisson-Boltzmann Equation

0.36 MB·English

by  Long Chen

THE FINITE ELEMENT APPROXIMATION OF THE
NONLINEAR POISSON–BOLTZMANN EQUATION
LONGCHEN,MICHAELHOLST,ANDJINCHAOXU
ABSTRACT. A widely used electrostatics model in the biomolecular modeling com-
munity, the nonlinear Poisson–Boltzmann equation, along with its finite element ap-
proximation, are analyzed in this paper. A regularized Poisson–Boltzmann equation is
introduced as an auxiliary problem, making it possible to study the original nonlinear
0 equation with delta distribution sources. A priori error estimates for the finite element
1
approximation are obtained for the regularized Poisson–Boltzmann equation based on
0
certain quasi-uniform grids in two and three dimensions. Adaptive finite element ap-
2
proximationthroughlocalrefinementdrivenbyanaposteriorierrorestimateisshown
n
toconverge. ThePoisson–Boltzmannequationdoesnotappeartohavebeenpreviously
a
J studiedindetailtheoretically,anditishopedthatthispaperwillhelpprovidemolecular
8 modelerswithabetterfoundationfortheiranalyticalandcomputationalworkwiththe
Poisson–Boltzmann equation. Note that this article apparently gives the first rigorous

] convergence result for a numerical discretization technique for the nonlinear Poisson–
A
Boltzmannequationwithdeltadistributionsources,anditalsointroducesthefirstprov-
N
ably convergent adaptive method for the equation. This last result is currently one of
. onlyahandfulofexistingconvergenceresultsofthistypefornonlinearproblems.
h
t
a
m
[

1 CONTENTS
v
0 1. Introduction 2
5
2. ThePoisson–Boltzmannequation 3
3
3. Regularizationofthecontinuousproblem 5
1
. 4. Existenceanduniqueness 7
1 5. ContinuousaprioriL∞-estimates 9
0
6. FiniteelementmethodsfortheregularizedPoisson–Boltzmannequation 11
0
1 6.1. Quasi-optimalapriorierrorestimate 12
: 6.2. DiscreteaprioriL∞-estimates 13
v
7. Convergenceofadaptivefiniteelementapproximation 14
i
X 7.1. Aposteriorierrorestimate 14
r 7.2. Markingandrefinementstrategy 16
a
7.3. Convergenceanalysis 17
8. Summaryandconcludingremarks 19
References 21
Date:November20,2006.
Keywordsandphrases. nonlinearPoisson–Boltzmannequation,finiteelementmethods,apriorianda
posteriorierrorestimate,convergenceofadaptivemethods.
The first author was supported in part by NSF awards 0411723 and 022560, in part by DOE awards
DE-FG02-04ER25620andDE-FG02-05ER25707,andinpartbyNIHawardP41RR08605.
ThesecondauthorwassupportedinpartbyNSFawards0411723and022560,inpartbyDOEawards
DE-FG02-04ER25620andDE-FG02-05ER25707,andinpartbyNIHawardP41RR08605.
The third author was supported in part by NSF DMS0308946, DMS-0619587, DMS-0609727, and
NSFC-10528102.
1
2 L.CHEN,M.HOLST,ANDJ.XU
1. INTRODUCTION
In this paper, we shall design and analyze finite element approximations of a widely
usedelectrostaticsmodelinthebiomolecularmodelingcommunity,thenonlinearPoisson–
Boltzmannequation(PBE):
(cid:88)Nm
−∇·(ε∇u˜)+κ¯2sinh(u˜) = q δ in Rd, d = 2,3, (1.1)
i i
i=1
where the dielectric ε and the modified Debye–Hu¨ckel parameter κ¯ are piecewise con-
stantsindomainsΩ (thedomainforthebiomoleculeofinterest)andΩ (thedomainfor
m s
asolventsurroundingthebiomolecule),andδ := δ(x−x )isaDiracdistributionatpoint
i i
x . The importance of (1.1) in biomolecular modeling is well-established; cf. [14, 43]
i
for thorough discussions. Some analytical solutions are known, but only for unrealistic
structure geometries, and usually only for linearizations of the equation; cf. [29] for a
collection of these solutions and for references to the large amount of literature on ana-
lytical solutions to the PBE and similar equations. The current technological advances
are more demanding and require the solution of highly nonlinear problems in compli-
cated geometries. To this end, numerical methods, including the finite element method,
arewidelyusedtosolvethenonlinearPBE[29,30,5,6,44,19,56].
The main difficulties for the rigorous analysis and provably good numerical approx-
imation of solutions to the nonlinear Poisson–Boltzmann equation include: (1) Dirac
distribution sources, (2) exponential rapid nonlinearities, and (3) discontinuous coeffi-
cients. We shall address these difficulties in this paper. To deal with the δ distribution
sources,wedecomposeu˜ asanunknownfunctioninH1 andaknownsingularfunction,
namely,
(cid:88)Nm
u˜ = u+G, with G = G ,
i
i=1
where G is the fundamental solution of −ε ∆G = q δ in Rd. Substituting this de-
i m i i i
composition into the PBE, we then obtain the so-called regularized Poisson–Boltzmann
equation(RPBE):
−∇·(ε∇u)+κ¯2sinh(u+G) = ∇·((ε−ε )∇G) in Rd, d = 2,3.
m
The singularities of the δ distributions are transferred to G, which then exhibits degen-
erate behavior at each {x } ⊂ Ω . At those points, both sinhG(x ) and ∇G(x ) exhibit
i m i i
blowup. However, since G is known analytically, one avoids having to build numerical
approximationstoG. Moreover,bothofthecoefficientsκ¯ andε−ε arezeroinsideΩ
m m
where the blowup behavior arises. Due to this cutoff nature of coefficients, we obtain a
well-defined nonlinear second-order elliptic equation for the regularized solution u with
a source term in H−1. We will show that it also admits a unique solution u ∈ H1, even
thoughtheoriginalsolutionu˜ ∈/ H1 duetothesingularitiespresentinG.
Singular function expansions are a common technique in applied and computational
mathematics for this type of singularity; this type of expansion has been previously pro-
posedforthePoisson–Boltzmannequationin[58]andwasshown(empirically)toallow
for more accurate finite difference approximations. In their work, the motivation for the
technique was the poor discrete approximation of arbitrarily placed delta distributions
usingonlythefixedcornersofuniformfinitedifferencemeshes. Inthepresentwork,our
interestisindevelopingfiniteelementmethodsusingcompletelyunstructuredmeshes,so
weareabletoplacethedeltadistributionspreciselywheretheyshouldbeanddonothave
FINITEELEMENTAPPROXIMATIONOFTHENONLINEARPOISSON-BOLTZMANNEQUATION 3
thisproblemwithapproximatedeltafunctionplacement. Ourmotivationhereforconsid-
ering a singular function expansion is rather that the solution to the Poisson–Boltzmann
equation is simply not smooth enough to either analyze or approximate using standard
methodswithoutusingsomesortoftwo-scaleormultiscaleexpansionthatrepresentsthe
nonsmooth part of the solution analytically. In fact, it will turn out that expanding the
solutionintothesumofthreefunctions,namely,aknownsingularfunction,anunknown
solution to a linear auxiliary problem, and an unknown solution to a second nonlinear
auxiliary problem, is the key to establishing some fundamental results and estimates for
thecontinuousproblemandisalsothekeytodevelopingacompleteapproximationthe-
ory for the discrete problem as well as provably convergent nonadaptive and adaptive
numericalmethods.
Startingwithsomebasicresultsonexistence,uniqueness,andaprioriestimatesforthe
continuousproblem,weanalyzethefiniteelementdiscretizationandderivediscreteana-
loguesofthecontinuousresultstoshowthatdiscretizationleadstoawell-poseddiscrete
problem. Usingmaximumprinciplesforthecontinuousanddiscreteproblems,wederive
a priori L∞-estimates for the continuous and discrete solutions to control the nonlinear-
ity, allowing us to obtain a priori error estimates for our finite element approximation of
theform
(cid:107)u−u (cid:107) (cid:46) inf (cid:107)u−v (cid:107) ,
h 1 h 1
v ∈Vh
h D
whereVh isthelinearfiniteelementsubspacedefinedoverquasi-uniformtriangulations
D
withacertainboundarycondition,andu isthefiniteelementapproximationofuinVh.
h D
The result is quasi-optimal in the sense it implies that the finite element approximation
totheRPBEiswithinaconstantofbeingthebestapproximationfromthesubspaceVh.
D
After establishing these results for finite element approximations, we describe an adap-
tive approximation algorithm that uses mesh adaptation through local refinement driven
byaposteriorierrorestimates. Theadaptivealgorithmcanbeviewedasamechanismfor
dealing with the primary remaining difficulty in the RPBE, namely, the discontinuities
of the coefficients across the interface between the solvent and the molecular regions.
Finally, we shall prove that our adaptive finite element method will produce a sequence
of approximations that converges to the solution of the continuous nonlinear PBE. This
last result is one of only a handful of existing results of this type for nonlinear elliptic
equations(theothersbeing[24,48,15]).
The outline of this paper is as follows. In section 2, we give a brief derivation and
overviewofthePoisson–Boltzmannequation. Insection3,wederivearegularizedform
ofthePoisson–Boltzmannequationbyusingasingularfunctionexpansion. Insection4,
we give some basic existence and uniqueness results for the RPBE. In section 5, we de-
riveanaprioriL∞-estimateforthecontinuousproblem. Afterintroducingfiniteelement
methods for the RPBE, in section 6 we derive an analogous a priori L∞-estimate for
the discrete problem, and based on this we obtain a quasi-optimal a priori error estimate
for the finite element approximation. In section 7, we describe the adaptive algorithm,
present an a posteriori error estimate, and prove a general convergence result for the al-
gorithm. In the last section, we summarize our work and give further remarks on the
practicalaspectsusingresultsinthepresentpaper.
2. THE POISSON–BOLTZMANN EQUATION
In this section we shall give a brief introduction to the nonlinear Poisson–Boltzmann
equation. Adetailedderivationcanbefoundin[47,29].
4 L.CHEN,M.HOLST,ANDJ.XU
The nonlinear PBE, a second-order nonlinear partial differential equation, is funda-
mental to Debye–Hu¨ckel continuum electrostatic theory [22]. It determines a dimen-
sionless potential around a charged biological structure immersed in a salt solution. The
PBE arises from the Gauss law, represented mathematically by the Poisson equation,
whichrelatestheelectrostaticpotentialΦinadielectrictothechargedensityρ:
−∇·(ε∇Φ) = ρ,
whereεisthedielectricconstantofthemediumandhereistypicallypiecewiseconstant.
Usuallyitjumpsbyoneortwoordersofmagnitudeattheinterfacebetweenthecharged
structure (a biological molecular or membrane) and the solvent (a salt solution). The
charge density ρ consist of two components: ρ = ρ +ρ . For the macromolecule,
macro ion
the charge density is a summation of δ distributions at N point charges in the point
m
chargebehavior,i.e.,
(cid:88)Nm 4πe2
ρ (x) = q δ(x−x ), q = cz ,
macro i i i i
κ T
B
i=1
where κ > 0 is the Boltzmann constant, T is the temperature, e is the unit of charge,
B c
andz istheamountofcharge.
i
For the mobile ions in the solvent, the charge density ρ cannot be given in a de-
ion
terministic way. Instead it will be given by the Boltzmann distribution. If the solvent
contains N types of ions, of valence Z and of bulk concentration c , then a Boltzmann
i i
assumptionabouttheequilibriumdistributionoftheionsleadsto
N (cid:18) (cid:19)
(cid:88) e Φ
c
ρ = c Z e exp −Z .
ion i i c i
κ T
B
i=1
Forasymmetric1 : 1electrolyte,N = 2,c = c ,andZ = (−1)i,whichyields
i 0 i
(cid:18) (cid:19)
e Φ
c
ρ = −2c e sinh .
ion 0 c
κ T
B
We can now write the PBE for modeling the electrostatic potential of a solvated bi-
ological structure. Let us denote the molecule region by Ω ⊂ Rd and consider the
m
solvent region Ω = Rd\Ω¯ . We use u˜ to denote the dimensionless potential and κ¯2 to
s m
denote the modified Debye–Hu¨ckel parameter (which is a function of the ionic strength
ofthesolvent). ThenonlinearPoisson–Boltzmannequationisthen
(cid:88)Nm
−∇·(ε∇u˜)+κ¯2sinh(u˜) = q δ inRd, (2.1)
i i
i=1
u˜(∞) = 0, (2.2)
where
(cid:26) (cid:26)
ε if x ∈ Ω , 0 if x ∈ Ω ,
ε = m m and κ¯ = √ m
ε if x ∈ Ω , ε κ > 0 if x ∈ Ω .
s s s s
It has been determined empirically that ε ≈ 2 and ε ≈ 80. The structure itself (e.g.,
m s
a biological molecule or a membrane) is represented implicitly by ε and κ¯, as well as
explicitly by the N point charges q = z e at the positions x . The charge positions
m i i c i
are located in the strict interior of the molecular region Ω . A physically reasonable
m
FINITEELEMENTAPPROXIMATIONOFTHENONLINEARPOISSON-BOLTZMANNEQUATION 5
mathematical assumption is that all charge locations obey the following lower bound on
theirdistancetothesolventregionΩ forsomeσ > 0:
s
|x−x | ≥ σ ∀x ∈ Ω , i = 1,...,N . (2.3)
i s m
In some models employing the PBE, there is a third region Ω (the Stern layer [11]), a
l
layer between Ω and Ω . In the presence of a Stern layer, the parameter σ in (2.3)
m s
increasesinvalue. Ouranalysisandresultscanbeeasilygeneralizedtothiscaseaswell.
Some analytical solutions of the nonlinear PBE are known, but only for unrealistic
structure geometries and usually only for linearizations of the equation; cf. [29] for a
collection of these solutions and for references to the large amount of literature on an-
alytical solutions to the PBE and similar equations. However, the problem is highly
nonlinear. Surface potentials of the linear and the nonlinear PBE differ by over an order
of magnitude [44]. Hence, using the nonlinear version of the PBE model is fundamen-
tallyimportanttoaccuratelydescribephysicaleffects,andaccesstoreliableandaccurate
numerical approximation techniques for the nonlinear PBE is critically important in this
researcharea.
Wefinishthissectionbymakingsomeremarksaboutanalternativeequivalentformu-
lation of the PBE. It is well known (cf. [47, 29]) that the PBE is formally equivalent to
acouplingoftwoequationsfortheelectrostaticpotentialindifferentregionsΩ andΩ
m s
through the boundary interface. This equivalence can be rigorously justified. Inside Ω ,
m
therearenoions. ThustheequationissimplythePoissonequation
(cid:88)Nm
−∇·(ε ∇u˜) = q δ inΩ .
m i i m
i=1
In the solvent region Ω , there are no atoms. Thus the density is given purely by the
s
Boltzmanndistribution
−∇·(ε ∇u˜)+κ¯2sinh(u˜) = 0 inΩ .
s s
These two equations are coupled together through the boundary conditions on the inter-
faceΓ := ∂Ω = ∂Ω ∩Ω :
m s m
(cid:20) (cid:21)
∂u˜
[u˜] = 0, and ε = 0,
Γ
∂n
Γ Γ
where [f]| = lim f(x+tn )−f(x−tn ), with n being the unit outward normal
Γ t→0 Γ Γ Γ
directionofinterfaceΓ. WewillassumeΓtobesufficientlysmooth,say,ofclassC2.
Solving the individual subdomain systems and coupling them through the boundary,
in the spirit of a nonoverlapping domain decomposition method, is nontrivial due to the
complicated boundary conditions and subdomain shapes. Approaches such as mortar-
based finite element methods to solve the coupled equations for linear or nonlinear PBE
canbefoundin[19,51].
3. REGULARIZATION OF THE CONTINUOUS PROBLEM
In this section, we shall introduce a regularized version of the nonlinear PBE for both
analysis and discretization purposes. We first transfer the original equation posed on
the whole space to a truncated domain using an artificial boundary condition taken from
an approximate analytical solution. Then we use the fundamental solution in the whole
space to get rid of the singularities caused by δ distributions. We shall mainly focus on
more difficult problems in three dimensions. Formulation and results in two dimensions
aresimilarandrelativelyeasy.
6 L.CHEN,M.HOLST,ANDJ.XU
Let Ω ⊂ R3 with a convex and Lipschitz-continuous boundary ∂Ω, and Ω ⊂ Ω.
m
In the numerical simulation, for simplicity, we usually choose Ω to be a ball or cube
containing a molecule region. The solvent region is chosen as Ω ∩ Ω and will be still
s
denotedbyΩ . On∂Ωwechoosetheboundaryconditionu˜ = g,with
s
(cid:18) e2 (cid:19)(cid:88)Ni e−κ|x−xi|
g = c . (3.1)
k T ε |x−x |
B s i
i=1
The boundary condition is usually taken to be induced by a known analytical solution
to one of several possible simplifications of the linearized PBE. Far from the molecule,
such analytical solutions provide a highly accurate boundary condition approximation
for the general nonlinear PBE on a truncation of R3. For example, (3.1) arises from the
use of the Green’s function for the Helmholtz operator arising from linearizations of the
Poisson–Boltzmannoperator,whereasingleconstantglobaldielectricvalueofε isused
s
togeneratetheapproximateboundarycondition. (Thisisthecaseofarod-likemolecule
approximation; cf. [29].) Another approach to handling the boundary condition more
accuratelyistosolvethePBEwithboundaryconditionssuchas(3.1)onalargeΩ(with
a coarse mesh) and then solve it in a smaller Ω (with a fine mesh) with the boundary
condition provided by the earlier coarse mesh solution. The theoretical justification of
this approach can be found at [28] using the two-grid theory [53]. We are not going to
discussmoreonthechoiceoftheboundaryconditioninthispaper.
Employing(3.1)weobtainthenonlinearPBEonatruncateddomain:
(cid:88)Nm
−∇·(ε∇u˜)+κ¯2sinh(u˜) = q δ in Ω, (3.2)
i i
i=1
u˜ = g on ∂Ω. (3.3)
Thisis,inmostrespects,astandardboundary-valueproblemforanonlinearsecond-order
elliptic partial differential equation. However, the right side contains a linear combina-
tionofδdistributions,whichindividuallyandtogetherarenotinH−1(Ω);thuswecannot
apply standard techniques such as classical potential theory. This has at times been the
sourceofsomeconfusioninthemolecularmodelingcommunity,especiallywithrespect
to the design of convergent numerical methods. More precisely, we will see shortly that
thesolutiontothenonlinearPoisson–Boltzmannequationissimplynotgloballysmooth
enough to expect standard numerical methods (currently used by most PBE simulators)
to produce approximations that converge to the solution to the PBE in the limit of mesh
refinement.
In order to gain a better understanding of the properties of solutions to the nonlinear
PBE, primarily so that we can design new provably convergent numerical methods, we
shall propose a decomposition of the solution to separate out the singularity caused by
the δ distributions. This decomposition will turn out to be the key idea that will allow
us to design discretization techniques for the nonlinear PBE which have provably good
approximationpropertiesand,basedonthis,alsodesignanewtypeofadaptivealgorithm
whichisprovablyconvergentforthenonlinearPBE.
Wenowgivethisdecomposition. Itiswellknownthatthefunction
q 1
i
G =
i
ε |x−x |
m i
solvestheequation
−∇·(ε ∇G ) = q δ in R3.
m i i i
FINITEELEMENTAPPROXIMATIONOFTHENONLINEARPOISSON-BOLTZMANNEQUATION 7
We thus decompose the unknown u˜ as an unknown smooth function u and a known
singularfunctionG:
u˜ = u+G,
with
(cid:88)Nm
G = G . (3.4)
i
i=1
Substitutingthedecompositioninto(3.2),wethenobtain
−∇·(ε∇u)+κ¯2sinh(u+G) = ∇·((ε−ε )∇G) in Ω, (3.5)
m
u = g −G on ∂Ω, (3.6)
and call it the RPBE. The singularities of the δ distribtuions are transferred to G, which
then exhibits degenerate behavior at each {x } ⊂ Ω . At those points, both sinhG(x )
i m i
and∇G(x )exhibitblowup. However,sinceGisknownanalytically,oneavoidshaving
i
to build numerical approximations to G. Moreover, both of the coefficients κ¯ and ε −
ε are zero inside Ω , where the blowup behavior arises. Due to this cutoff nature of
m m
coefficients, the RPBE is a mathematically well defined nonlinear second-order elliptic
equation for the regularized solution u with the source term in H−1. We give a fairly
standard argument in the next section to show that it also admits a unique solution u ∈
H1, even though the original solution u˜ ∈/ H1 due to the singularities present in G. In
the remainder of the paper we shift our focus to establishing additional estimates and
developing an approximation theory to guide the design of convergent methods, both
nonadaptiveandadaptive.
Beforemovingon,itisusefultonotethat,awayfrom{x },thefunctionGissmooth.
i
In particular, we shall make use of the fact that G ∈ C∞(Ω ) ∩ C∞(Γ) ∩ C∞(∂Ω)
s
in the later analysis. Also, a key technical tool will be a further decomposition of the
regularizedsolutionuintolinearandnonlinearparts,u = ul +un,whereul satisfies
−∇·(ε∇ul) = ∇·((ε−ε )∇G) inΩ, (3.7)
m
ul = 0 on∂Ω, (3.8)
andwhereun satisfies
−∇·(ε∇un)+κ¯2sinh(un +ul +G) = 0 inΩ, (3.9)
un = g −G on∂Ω. (3.10)
4. EXISTENCE AND UNIQUENESS
In this section we shall discuss the existence and uniqueness of the solution of the
continuous RPBE. The arguments we use in this section appear essentially in [29], ex-
cept there the PBE was artificially regularized by replacing the delta distributions with
H−1-approximations directly rather than being regularized through a singular function
expansion.
We first write out the weak formulation. Since ∆G = 0 away from {x }, through
i
integrationbypartswegettheweakformulationofRPBE:Find
u ∈ M := {v ∈ H1(Ω)|ev,e−v ∈ L2(Ω ), and v = g −G on ∂Ω}
s
suchthat
A(u,v)+(B(u),v)+(cid:104)f ,v(cid:105) = 0 ∀v ∈ H1(Ω), (4.1)
G 0
where
• A(u,v) = (ε∇u,∇u),
8 L.CHEN,M.HOLST,ANDJ.XU
• (B(u),v) = (κ¯2sinh(u+G),v),and
(cid:82)
• (cid:104)f ,v(cid:105) = (ε−ε )∇G·∇v.
G Ω m
LetusdefinetheenergyonM:
(cid:90)
ε
E(w) = |∇w|2 +κ¯2cosh(w+G)+(cid:104)f ,w(cid:105).
G
2
Ω
Itiseasytocharacterizethesolutionof(4.1)astheminimizeroftheenergy.
Lemma4.1. Ifuisthesolutionoftheoptimizationproblem,i.e.,
E(u) = inf E(w),
w∈M
thenuisthesolutionof(4.1).
Proof. For any v ∈ H1(Ω) and any t ∈ R, the function F(t) = E(u + tv) attains the
0
minimalpointatt = 0,andthusF(cid:48)(0) = 0,whichgivesthedesiredresult. (cid:3)
We now recall some standard variational analysis on the existences of the minimizer.
InwhatfollowswesupposeS isasetinsomeBanachspaceV withnorm(cid:107)·(cid:107),andJ(u)
is a functional defined on S. S is called weakly sequential compact if, for any sequence
{u } ⊂ S, there exists a subsequence {u } such that u (cid:42) u ∈ S, where (cid:42) stands for
k ki ki
the convergence in the weak topology. For any u (cid:42) u, if J(u ) → J(u), we say J is
k k
weaklycontinuousatu;if
J(u) ≤ liminfJ(u ),
k
k→∞
we say J is weakly lower semicontinuous (w.l.s.c.) at u. The following theorem can be
provedbythedefinitioneasily.
Theorem4.2. If
1. S isweaklysequentialcompact,and
2. J isweaklylowersemicontinuousonS,
thenthereexistsu ∈ S suchthat
J(u) = inf J(w).
w∈S
We shall give conditions for the weakly sequential compactness and weakly lower
semicontinuity. First we use the fact that a bounded set in a reflexive Banach space is
weaklysequentialcompact.
Lemma4.3. Onehasthefollowingresults:
1. TheclosedunitballinareflexiveBanachspaceV isweaklysequentialcompact.
2. Iflim J(v) = ∞,then
(cid:107)v(cid:107)→∞
inf J(w) = inf J(w).
w∈V w∈S
The next lemma concerns when the functional is w.l.s.c. The proof can be found
at[57].
Lemma4.4. IfJ isaconvexfunctionalonaconvexsetS andJ isGaˆteauxdifferentiable,
thenJ isw.l.s.c.onS.
Now we are in the position to establish the existence and uniqueness of solutions to
theRPBE.
Theorem4.5. Thereexistsauniqueu ∈ M ⊂ H1(Ω)suchthat
E(u) = inf E(w).
w∈M
FINITEELEMENTAPPROXIMATIONOFTHENONLINEARPOISSON-BOLTZMANNEQUATION 9
Proof. ItiseasytoseeE(w)isdifferentiableinM with
(cid:104)DE(u),v(cid:105) = A(u,v)+(B(u),v)+(cid:104)f ,v(cid:105).
G
Toprovetheexistenceoftheminimizer,weneedonlytoverifythat
(1) M isaconvexset,
(2) E isconvexonM,and
(3) lim E(v) = ∞.
(cid:107)v(cid:107)1→∞
The verification of (1) is easy and thus skipped here. (2) comes from the convexity
of functions x2 and cosh(x). Indeed E is strictly convex. (3) is a consequence of the
inequality
E(v) ≥ C(ε,κ¯)(cid:107)v(cid:107)2 +C(G,g), (4.2)
1
whichcanbeprovedasfollowing. First,byYoung’sinequalitywehaveforanyδ > 0
1
(cid:104)f ,v(cid:105) ≤ ε (cid:107)∇G(cid:107) (cid:107)∇v(cid:107) ≤ (cid:107)∇G(cid:107)2 +δε2(cid:107)∇ (cid:107)2 .
G s Ωs Ωs δ Ωs s v Ωs
Sincecosh(x) ≥ 0,wehavethenE(v) ≥ C(ε,κ¯)(cid:107)∇v(cid:107)2−(1/δ)(cid:107)∇G(cid:107)2 ,wherewecan
Ωs
ensure C(ε,κ¯) > 0 if δ is chosen to be sufficiently small. Then using norm equivalence
onM,weget(4.2). Theuniquenessoftheminimizercomesfromthestrictconvexityof
E. (cid:3)
5. CONTINUOUS A PRIORI L∞-ESTIMATES
Inthissection,weshallderiveaprioriL∞-estimatesofthesolutionoftheRPBE.The
mainresultofthissectionisthefollowingtheorem.
Theorem5.1. LetubetheweaksolutionofRPBEinH1(Ω). ThenuisalsoinL∞(Ω).
Note that we cannot apply the analysis of [31, 32] directly to the RPBE, since the
right side f ∈ H−1(Ω) and does not lie in L∞(Ω) as required for use of these results.
G
We shall overcome this difficulty through further decomposition of u into linear and
nonlinearparts.
Let u = ul + un, where ul ∈ H1(Ω) satisfies the linear elliptic equation (the weak
0
formof(3.7)–(3.8))
A(ul,v)+(cid:104)f ,v(cid:105) = 0 ∀v ∈ H1(Ω) (5.1)
G 0
and where un ∈ M satisfies the nonlinear elliptic equation (the weak form of (3.9)–
(3.10))
A(un,v)+(B(un +ul),v) = 0 ∀v ∈ H1(Ω). (5.2)
0
Theorem 5.1 then follows from the estimates of ul and un in Lemmas 5.2 and 5.3; cf.
(5.3)and(5.4).
Lemma5.2. Letul betheweaksolutionof(5.1). Then
ul ∈ L∞(Ω). (5.3)
Proof. Since∆G = 0inΩ ,usingintegralbypartswecanrewritethefunctionalf as
s G
(cid:18) (cid:19)
∂G
(cid:104)f ,v(cid:105) = ((ε−ε )∇G,∇v) = [ε] ,v ,
G m
∂n
Γ Γ
where [ε] = ε −ε is the jump of ε at the interface. We shall still use f to denote the
s m G
smoothfunction[ε]∂G onΓ.
∂nΓ
10 L.CHEN,M.HOLST,ANDJ.XU
It is easy to see that the linear equation (5.1) is the weak formulation of the elliptic
interfaceproblem
(cid:20) ∂ul(cid:21)
−∇·(ε∇ul) = 0inΩ [ul] = 0, ε = f onΓ, and u = 0on∂Ω.
G
∂n
Since f ∈ C∞(Γ) and Γ ∈ C2, by the regularity result of the elliptic interface problem
G
[4,12,20,41],wehaveul ∈ H2(Ω )∩H2(Ω )∩H1(Ω). Inparticularbytheembedding
m s 0
theoremweconcludethatul ∈ L∞(Ω). (cid:3)
Toderiveasimilarestimateforthenonlinearpartun,wedefine
(cid:16) (cid:17) (cid:16) (cid:17)
α(cid:48) = argmax κ¯2sinh(c+ sup(ul +G)) ≤ 0 , α = min α(cid:48),inf(g −G) ,
c x∈Ωs ∂Ω
(cid:16) (cid:17) (cid:16) (cid:17)
β(cid:48) = argmin κ¯2sinh(c+ inf (ul +G)) ≥ 0 , β = max β(cid:48),sup(g −G) .
c x∈Ωs ∂Ω
ThenextlemmagivestheaprioriL∞-estimateofun.
Lemma5.3. Letun betheweaksolutionof(5.2). Thenα ≤ un ≤ β,andthus
un ∈ L∞(Ω). (5.4)
Proof. Weuseacutoff-functionargumentsimilartothatusedin[31]. Sincetheboundary
conditiong−G ∈ C∞(∂Ω),wecanfindau ∈ H1(Ω)suchthatu = g−Gon∂Ωin
D D
thetracesense,ormoreprecisely
Tu = g −G,
D
whereT : Ω (cid:55)→ ∂Ωisthetraceoperator. Thenthesolutioncanbewrittenun = u +u ,
D 0
with u ∈ H1(Ω). Let φ = (un − β)+ = max(un − β,0) and φ = (un − α)− =
0 0
min(un −α,0). Thenfrom
0 ≤ φ = (un −β)+ = (u +u −β)+ ≤ (u −β)+ +u+,
D 0 D 0
0 ≥ φ = (un −α)− = (u +u −α)− ≥ (u −α)− +u−,
D 0 D 0
and
0 ≤ Tφ ≤ T(u −β)+ +Tu+ = 0,
D 0
0 ≥ Tφ ≥ T(u −α)− +Tu− = 0,
D 0
weconcludethatbothφ,φ ∈ H1(Ω). Thusforeitherφ = φorφ = φ,wehave
0
(ε∇un,∇φ)+(κ¯2sinh(un +ul +G),φ) = 0.
Notethatφ ≥ 0inΩanditssupportisthesetY = {x ∈ Ω¯ |un(x) ≥ β}. OnY,wehave
(cid:16) (cid:17)
κ¯2sinh(un +ul +G) ≥ κ¯2sinh β(cid:48) + inf (ul +G) ≥ 0.
x∈Ωs
Similarly,φ ≤ 0inΩwithsupportsetY = {x ∈ Ω¯ |un(x) ≤ α}. OnY,wenowhave
(cid:16) (cid:17)
κ¯2sinh(un +ul +G) ≤ κ¯2sinh α(cid:48) + inf (ul +G) ≤ 0.
x∈Ωs
Togetherthisimplies
0 ≥ (ε∇un,∇φ) = (ε∇(un −β),∇φ) = ε(cid:107)∇φ(cid:107)2 ≥ 0
foreitherφ = φorφ = φ. UsingthePoincareinequalitywehavefinally
0 ≤ (cid:107)φ(cid:107) (cid:46) (cid:107)∇φ(cid:107) ≤ 0,