Abstract
A second order accurate (in time) numerical scheme is proposed and analyzed for the Poisson–Nernst–Planck equation (PNP) system, reformulated as a non-constant mobility \(H^{-1}\) gradient flow in the Energetic Variational Approach (EnVarA). The centered finite difference is taken as the spatial discretization. Meanwhile, the highly nonlinear and singular nature of the logarithmic energy potentials has always been the essential difficulty to design a second order accurate scheme in time, while preserving the variational energetic structures. The mobility function is updated with a second order accurate extrapolation formula, for the sake of unique solvability. A modified Crank–Nicolson scheme is used to approximate the logarithmic term, so that its inner product with the discrete temporal derivative exactly gives the corresponding nonlinear energy difference; henceforth the energy stability is ensured for the logarithmic part. In addition, nonlinear artificial regularization terms are added in the numerical scheme, so that the positivity-preserving property could be theoretically proved, with the help of the singularity associated with the logarithmic function. Furthermore, an optimal rate convergence analysis is provided in this paper, in which the higher order asymptotic expansion for the numerical solution, the rough error estimate and refined error estimate techniques have to be included to accomplish such an analysis. This work combines the following theoretical properties for a second order accurate numerical scheme for the PNP system: (i) second order accuracy in both time and space, (ii) unique solvability and positivity, (iii) energy stability, and (iv) optimal rate convergence. A few numerical results are also presented.
Similar content being viewed by others
Data Availibility
Data Availability Enquiries about data availability should be directed to the authors.
References
Abels, H., Wilke, M.: Convergence to equilibrium for the Cahn–Hilliard equation with a logarithmic free energy. Nonlinear Anal. 67, 3176–3193 (2007)
Baskaran, A., Hu, Z., Lowengrub, J., Wang, C., Wise, S.M., Zhou, P.: Energy stable and efficient finite-difference nonlinear multigrid schemes for the modified phase field crystal equation. J. Comput. Phys. 250, 270–292 (2013)
Baskaran, A., Lowengrub, J., Wang, C., Wise, S.: Convergence analysis of a second order convex splitting scheme for the modified phase field crystal equation. SIAM J. Numer. Anal. 51, 2851–2873 (2013)
Bazant, M.Z., Thornton, K., Ajdari, A.: Diffuse-charge dynamics in electrochemical systems. Phys. Rev. E 70(2), 021506 (2004)
Ben, Y., Chang, H.C.: Nonlinear Smoluchowski slip velocity and micro-vortex generation. J. Fluid Mech. 461, 229–238 (2002)
Chen, W., Wang, C., Wang, X., Wise, S.M.: A linear iteration algorithm for energy stable second order scheme for a thin film model without slope selection. J. Sci. Comput. 59, 574–601 (2014)
Chen, W., Feng, W., Liu, Y., Wang, C., Wise, S.M.: A second order energy stable scheme for the Cahn–Hilliard–Hele–Shaw equation. Discrete Contin. Dyn. Syst. Ser. B 24(1), 149–182 (2019)
Chen, W., Wang, C., Wang, X., Wise, S.M.: Positivity-preserving, energy stable numerical schemes for the Cahn–Hilliard equation with logarithmic potential. J. Comput. Phys.: X 3, 100031 (2019)
Cheng, K., Wang, C., Wise, S.M., Yue, X.: A second-order, weakly energy-stable pseudo-spectral scheme for the Cahn–Hilliard equation and its solution by the homogeneous linear iteration method. J. Sci. Comput. 69, 1083–1114 (2016)
Debussche, A., Dettori, L.: On the Cahn–Hilliard equation with a logarithmic free energy. Nonlinear Anal. 24, 1491–1514 (1995)
Diegel, A., Wang, C., Wise, S.M.: Stability and convergence of a second order mixed finite element method for the Cahn–Hilliard equation. IMA J. Numer. Anal. 36, 1867–1897 (2016)
Diegel, A., Wang, C., Wang, X., Wise, S.M.: Convergence analysis and error estimates for a second order accurate finite element method for the Cahn–Hilliard–Navier–Stokes system. Numer. Math. 137, 495–534 (2017)
Ding, J., Wang, C., Zhou, S.: Optimal rate convergence analysis of a second order numerical scheme for the Poisson–Nernst–Planck system. Numer. Math. Theor. Meth. Appl. 12, 607–626 (2019)
Dong, L., Feng, W., Wang, C., Wise, S.M., Zhang, Z.: Convergence analysis and numerical implementation of a second order numerical scheme for the three-dimensional phase field crystal equation. Comput. Math. Appl. 75(6), 1912–1928 (2018)
Dong, L., Wang, C., Zhang, H., Zhang, Z.: A positivity-preserving, energy stable and convergent numerical scheme for the Cahn–Hilliard equation with a Flory–Huggins–deGennes energy. Commun. Math. Sci. 17, 921–939 (2019)
Dong, L., Wang, C., Zhang, H., Zhang, Z.: A positivity-preserving second-order BDF scheme for the Cahn–Hilliard equation with variable interfacial parameters. Commun. Comput. Phys. 28, 967–998 (2020)
Duan, C., Liu, C., Wang, C., Yue, X.: Convergence analysis of a numerical scheme for the porous medium equation by an energetic variational approach. Numer. Math. Theor. Methods Appl. 13, 1–18 (2020)
Eisenberg, R.S.: Computing the field in proteins and channels. J. Mem. Biol. 150, 1–25 (1996)
Eisenberg, B., Hyon, Y., Liu, C.: 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)
Elliott, C.M., Garcke, H.: On the Cahn–Hilliard equation with degenerate mobility. SIAM J. Math. Anal. 27, 404–423 (1996)
Flavell, A., Machen, M., Eisenberg, R., Kabre, J., Liu, C., Li, X.: A conservative finite difference scheme for Poisson–Nernst–Planck equations. J. Comput. Electron. 13, 235–249 (2014)
Gavish, N., Yochelis, A.: Theory of phase separation and polarization for pure ionic liquids. J. Phys. Chem. Lett. 7, 1121–1126 (2016)
Guan, Z., Wang, C., Wise, S.M.: A convergent convex splitting scheme for the periodic nonlocal Cahn–Hilliard equation. Numer. Math. 128, 377–406 (2014)
Guan, Z., Lowengrub, J.S., Wang, C.: Convergence analysis for second order accurate schemes for the periodic nonlocal Allen–Cahn and Cahn–Hilliard equations. Math. Methods Appl. Sci. 40(18), 6836–6863 (2017)
Guo, J., Wang, C., Wise, S.M., Yue, X.: An \(H^2\) convergence of a second-order convex-splitting, finite difference scheme for the three-dimensional Cahn–Hilliard equation. Commun. Math. Sci. 14, 489–515 (2016)
Han, D., Wang, X.: A second order in time, uniquely solvable, unconditionally stable numerical scheme for Cahn–Hilliard–Navier–Stokes equation. J. Comput. Phys. 290, 139–156 (2015)
Hu, J., Huang, X.: A fully discrete positivity-preserving and energy-dissipative finite difference scheme for Poisson–Nernst–Planck equations. Numer. Math. 145, 77–115 (2020)
Hu, Z., Wise, S.M., Wang, C., Lowengrub, J.S.: Stable and efficient finite-difference nonlinear-multigrid schemes for the phase-field crystal equation. J. Comput. Phys. 228, 5323–5339 (2009)
Hunter, R.J.: Foundations of Colloid Science. Oxford University Press, Oxford (2001)
Jerome, J.W.: Analysis of Charge Transport. Mathematical Theory and Approximation of Semi-conductor Models. Springer, New York (1995)
Li, X., Qiao, Z., Wang, C.: Convergence analysis for a stabilized linear semi-implicit numerical scheme for the nonlocal Cahn–Hilliard equation. Math. Comput. 90, 171–188 (2021)
Liu, C., Wang, C., Wise, S., Yue, X., Zhou, S.: A positivity-preserving, energy stable and convergent numerical scheme for the Poisson–Nernst–Planck system. Math. Comput. 90(331), 2071–2106 (2021)
Lyklema, J.: Fundamentals of Interface and Colloid Science. Volume II Solid–Liquid Interfaces. Academic Press Limited, San Diego (1995)
Markowich, P.A.: The Stationary Seminconductor Device Equations. Springer, Vienna (1986)
Markowich, P.A., Ringhofer, C.A., Schmeiser, C.: Seminconductor Equations. Springer, New York (1990)
Miranville, A., Zelik, S.: Robust exponential attractors for Cahn–Hilliard type equations with singular potentials. Math. Methods Appl. Sci. 27, 545–582 (2004)
Nazarov, I., Promislow, K.: The impact of membrane constraint on PEM fuel cell water management. J. Electrochem. Soc. 154(7), 623–630 (2007)
Nonner, W., Chen, D.P., Eisenberg, B.: Progress and prospects in permeation. J. Gen. Physiol. 113, 773–782 (1999)
Promislow, K., Stockie, J.M.: Adiabatic relaxation of convective–diffusive gas transport in a porous fuel cell electrode. SIAM J. Appl. Math. 62(1), 180–205 (2001)
Samelson, R., Temam, R., Wang, C., Wang, S.: Surface pressure Poisson equation formulation of the primitive equations: numerical schemes. SIAM J. Numer. Anal. 41, 1163–1194 (2003)
Samelson, R., Temam, R., Wang, C., Wang, S.: A fourth order numerical method for the planetary geostrophic equations with inviscid geostrophic balance. Numer. Math. 107, 669–705 (2007)
Shen, J., Xu, J.: Unconditionally positivity preserving and energy dissipative schemes for Poisson–Nernst–Planck equations. Numer. Math. 148, 671–697 (2021)
Shen, J., Wang, C., Wang, X., Wise, S.M.: Second-order convex splitting schemes for gradient flows with Ehrlich–Schwoebel type energy: application to thin film epitaxy. SIAM J. Numer. Anal. 50, 105–125 (2012)
Wang, C., Liu, J.-G.: Convergence of gauge method for incompressible flow. Math. Comp. 69, 1385–1407 (2000)
Wang, C., Liu, J.-G.: Analysis of finite difference schemes for unsteady Navier–Stokes equations in vorticity formulation. Numer. Math. 91, 543–576 (2002)
Wang, C., Wise, S.M.: An energy stable and convergent finite-difference scheme for the modified phase field crystal equation. SIAM J. Numer. Anal. 49, 945–969 (2011)
Wang, C., Liu, J.-G., Johnston, H.: Analysis of a fourth order finite difference method for incompressible Boussinesq equations. Numer. Math. 97, 555–594 (2004)
Wang, L., Chen, W., Wang, C.: An energy-conserving second order numerical scheme for nonlinear hyperbolic equation with an exponential nonlinear term. J. Comput. Appl. Math. 280, 347–366 (2015)
Weinan, E., Liu, J.-G.: Projection method I. Convergence and numerical boundary layers. SIAM J. Numer. Anal. 32, 1017–1057 (1995)
Weinan, E., Liu, J.-G.: Projection method III. Spatial discretization on the staggered grid. Math. Comput. 71, 27–47 (2002)
Wise, S.M.: Unconditionally stable finite difference, nonlinear multigrid simulation of the Cahn–Hilliard–Hele–Shaw system of equations. J. Sci. Comput. 44, 38–68 (2010)
Wise, S.M., Wang, C., Lowengrub, J.: An energy stable and convergent finite-difference scheme for the phase field crystal equation. SIAM J. Numer. Anal. 47, 2269–2288 (2009)
Funding
This work is supported in part by the National Science Foundation (USA) grants NSF DMS-1759535, NSF DMS-1759536 (C. Liu), NSF DMS-2012669 and DMS-2309548 (C. Wang), NSF DMS-1719854, DMS-2012634 (S. Wise), National Natural Science Foundation of China 11971342 (X. Yue), and 12171319 (S. Zhou).
Author information
Authors and Affiliations
Corresponding author
Ethics declarations
Conflict of interest
The authors declare that they have no conflict of interest.
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Appendices
Appendix A: Proof of Proposition 4.1
In terms of the temporal discretization, the following local truncation error can be derived by a Taylor expansion in time, combined with the projection estimate (4.3):
with the projection accuracy order \(m_0 \ge 4\). In fact, the spatial functions \(G_n^{(0)}\) \(G_p^{(0)}\) are smooth enough in the sense that their derivatives are bounded.
Subsequently, the leading order temporal correction function \(({{\textsf{N}}}_{{\Delta t}, 1}, {{\textsf{P}}}_{{\Delta t}, 1})\) turns out to be the solution of the following linear equations:
In fact, existence of a solution of the above linear and parabolic PDE system is straightforward. It depends only on the projection solution \(({{\textsf{N}}}_N, {{\textsf{P}}}_N)\). And also, the derivatives of \(({{\textsf{N}}}_{{\Delta t}, 1}, {{\textsf{P}}}_{{\Delta t}, 1})\) in various orders are bounded. In turn, an application of the semi-implicit discretization to (A.3)–(A.4) gives
A combination of (A.1)–(A.2) and (A.5)–(A.6) results in the third order temporal truncation error for \(\hat{{\textsf{N}}}_1:= {\mathsf N}_N + {\Delta t}^2 {{{\mathcal {P}}}}_N {{\textsf{N}}}_{{\Delta t}, 1}\), \(\hat{{\textsf{P}}}_1:= {{\textsf{P}}}_N + {\Delta t}^2 {{{\mathcal {P}}}}_N {{\textsf{P}}}_{{\Delta t}, 1}\):
In the derivation of (A.7)–(A.8), the following linearized expansions have been applied:
in which property (3) of \(G_a^1 (x)\) (as stated in Lemma 3.1) is recalled. The corresponding expansions for \(G_{\hat{{\textsf{P}}}_1^m}^1 ( \hat{{\textsf{P}}}_1^{m+1} )\) could be similarly derived, and the technical details are skipped for the sake of brevity.
Similarly, the next order temporal correction function \(\big ({\mathsf N}_{{\Delta t}, 2}, {{\textsf{P}}}_{{\Delta t}, 2}\big )\) turns out to be the solution of following linear equations:
Again, the solution depends only on the exact solution \(({\mathsf N}, {{\textsf{P}}})\), with derivatives of various orders stay bounded. Of course, an application of the semi-implicit discretization to (A.10)–(A.11) gives
As a result, a combination of (A.10)–(A.11) and (A.12)–(A.13) yields the fourth order temporal truncation error for \(\hat{{\textsf{N}}}_2:= \hat{{\textsf{N}}}_1 + {\Delta t}^3 {{{\mathcal {P}}}}_N {{\textsf{N}}}_{{\Delta t}, 2}\), \(\hat{{\textsf{P}}}_2:= \hat{{\textsf{P}}}_1 + {\Delta t}^3 {{{\mathcal {P}}}}_N {\mathsf P}_{{\Delta t}, 2}\):
in which similar linearized expansions (as in (A.9)) have been used in the derivation.
In terms of spatial discretization, we construct the spatial correction term \(({{\textsf{N}}}_{h, 1}, {{\textsf{P}}}_{h, 1})\) to improve the spatial accuracy order. The following truncation error estimate for the spatial discretization is available, by using a straightforward Taylor expansion for the constructed profile \((\hat{{\textsf{N}}}_2, \hat{{\textsf{P}}}_2)\):
in which the average operator is taken in a similar form as (2.9). Similarly, the spatially discrete functions \(H_n^{(0)}\), \(H_p^{(0)}\) are smooth enough in the sense that their discrete derivatives are bounded. Because of the symmetry in the centered finite difference approximation, there is no \(O (h^3)\) truncation error term. In turn, the spatial correction function \(({{\textsf{N}}}_{h, 1}, {{\textsf{P}}}_{h, 1})\) is determined by solving the solution of the following linear PDE system:
Again, the solution depends only on the exact solution \(({\mathsf N}, {{\textsf{P}}})\), with the divided differences of various orders stay bounded. An application of a full discretization to (A.18)–(A.19) leads to
Finally, a combination of (A.18)–(A.19) and (A.20)–(A.21) yields the higher order truncation error for \((\hat{{\textsf{N}}}, \hat{\mathsf P})\), as given by (4.10)–(4.11). Of course, the linear expansions have been extensively utilized.
Moreover, we see that trivial initial data \({{\textsf{N}}}_{{\Delta t}, j} (\, \cdot \,, t=0), {{\textsf{P}}}_{{\Delta t}, j} (\, \cdot \,, t=0) \equiv 0\) could be taken (\(j=1, 2\)) as in (A.3)–(A.4), (A.10)–(A.11), respectively, as well as \(({{\textsf{N}}}_{h, 1}, {{\textsf{P}}}_{h, 1})\) in (A.18)–(A.19). Consequently, using similar arguments as in (4.4)–(4.5), we arrive at the mass conservative identities (4.12), (4.13). Notice that the first step of (4.13) is based on the fact that \(\hat{{\textsf{N}}} \in {{{\mathcal {B}}}}^K\), and the second step comes from the mass conservative property of \(\hat{{\textsf{N}}}\) at the continuous level. And also, the mass conservative property of \((\hat{{\textsf{N}}}, \hat{{\textsf{P}}})\) is stated in (4.13), and we conclude that the local truncation error \(\tau _n\), \(\tau _p\) has a similar property, so that (4.14) is proved.
Based on the fact that the temporal and spatial correction functions \(({{\textsf{N}}}_{{\Delta t}, j}, {{\textsf{P}}}_{{\Delta t}, j})\), \(({{\textsf{N}}}_{h, 1}, {{\textsf{P}}}_{h, 1})\) are bounded, we recall the separation property (4.2) for the exact solution. In turn, a similar property (4.15) becomes available for the constructed profile \((\hat{{\textsf{N}}}, \hat{{\textsf{P}}})\), in which the projection estimate (4.3) has been recalled. Of course, \({\Delta t}\) and h could be taken sufficiently small so that (4.15) is valid for a modified value \(\epsilon _0^\star \), such as \(\epsilon _0^\star = \frac{1}{4} \epsilon _0\).
Furthermore, we recall the fact that the correction functions stay bounded, in terms of both the spatial and temporal derivatives, since they only depend on \(({{\textsf{N}}}_N, {{\textsf{P}}}_N)\) and the exact solution. Therefore, a discrete \(W^{1,\infty }\) bound for \((\hat{{\textsf{N}}}, \hat{{\textsf{P}}})\) could be derived as in (4.16), as well as a bound (4.17) for its discrete temporal derivative. This completes the proof of Proposition 4.1.
Appendix B: Proof of Lemma 4.1
For the nonlinear inner product associated with the artificial regularization, the following fact is observed:
in which the positivity-preserving property of \(n^{m+1}\) and \(\check{{\textsf{N}}}^{m+1}\) has been applied. For the additional term in the artificial regularization part, we have the following estimates:
For the error term \(\langle {\tilde{n}}^{m+1}, G_{\hat{\mathsf N}^m}^1 ( \hat{{\textsf{N}}}^{m+1} ) - G_{n^m}^1 ( n^{m+1} ) \rangle \), we begin with the following decomposition:
in which \(G_a^1 (x)\) has been introduced in (3.1). For the \(\mathcal{NLE}_{2}^m\) error term, we apply the intermediate value theorem and obtain
Meanwhile, by property (3) in Lemma 3.1, we conclude that
By the combined fact that, \(\xi _n^{(m)}\) is between \(\hat{\mathsf N}^{m+1}\) and \(\eta _n^{(m)}\), \(\eta \) is between \(\hat{{\textsf{N}}}^m\) and \(n^m\), the following bound is available
in which the phase separation properties (4.15), (4.30) have been recalled. Then we arrive at
Based on this point-wise bound, the following inequality is available
For the \(\mathcal{NLE}_{1}^m\) error term, a similar nonlinear analysis could be performed:
This in turn leads to
Subsequently, a combination of (B.10) and (B.12) yields
Moreover, its combination with (B.1)–(B.4) indicates that
At each fixed grid point (i, j, k), if (i, j, k) is not in \(\Lambda _n\), i.e., \(0< n^{m+1}_{i,j,k} < 2C^* +1\), the following estimates are available:
In turn, the following inequality is valid for \(0< n^{m+1}_{i,j,k} < 2C^* +1\):
At each fixed grid point (i, j, k), if \((i, j, k) \in \Lambda _n\), i.e., \(n^{m+1}_{i,j,k} \ge 2C^* +1\), we see that \(n^{m+1}_{i,j,k} > \max ( \hat{{\textsf{N}}}^{m+1}_{i,j,k}, n^m_{i,j,k} )\), so that
Meanwhile, since \(n^{m+1}_{i,j,k} > \hat{{\textsf{N}}}^{m+1}_{i,j,k}\), we see that \({\tilde{n}}^{m+1}_{i,j,k} = \hat{{\textsf{N}}}^{m+1}_{i,j,k} - n^{m+1}_{i,j,k} < 0\), and the following fact is observed:
Subsequently, the following inequality could be derived, for \(n^{m+1}_{i,j,k} \ge 2 C^* +1\):
in which the \(\Vert \cdot \Vert _\infty \) estimate (4.28) has been applied in the second step, and the fact the \(\frac{\frac{1}{2} (C^* + 1) }{2 C^* +1} - \frac{3}{2} ( \epsilon _0^*)^{-1} {\Delta t}^2 \ge \frac{1}{4}\) has been used in the fourth step.
Consequently, a substitution of the point-wise inequalities (B.16), (B.19) into (B.14) results in the desired estimate (4.33), by taking \({\tilde{C}}_4 = \frac{9 (2 C^* +1)}{16} ( \epsilon _0^*)^{-2}\). The other inequality in (4.33) could be derived in the same manner; the technical details are skipped for the sake of brevity.
Finally, if \(K_n^* =0\) and \(K_p^*=0\), i.e, both \(\Lambda _n\) and \(\Lambda _p\) are empty sets, we see that inequality (B.16) is valid at every grid points. In turn, a summation in space results in the improved nonlinear estimate (4.34). This finishes the proof of Lemma 4.1.
Appendix C: Proof of Lemma 4.2
The Taylor expansions (B.6)–(B.8), (B.11) are still valid, which in turn imply that
In particular, \(\xi _n^{(m)}\) could be analyzed in a more precise way:
in which the regularity requirement (4.16) for the constructed profile \(\hat{{\textsf{N}}}\), as well as the \(\Vert \cdot \Vert _\infty \) a-priori error estimate for \({\tilde{n}}^m\), has been applied. A similar error bound could also be derived for \(\xi _n^{(m+1)}\), with the technical details skipped for the sake of brevity:
Then we arrive at the point-wise error bounds:
with \(M^{(0)} = 4 ( \epsilon _0^* )^{-2} \cdot (C^* +1 )\), and the phase separation estimates (4.30), (4.31), (4.51) have been applied. As a further consequence, the following inequality could be derived
On the other hand, the following estimates are straightforward:
with \(M^{(1)} = 4 C^* ( \epsilon _0^* )^{-2}\), in which the temporal regularity assumptions (4.17) for the constructed profile \(\hat{{\textsf{N}}}\), as well as the phase separation estimate (4.30), have been applied. Subsequently, the following inequalities are obvious:
Finally, a combination of (C.5) and (C.7) results in the refined estimate (4.53), by taking \({\tilde{C}}_6 = 2\,M^{(0)} + M^{(1)}\). Since both \(M^{(0)}\) and \(M^{(1)}\) only depend on \(C^*\), \(\epsilon _0^*\), the same dependence preserves for \({\tilde{C}}_6\).
Inequalities (4.54), (4.55), and (4.56) could be derived in a similar manner; the technical details are skipped for the sake of brevity.
Rights and permissions
Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.
About this article
Cite this article
Liu, C., Wang, C., Wise, S.M. et al. A Second Order Accurate, Positivity Preserving Numerical Method for the Poisson–Nernst–Planck System and Its Convergence Analysis. J Sci Comput 97, 23 (2023). https://doi.org/10.1007/s10915-023-02345-9
Received:
Revised:
Accepted:
Published:
DOI: https://doi.org/10.1007/s10915-023-02345-9
Keywords
- Poisson–Nernst–Planck (PNP) system
- Second order accuracy
- Positivity preserving
- Energy stability
- Optimal rate convergence analysis
- Rough error estimate and refined estimate