# Ill-conditioned Jacobian matrix from Nernst-Planck equation with Butler-Volmer reactions

The governing equations are listed here of my notes on page 4. It's a reproduction of other's paper which solves the equations with COMSOL. The problems arise when I want to solve for the consistent initial condition of the algebraic components, $\phi_1$ and $\phi_2$. They consist of a kinetic reaction term for five reactants

\begin{equation} a \sum_j i_j \end{equation} \begin{equation} i_j = i_{o,j}\biggl[ \prod_i \biggl(\frac{C_i}{C_{i,o}}\biggr)^{p_{i,j}}\exp\biggl(\frac{\alpha_{aj} F}{RT}(\eta-U_j)\biggr) - \prod_i \biggl(\frac{C_i}{C_{i,o}}\biggr)^{q_{i,j}}\exp\biggl(-\frac{\alpha_{cj} F}{RT}(\eta-U_j)\biggr) \biggr] \end{equation} where $\eta=\phi_1-\phi_2$.

First try for these two algebraic equations on MINPACK gives $l^2$ norm of the residual around $1.0E-2$ for one dimension case $=500$. Then I twisted the equations so they become one equation of variable $\eta$ (not two coupled equations), the equations listed here on page 4. Remarkably the $l^2$ norm drops to $1.0E-5$. Nevertheless, when I integrate over time they failed on the very first step. It seems that when the equations are all coupled together with this kinetic term, the convergence of the Newton's method failed (I've only tried on the solver DASKR (DASSL's variant) and RADAU5). I'm using method of line and the space discretization is second order. And the Neumann boundary conditions are all coupled back to the governing equations with the ghost points. The Jacobian matrices are generated by Tapenade. The problem is very close to this one.

Furthermore I think (with my little numerical knowledge) the equation of this kind \begin{equation} \sigma \nabla^2 \phi_1 = a\sum_j i_j \end{equation} resembles a singular perturbed equation for the coefficient $\sigma/a$ can be of order $1.0E-10$ on the sulfur cathode. $a$ is the surface area of the pore walls per unit volume of the total electrode. From the output data at the boundary there is a one point drop from a relatively smooth inner domain. I'm thinking pseudo-spectral method might solve the problem?

Another question is that is it perfectly fine to impose boundary conditions (specifically Neumann type) on PDE as a discretized algebraic equations without putting back into the governing equation with ghost point? I feel like somehow the boundary point won't satisfy the governing equation this way. And when implementing in the IDE or DAE solver I need to explicitly specify they are algebraic term though.

One last question: do the initial conditions need to satisfy the boundary conditions? I think most system don't. Or else I have to let the applied current to be zero first and then at $t>0$ crank up to some value, making a discontinuity at the first step. I've found one paper describe this discontinuous initial condition, DAEs not ODEs page 371, but I still don't know the technique in solving them.

I've been working this whole summer vacation from scratch and this is my very first time solving a real (non-toy) numerical problems. Any suggestion would be greatly appreciated.

EDIT

Quad precision without precondition just solves my first problem regarding ill-conditioned matrix. The third problem (inconsistent initial and boundary condition) is also solved by making initial guessing on the algebraic variables sufficiently close to the values when the source term jumps. But the problem is quad precision implementation is more than 10x slower than 8 bytes double precision. Will it help if I reorder the Jacobian matrix so that it is banded (currently the linear solver is LAPACK with dense matrix)? Or I just need to learn PETSc and make use of the parallel feature and sparse solver? Also, how fined should it be so that the solution can considered fine-grid solution (and hence make error estimation from that)? The second problem above in the main context still puzzles me though.

EDIT.II

Reaction current $\nabla \cdot i_2$ figure from a review paper on porous electrode. The author define a dimensionless current density $\delta = \frac{\alpha_a FIL}{RT}(\frac{1}{\kappa}+\frac{1}{\sigma})$ which monitor how non-uniform the current distribution can be. When $\delta=100$, the current at the boundary is very steep already. The system I studied has even higher of this parameter though.

EDIT.III

I solve the discontinuity problem regarding the Poisson equation. From the dimensionless parameter above, there is a term $L$. I put $1$ in my code which is non-physical for it means 1m electrode length for consistent unit. After lowering it down to around $1E-5$, the source term becomes smooth. Thanks for Geoff's a priori thought that made me look closer to the physical quantity.

I'm thinking pseudo-spectral method might solve the problem?

A priori, I don't see why it would. As a continuous expression, the Butler-Volmer equation varies tremendously with small changes in potential. Any discretization inherits this behavior, so I think you're stuck with ill-conditioning no matter what. Maybe you could nondimensionalize the expressions, which would help somewhat. Nondimensionalization doesn't seem to be common practice in the battery simulation literature.

Another question is that is it perfectly fine to impose boundary conditions (specifically Neumann type) on PDE as a discretized algebraic equations without putting back into the governing equation with ghost point?

I don't see why you couldn't express boundary conditions as explicit algebraic equations, but it would behoove you to eliminate any algebraic equations if possible, since ODE solvers are generally faster than DAE solvers unless the DAE system has significantly sparser Jacobian (or, more generally, the linear systems being solved during the DAE solution algorithm are solved faster than the corresponding linear systems being solved during the ODE solution algorithm).

One last question: do the initial conditions need to satisfy the boundary conditions?

No, they don't, although when this happens, there is normally a rapid transient associated with the solution satisfying the boundary conditions immediately after time $t = 0$.

But the problem is quad precision implementation is more than 10x slower than 8 bytes double precision.

Yes, this issue is the drawback in using quad precision.

Will it help if I reorder the Jacobian matrix so that it is banded (currently the linear solver is LAPACK with dense matrix)?

Well, sure, it would be faster, but you're essentially solving a Poisson equation with a source term, so a preconditioned iterative solver would be faster.

Or I just need to learn PETSc and make use of the parallel feature and sparse solver?

Forget the parallel feature for a second. Mostly, you need to learn how to use PETSc for the sparse solvers. An algebraic multigrid preconditioned GMRES solver will probably work for the linear systems; a similar solver combination worked well for me when solving a Poisson equation with source terms from Butler-Volmer.

Also, how fined should it be so that the solution can considered fine-grid solution (and hence make error estimation from that)?

Run a convergence study and refine your mesh repeatedly until the difference between solutions on successively refined meshes is small enough for you.

• The solution ($\phi_1$ and $\phi_2$) of the Poisson equations are smooth, but the problem is that when I plug them back into the source there is a discontinuous boundary point. It happens when the applied current is large and conductivity is relatively small (on sulfur system). It becomes smooth in the opposite setting though. – Witty Viper of Hidden Glen Aug 24 '15 at 9:42
• I think the case is close right? I have to use PETSc after all, or my science will be limited. What if preconditioned iteration failed? Then quad precision will seemingly be the the only solution, but it is slow on software... Nevertheless sparse solver on direct method is sufficient for 1D problems I think. For this kind of system, are there needs to implement in 2D or even 3D? If needed, why? Equation $\nabla \cdot i_1 + \nabla \cdot i_2 = 0$ becomes $i_1 + i_2 = I$ in 1D. What's the impact? If the question is nontrivial I will post a question on this. Thanks a lot! – Witty Viper of Hidden Glen Aug 24 '15 at 17:29

This is a part of a much bigger problem.

Look at

http://arxiv.org/abs/1502.07251

and

http://arxiv.org/ftp/arxiv/papers/1502/1502.07251.pdf

Ever yours Bob Eisenberg Chairman emeritus Molecular Biophysics Rush University Medical Center Chicago

• Welcome to SciComp! Stack Exchange tends to discourage link-only answers in favor of answers that use links as supporting materials. For an answer like yours, I'd advise you to summarize the main results of your paper, and explain why, for instance, the inconsistency between current conservation and mass action kinetics results in ill-conditioned systems. With a compelling summary acting as an abstract, I think you'll catch people's interest, and they will then click on the link you've provided to read your paper. – Geoff Oxberry Aug 23 '15 at 21:09
• Having skimmed this paper, and having read the first 8-10 pages in detail, I don't think it addresses the question. This paper does not discuss ill-conditioning, nor does it use the phrase "condition number" anywhere. – Geoff Oxberry Aug 24 '15 at 8:20