# Finite difference for nonlinear system of equation

$$\frac{\partial C_i}{\partial t} = D_i \nabla^2 C_i - \frac{I \cdot \nabla t_i}{z_i F} - \sum_{i'} \frac{z_{i'}}{z_i} D_{i'}\nabla \cdot (t_i\nabla C_{i'})$$

$$\nabla \eta = \frac{I}{\kappa} + \frac{F}{\kappa} \sum_i z_i D_i \nabla C_i$$ The dependent variables are $C_i$ and $\eta$ ($t_i$ is a function of $C_i$). In the 1D case.

I have some thought on difference scheme (TRBDF2) to try on the first equation. The second equation looks simple, but I'm afraid a simple centered difference will not work. I've try something like

$$\left. D_i\frac{\partial C_i}{\partial x}\right \vert_{k-1/2} = D_{i, k-1/2}\frac{C_{i,k}-C_{i,k-1}}{\Delta x} \approx \frac{D_{i, k}+D_{i,k-1}}{2}\frac{C_{i,k}-C_{i,k-1}}{\Delta x}$$

but with no luck when implementing in the nonlinear solver MINPACK. What will be a good scheme on differencing this two equations? Also, should I run analysis on stability issue every time I came up with a difference scheme before I compute them?

Another set of equation look like this, where the dependent variables are $C_i$, $\eta$, $i_2$, $\epsilon$ and $\epsilon_k$. I've try difference scheme and put on some efforts but in vain. Is there any general guideline as to the way of discretization? I'm using the method of line by the way.

$$\frac{\partial \epsilon C_i}{\partial t} = \nabla \cdot \biggl( \epsilon D_i \nabla C_i \biggr) t-\frac{i_2\cdot \nabla t_i}{z_i F} - \sum_j ai_j\biggl( \frac{t_i}{z_i F}+\frac{s_{i,j}}{n_j F} \biggr) - \sum_{i'} \frac{z_{i'}}{z_i}\epsilon D_{i'}\nabla \cdot (t_i\nabla C_{i'}) - R_i$$

$$\nabla \eta = -\biggl( \frac{I-i_2}{\sigma} \biggr) + \frac{i_2}{\epsilon \kappa} + \frac{F}{\kappa} \sum_i z_i D_i \nabla C_i$$

$$\nabla \cdot i_2 = a \sum_j i_j$$

$$\frac{\partial \epsilon}{\partial t} = -\sum_k \overset{\sim}{V_k} k_k\epsilon_k\biggl(\prod_i C_i^{\gamma_i,k}-K_{sp,k}\biggr)$$

$$\frac{\partial \epsilon_k}{\partial t} = \overset{\sim}{V_k} k_k\epsilon_k\biggl(\prod_i C_i^{\gamma_i,k}-K_{sp,k}\biggr)$$

EDIT To be more specific I'm now trying to solve the initial condition for $\eta$ and $i_2$ which is not given. Since there is no concentration gradient, in the separator (of the battery)

$$\nabla \eta = \frac{I}{\kappa}$$ $$i_2=I$$

In the porous cathode

$$\nabla \eta = -\biggl( \frac{I-i_2}{\sigma} \biggr) + \frac{i_2}{\epsilon \kappa}$$

$$\nabla \cdot i_2 = a \sum_j i_j$$

$i_j$ is Butler-Volmer equation, which in the Newton's method tends to explode when initial guess is bad.

$$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_j-U_j)\biggr) - \prod_i \biggl(\frac{C_i}{C_{i,o}}\biggr)^{q_{i,j}}\exp\biggl(-\frac{\alpha_{cj} F}{RT}(\eta_j-U_j)\biggr) \biggr]$$

I discretize the spatial derivative with second order difference, and the output shows a jump every other mesh point (I did not encounter this when I wrote my own code with Newton's method and analytic Jacobian, but I am afraid later on I will make another error when writing down analytic Jacobian since there are so many variables. So I try the MINPACK with numerical Jacobian but with no luck).

I felt like something is wrong in the boundary. Since $\eta$ is not continuous in the separator/porous cathode interface, so I imposed $\nabla \eta_s = \nabla \eta_c$. This is correct because $\eta=\phi_1-\phi_2$, and I assumed $\phi_2$ is continuous in the interface, $\phi_1$ does not exist in the separator but $\nabla \phi_1=0$ in the interface. Particularly in the interface mesh points I wrote

$$\frac{-3\eta_{s,n}+4\eta_{s,n-1}-\eta_{s,n-2}}{2\Delta x} = \frac{3\eta_{c,1}-4\eta_{c,2}+\eta_{c,3}}{2\Delta x}$$

Thanks!

• I'm still new to the forum, so please tell me if there is anything against the rule. – Henry Jul 22 '15 at 15:55

I'm guessing based on the notation that the application is electrolyte-related (batteries?), because the equations look like Nernst-Planck. Regarding the last system of equations: you're probably going to have to clarify your notation here. The variable $\eta$ looks like a potential, but I don't see how to close the system of equations you have. Usually, you're looking to solve for concentrations and electric potential.

I have some thought on difference scheme (TRBDF2) to try on the first equation...

TRBDF is usually used for time-stepping discretization. What you have depicted is a spatial discretization; don't use TRBDF for that.

Typically, these equations have only diffusive terms (that is, things that look like scaled second-order derivatives), in which case, any stencil that typically works for second-order derivatives is likely to work here. A collaborator of mine uses second-order centered differences for your spatial derivatives, which I would recommend to start unless you need high accuracy.

...but with no luck when implementing in the nonlinear solver MINPACK.

It's not clear to me what system of equations you are solving in MINPACK.

What I would do here is semi-discretize your PDE system by choosing a discretization in space. The result will be a set of ODEs, which can then be solved using standard software. Since your system of equations is diffusive in nature, choice of time-step and scheme is crucial for avoiding instabilities and achieving accurate results. You will want to choose an L-stable scheme to avoid ringing effects. (See Hairer and Wanner, Vol. 2 for details.) For example, TRBDF2 is an L-stable scheme.

In practice for your problem, an almost L-stable scheme might be good enough, and an A-stable scheme that lacks L-stability might also work if the spatial discretization is not too fine. (L-stability is really best here.)

In terms of software, I would use an ODE solver such as the schemes built into PETSc or Trilinos, or maybe use the implicit BDF methods in CVODE in SUNDIALS (or something else implicit, maybe an implicit Rosenbrock method like RODAS3). Using MINPACK will require you to implement your own time-stepping discretization and error control, which I think is usually a waste of effort, given the availability of good libraries that implement state-of-the-art versions of these features (and are well-tested!).

...is [sic] Butler-Volmer equation, which in the Newton's method tends to explode when initial guess is bad.

Yes, the Butler-Volmer equation has a phenomenally ill-conditioned Jacobian matrix.

I discretize the spatial derivative with second order difference, and the output shows a jump every other mesh point (I did not encounter this when I wrote my own code with Newton's method and analytic Jacobian, but I am afraid later on I will make another error when writing down analytic Jacobian since there are so many variables. So I try the MINPACK with numerical Jacobian but with no luck).

Don't use the numerical Jacobian. The numerical error in finite differencing is going to give you garbage. Use an analytical Jacobian. If you're afraid of making an error, derive it using a computer algebra system (e.g., Maple, Mathematica, SymPy, SageMath...).

Since $\eta$ is not continuous in the separator/porous cathode interface, so I imposed $\nabla\eta_{s} = \nabla_{c}$.

This statement is incorrect, because:

1. these quantities are both vectors (except in the 1-D case, in which case it's just misleading)
2. these should be current fluxes, and thus have the wrong units (so they should probably look something like $-\kappa_{s}\nabla\eta_{s} \cdot n$ instead of $\nabla\eta_{s}$, where $\kappa_{s}$ is the electrical conductivity of the separator, and $n$ is the interface outward normal)
3. you really need to explain your notation because most readers aren't going to be familiar with it, and even though I am, I'm largely just guessing based on experience -- people use different conventions all the time
4. assuming phase 1 is the solid phase and phase 2 is the liquid phase, and their electrical potentials are $\phi_{1}$ and $\phi_{2}$, then $\phi_{2}$ is going to be continuous across the separator-electrode interface (doesn't matter which one), but $\phi_{1}$ is not, and experiences a zero flux condition across the interface ($\nabla\phi_{1} \cdot n = 0$, where $n$ is the interface outward normal).

This combination of incorrect boundary conditions and general confusion probably contributes to your weird results.

I think PETSc is too much for me in the meantime. I am wondering which Fortran ode solver (subroutine) is the best alternative. Besides RODAS (Rosenbrock method), there are also RADAU5 (IRK) and the ones in netlib.org/ode (mostly BDF), and also DVODE which you don't recommend. Which is the most stable and efficient one?

I still recommend using PETSc. You are not going to know which method is most efficient without experimenting; PETSc makes it easy to experiment. Provided you use the library according to the documentation, changing numerical methods should be as easy as changing flags at the command line, plus PETSc includes features that make it easier to diagnose numerical problems (command-line switches that output residuals and performance information). Assuming you implement everything correctly, you should also be able to get decent parallel performance, since PETSc is portable and includes data structures you can use to write your program to take advantage of distributed-memory parallelism.

If you use any of those other libraries, most of them are written without any parallelism features at all (certainly RODAS and RADAU, probably also DVODE and LSODE; the only one I can think of that includes any parallel features would be FCVODE, which you are also avoiding). Thus, you will be limited in the size and dimension of problem you can solve, which will likely limit your science as well. Other developer efficiencies are also an issue. Both PETSc and Trilinos make it easy to experiment with different numerical methods through use of things like command-line switches. If, on the other hand, you want to switch time-steppers from RODAS to DVODE, you are likely to have to write whole new drivers, which takes time, and is more difficult to manage. Many of those old Fortran codes do not include any options for iterative linear solvers or preconditioners; again, lack of such features in these old codes will limit the size and dimension of problem you can solve, which again, limits your science.

There are problems for which RODAS, or RADAU, or other old Fortran codes are certainly appropriate (even if these libraries are dated); however, for this particular problem, I strongly recommend you use a framework like PETSc (or Trilinos), since it will not limit your work in the long run in a way that those other solvers will. In the case of PETSc, again: using Fortran (or Trilinos) is not an issue. C (or in the case of Trilinos, C++) is probably a better choice -- I would argue that knowing how to write C (respectively, C++) programs will make you more employable -- but you can use PETSc's (respectively, Trilinos') Fortran wrappers without knowing C (respectively, C++) , and PETSc (respectively, Trilinos) includes Fortran examples.

• Thanks for the reply! I have one semester course in C++ data structure, but I am more comfortable writing in Fortran for numerical simulation. Any good ODE solver you recommend in Fortran? I am a physics major in NTU by the way. – Henry Jul 23 '15 at 3:45
• RODAS3 is written in Fortran, and is on Ernst Hairer's web site. CVODE has a Fortran interface (called FCVODE); it's probably better to use FCVODE than use the predecessor to CVODE, a Fortran code called VODE (or DVODE, for the double-precision version). PETSc has Fortran wrappers for most functions, and if they don't, they usually write them upon request (because they have an automatic script that generates the wrappers). PETSc and CVODE are written in C; the most complicated thing about these libraries is their use of opaque data types (particularly opaque pointers). – Geoff Oxberry Jul 23 '15 at 4:32
• Actually, I want to give a try to C numerical coding with PETSc. It felt like after I fully understand the package, I will never write in Fortran. – Henry Jul 23 '15 at 5:47
• I think PETSc is too much for me in the meantime. I am wondering which Fortran ode solver (subroutine) is the best alternative. Besides RODAS (Rosenbrock method), there are also RADAU5 (IRK) and the ones in netlib.org/ode (mostly BDF), and also DVODE which you don't recommend. Which is the most stable and efficient one? Thanks! – Henry Jul 23 '15 at 11:46
• I am very interesting in DASKR since it includes a procedure for calculating consistent initial conditions. And also LSODE which the website lh3lh3.users.sourceforge.net/solveode.shtml#st-rk claims the best. – Henry Jul 23 '15 at 12:21