An Error Analysis for the Finite Element Approximation to the Steady-State Poisson-Nernst-Planck Equations

Poisson-Nernst-Planck equations are a coupled system of nonlinear partial differential equations consisting of the Nernst-Planck equation and the electrostatic Poisson equation with delta distribution sources, which describe the electrodiffusion of ions in a solvated biomolecular system. In this paper, some error bounds for a piecewise finite element approximation to this problem are derived. Several numerical examples including biomolecular problems are shown to support our analysis.     

Multiscale molecular dynamics using the matched interface and boundary method

The Poisson-Boltzmann (PB) equation is an established multiscale model for electrostatic analysis of biomolecules and other dielectric systems. PB based molecular dynamics (MD) approach has a potential to tackle large biological systems. Obstacles that hinder the current development of PB based MD methods are concerns in accuracy, stability, efficiency and reliability. The presence of complex solvent-solute interface, geometric singularities and charge singularities leads to challenges in the numerical solution of the PB equation and electrostatic force evaluation in PB based MD methods. Recently, the matched interface and boundary (MIB) method has been utilized to develop the first second order accurate PB solver that is numerically stable in dealing with discontinuous dielectric coefficients, complex geometric singularities and singular source charges. The present work develops the PB based MD approach using the MIB method. New formulation of electrostatic forces is derived to allow the use of sharp molecular surfaces. Accurate reaction field forces are obtained by directly differentiating the electrostatic potential. Dielectric boundary forces are evaluated at the solvent-solute interface using an accurate Cartesian-grid surface integration method. The electrostatic forces located at reentrant surfaces are appropriately assigned to related atoms. Extensive numerical tests are carried out to validate the accuracy and stability of the present electrostatic force calculation. The new PB based MD method is implemented in conjunction with the AMBER package. MIB based MD simulations of biomolecules are demonstrated via a few example systems.     

Second-order Poisson Nernst-Planck solver for ion channel transport

The Poisson Nernst-Planck (PNP) theory is a simplified continuum model for a wide variety of chemical, physical and biological applications. Its ability of providing quantitative explanation and increasingly qualitative predictions of experimental measurements has earned itself much recognition in the research community. Numerous computational algorithms have been constructed for the solution of the PNP equations. However, in the realistic ion-channel context, no second order convergent PNP algorithm has ever been reported in the literature, due to many numerical obstacles, including discontinuous coefficients, singular charges, geometric singularities, and nonlinear couplings. The present work introduces a number of numerical algorithms to overcome the abovementioned numerical challenges and constructs the first second-order convergent PNP solver in the ion-channel context. First, a Dirichlet to Neumann mapping (DNM) algorithm is designed to alleviate the charge singularity due to the protein structure. Additionally, the matched interface and boundary (MIB) method is reformulated for solving the PNP equations. The MIB method systematically enforces the interface jump conditions and achieves the second order accuracy in the presence of complex geometry and geometric singularities of molecular surfaces. Moreover, two iterative schemes are utilized to deal with the coupled nonlinear equations. Furthermore, extensive and rigorous numerical validations are carried out over a number of geometries, including a sphere, two proteins and an ion channel, to examine the numerical accuracy and convergence order of the present numerical algorithms. Finally, application is considered to a real transmembrane protein, the Gramicidin A channel protein. The performance of the proposed numerical techniques is tested against a number of factors, including mesh sizes, diffusion coefficient profiles, iterative schemes, ion concentrations, and applied voltages. Numerical predictions are compared with experimental measurements.     

Diffusion front capturing schemes for a class of Fokker–Planck equations: Application to the relativistic heat equation

In this research work we introduce and analyze an explicit conservative finite difference scheme to approximate the solution of initial-boundary value problems for a class of limited diffusion Fokker–Planck equations under homogeneous Neumann boundary conditions. We show stability and positivity preserving property under a Courant–Friedrichs–Lewy parabolic time step restriction. We focus on the relativistic heat equation as a model problem of the mentioned limited diffusion Fokker–Planck equations. We analyze its dynamics and observe the presence of a singular flux and an implicit combination of nonlinear effects that include anisotropic diffusion and hyperbolic transport. We present numerical approximations of the solution of the relativistic heat equation for a set of examples in one and two dimensions including continuous initial data that develops jump discontinuities in finite time. We perform the numerical experiments through a class of explicit high order accurate conservative and stable numerical schemes and a semi-implicit nonlinear Crank–Nicolson type scheme.     

Electrostatic potential of point charges inside dielectric prolate spheroids

The exact solution to an electrostatic problem of finding the electric potential of point charges inside a dielectric prolate spheroid is discussed in this note by using the classical electrostatic theory, where the prolate spheroid is embedded in a dissimilar dielectric medium. Such a problem may find its application in hybrid solvent biomolecular simulations, in which biomolecules and a part of solvent molecules within a dielectric cavity are explicitly modeled while a surrounding dielectric continuum is used to model bulk effects of the solvent beyond the cavity. Numerical experiments have demonstrated the convergence of the proposed series solutions.     

Navier-Stokes方程的最小二乘等几何方法

基于Hermite多项式的C¹型单元构造复杂,限制了最小二乘有限元法的应用.引入高阶光滑的非均匀有理B样条作为基函数简化C¹型单元构造,提出求解黏性不可压流动Navier-Stokes方程的最小二乘等几何方法.用Newton法或Picard法对Navier-Stokes方程线性化,用线性化偏微分方程的余量定义最小二乘泛函,导出最小二乘变分方程,用NURBS构造高阶光滑的有限维空间来近似速度场和压力场.计算表明:本文方法计算的二维顶盖驱动流数值解能准确描述流动状况,计算的二维通道内圆柱绕流全局质量损失由最小二乘有限元法的6%降为0.018%,该方法可用于Navier-Stokes方程的求解,并且具有较好的质量守恒性.     

二维Helmholtz边界超奇异积分方程解析研究

基于常规边界元法及超奇异边界积分方程复线性耦合的Burton-Miller方法应用于无限域声学问题的最大难点在于处理超奇异积分（二维问题）.目前,此类超奇异积分主要使用各种弱奇异/正则化方法求解,而这些弱奇异/正则化方法具有时间消耗大等弱点.基于围道积分定理,本文给出一种使用常值单元的二维Helmholtz边界超奇异积分的解析表达式.在有限部分积分意义下,所有的奇异和超奇异积分可以解析表达.数值算例表明该解析表达式是有效的.     

基于局部基本解法的静电场仿真分析

将局部基本解方法应用于静电场问题的模拟与分析。局部基本解方法是利用控制方程的基本解,基于局部理论和移动最小二乘原理提出的一种无网格算法。相比于有限元和有限差分等传统网格类方法,该方法仅需离散节点,避免了复杂的网格剖分难题。作为一种半解析数值技术,物理问题的基本解被作为插值基函数建立数值离散模型,从而保证了算法的较高精度。此外,与具有全局离散格式的无网格方法相比,局部基本解法更适用于高维复杂几何和大尺度模拟。二维和三维数值试验表明,该方法具有实施方便灵活,计算精度高和计算速度快等优势。为静电场仿真研究开辟新的途径,拓展了局部基本解方法的应用领域。     

Solution of the Electrodiffusion Equations in the Quasi-Stationary Approximation

The electrodiffusion equations are solved together with the equation for the current running through the boundary of the examined system with the environment. The complete system of equations allows the volt-ampere characteristic of the system to be obtained for an arbitrary time dependence of the external potential. The procedure of finding the first-order approximation is developed to solve the system of electrodiffusion equations.     

Discrete Charges on a Two Dimensional Conductor

We investigate the electrostatic equilibria of N discrete charges of size 1/N on a two dimensional conductor (domain). We study the distribution of the charges on symmetric domains including the ellipse, the hypotrochoid and various regular polygons, with an emphasis on understanding the distributions of the charges, as the shape of the underlying conductor becomes singular. We find that there are two regimes of behavior, a symmetric regime for smooth conductors, and a symmetry broken regime for “singular” domains. For smooth conductors, the locations of the charges can be determined, to within $O\left( {\sqrt {\log {N \mathord{\left/ {\vphantom {N {N^2 }}} \right. \kern-0em} {N^2 }}} } \right)$ by an integral equation due to Pommerenke [ Math. Ann., 179: 212–218, (1969)]. We present a derivation of a related (but different) integral equation, which has the same solutions. We also solve the equation to obtain (asymptotic) solutions which show universal behavior in the distribution of the charges in conductors with somewhat smooth cusps. Conductors with sharp cusps and singularities show qualitatively different behavior, where the symmetry of the problem is broken, and the distribution of the discrete charges does not respect the symmetry of the underlying domain. We investigate the symmetry breaking both theoretically, and numerically, and find good agreement between our theory and the numerics. We also find that the universality in the distribution of the charges near the cusps persists in the symmetry broken regime, although this distribution is very different from the one given by the integral equation.     

半解析对偶棱边元及其在波导不连续性问题中的应用

给出电磁波导的对偶变量变分原理,并采用对偶棱边元对波导的横截面进行半解析离散. 将波导中沿纵向均匀的区段视为子结构,运用基于Riccati方程的精细积分算法求出其出口刚度阵,然后与不均匀区段的常规有限元网格拼装即可对波导不连续性问题进行求解. 半解析对偶棱边元的采用可以在最大程度上对有限元网格进行缩减,并且能够在不增加计算量的前提下任意增加子结构的长度,从而可以将截断求解区域的人工边界设置在距离不均匀区段充分远的地方,极大地减少了近似边界条件所带来的误差. 数值算例证明这种方法具有很高的精度与效率. 关键词： 波导的不连续性 半解析辛分析 对偶棱边元 精细积分     

Homogenization of 3D finite photonic crystals with heterogeneous permittivity and permeability

We consider a heterogeneous magneto-dielectric photonic crystal and derive the so-called ‘homogenized Maxwell system’ via the multi-scale method and provide ad hoc proofs for the convergence of the electromagnetic field towards the homogeneous one using the notion of two-scale convergence. The homogenized medium is described by anisotropic matrices of permittivity and permeability, deduced from the resolution of two annex problems of electrostatic type on a periodic cell. Noteworthily, this asymptotic analysis also covers the case of photonic crystals with non-cuboidal periodic cells. We solve numerically the associated system of partial differential equations with a method of fictitious charges and a finite element method (FEM) in order to exhibit the matrices of effective permittivity and permeability for given magneto-dielectric periodic composites. We then compare our results in the 2D case against some Fourier expansion approach and provide duality relations in the case of magneto-dielectric checkerboards. We further compute some low-frequency eigenmodes of a photonic crystal fiber with metallic outer boundary and compare them with the eigenmodes of a corresponding effective anisotropic waveguide, thanks to the FEM. Finally, we derive the effective properties of a 3D photonic crystal both through classical homogenization (solving numerically two decoupled annex problems) and Bloch wave homogenization. In the case of spherical inclusions, the latter approach amounts to evaluating the slope of the first band around the origin on a Bloch diagram which we compute using finite edge elements.     

Extended Nernst–Planck Equation Incorporating Partial Dehydration Effect

Novel ionic transporting phenomena emerge as nanostructures approach the molecular scale.At the sub-2 nm scale,widely used continuum equations,such as the Nernst-Planck equation,break down.Here,we extend the Nernst-Planck equation by adding a partial dehydration effect.Our model agrees with the reported ion fluxes through graphene oxide laminates with sub-2 nm interlayer spacing,outperforming previous models.We also predict that the selectivity sequences of alkali metal ions depend on the geometries of the nanostructures.Our model opens a new avenue for the investigation of the underlying mechanisms in nanofluidics at the sub-2 nm scale.     

On the Theory of Electric-Current Flow through a Mixture of Charged Particles

Current conduction through a mixture made of two species of positively charged particles is considered where one of the latter species participates in the exchange with the surrounding medium. A solution to the electrodiffusion equations together with Poisson's equation is obtained in the first approximation in terms of the small parameter. A condition is determined where the distribution of charged particles involved in the exchange with the surrounding medium is derived using the diffusion equation for neutral particles. It is shown that the solution to the electrodiffusion equations contains a component decaying with time.     

A sweeping preconditioner for time-harmonic Maxwell’s equations with finite elements

This paper is concerned with preconditioning the stiffness matrix resulting from finite element discretizations of Maxwell’s equations in the high frequency regime. The moving PML sweeping preconditioner, first introduced for the Helmholtz equation on a Cartesian finite difference grid, is generalized to an unstructured mesh with finite elements. The method dramatically reduces the number of GMRES iterations necessary for convergence, resulting in an almost linear complexity solver. Numerical examples including electromagnetic cloaking simulations are presented to demonstrate the efficiency of the proposed method.     

A nonconforming combination for solving elliptic problems with interfaces

Nonconforming combinations are provided for solving interface problems of elliptic equations. In these approaches, the Ritz-Galerkin method with particular solutions is used for the part of a solution domain where there are interface singular points; and the conventional finite element method is used for the rest of the solution domain. In addition, admissible functions chosen are constrained to be continuous only at the element nodes on the common boundary of the subdomains. Error bounds are derived in the Sobolev norms, and numerical experiments are given for solving a model interface problem of the equation, −Δu + U = 0. Moreover, a significant coupling relation, L + 1 = O(|ln h|), is found for interface problems by using the nonconforming combinations, where (L + 1) is the total number of particular solutions used in the Ritz-Galerkin method, and h is the maximal boundary length of triangular elements in the finite element method.     

A method of cauchy integral equation for non-coherent transfer in half-space

A technique is presented which allows easy construction of solutions for various half-space problems arising in non-coherent radiative transfer with complete redistribution. By use of an inverse Laplace transform method, Wiener-Hopf integral equations are reduced to Cauchy-type singular integral equations. The factorization technique used by Case and Zweifel for coherent scattering can then be carried over to non-coherent transfer. The method is applied to the inhomogeneous integral equation for the source function of a two-level atom, previously solved by Ivanov. It is also applied to the conservative, homogeneous case and to singular Wiener-Hopf equations arising from asymptotic expansions in the limit of vanishing probability of collisional destruction ?. Consequences for the scaling laws in a finite slab are examined in a companion paper.     

Fast finite difference solvers for singular solutions of the elliptic Monge–Ampère equation

The elliptic Monge–Ampère equation is a fully nonlinear Partial Differential Equation which originated in geometric surface theory, and has been applied in dynamic meteorology, elasticity, geometric optics, image processing and image registration. Solutions can be singular, in which case standard numerical approaches fail.In this article we build a finite difference solver for the Monge–Ampère equation, which converges even for singular solutions. Regularity results are used to select a priori between a stable, provably convergent monotone discretization and an accurate finite difference discretization in different regions of the computational domain. This allows singular solutions to be computed using a stable method, and regular solutions to be computed more accurately. The resulting nonlinear equations are then solved by Newton’s method.Computational results in two and three-dimensions validate the claims of accuracy and solution speed. A computational example is presented which demonstrates the necessity of the use of the monotone scheme near singularities.     

The HN method for solving linear transport equation: theory and applications

The system of singular integral equations which is obtained from the integro-differential form of the linear transport equation using the Placzek lemma is solved. The exit distributions at the boundaries of the various media and the infinite medium Green's function are used. The process is applied to the half-space and finite slab problems. The neutron angular density in terms of singular eigenfunctions of the method of elementary solutions is also used to derive the same analytical expressions.     

球几何中子输运保正线性间断有限元格式

针对球几何中子输运方程线性间断有限元方法计算的负中子通量问题,构造了保正线性间断有限元格式,该格式保持中子角通量0阶矩和1阶矩。现有方法计算中子角通量非负时,采用传统的线性间断有限元方法,求解线性方程组;原方法计算出现负通量,则采用构造的保正格式,求解非线性方程组。编制了球几何中子输运问题保正格式程序模块,并集成到应用程序。数值算例表明构造的保正格式计算的中子通量非负,有效降低数值误差,提高数值计算的精度。