How could we solve coupled PDE with finite difference method and Newton-Raphson method?

Question

I'm trying to solve coupled PDE by Crank-Nicolson (CN) and Newton-Raphson method with MATLAB. I have used CN method but not for coupled problem. Please if someone could help let me know to add more details about the equation.

Thanks in advance

Detail of coupled equation:

I’m going to solve Poisson-Nernst-Planck system of equations. The Nernst-Planck equation is $$\frac{\partial n_i}{\partial t}=A_i\frac{\partial^2 n_i}{\partial x^2}+B_i\frac{\partial n_i}{\partial x}\frac{\partial \phi}{\partial x}+B_i(n_i\frac{\partial^2 \phi}{\partial x^2}),$$ where $n_i$ is the concentration of specific ion, $A_i$ and $B_i$ is the coefficient of equation for specific ion, $\phi$ is the electrical potential. The Poisson equation that couples by the previous equation is $$\epsilon \frac{\partial^2 \phi}{\partial x^2}=F\sum(n_i*z_i),$$ where $\epsilon$ is the permittivity of medium, $F$ is the Faraday constant, $n_i$ is the concentration of specific ion and $z_i$ is the ionic charge. So how could I implement Poisson equation into Nernst-Planck and after that implement Newton-Raphson method for linearize it. Or do we need Newton-Raphson method to solve this equation? or we could solve it just by CN method. Please let me know if anybody needs more details or if something is unclear.

Answer 1

As mentioned by Matt Knepley, this is naturally formulated as a system of partial differential algebraic equations. Because you're in Matlab, you could consider doing the spatial discretization yourself (e.g. finite difference, finite volume, finite element) to obtain a system of DAE's, then use the method of lines to step forward in time. You can discretize both the Nernst-Planck and Poisson equations the same way in space. This should give a list of equations like

$$\frac{\partial n_{i,x_j}}{\partial t} = f(\{n_{i,x_j}\}, \{\phi_{x_j}\})$$ $$0 = g(\{n_{i,x_j}\}, \{\phi_{x_j}\})$$

where the $x_j$ refers to the value of the variable at a particular discretized $x$ coordinate, and the functions $f$ and $g$ depend on your spatial discretization and will be different for each spatial location. If you have $N$ species, you should have $N$ of the first equation at each grid point and one of the second, leading to the same number of unknowns and equations (I don't know about your boundary conditions, so you can treat them as extra equations, substitute them in algebraically to the discretized equations at the edges, etc.).

If you're going to use CN time stepping, this will naturally lead to a non-linear system of equations you have to solve at each time step, which you could do using the Newton-Raphson method. However, unless you are particularly set on CN time stepping, I would suggest considering Matlab's ode15s (see documentation), which is reasonably good at time stepping through coupled, nonlinear (index-1) DAE's like this one. It uses variable order, variable size stepping, and I've found it to work well for coupled Nernst-Planck style equations. To couple the differential and algebraic equations in ode15s, you will need to pass a (singular) mass matrix to ode15s. It will take care of doing the Newton-Raphson iterations for you, and you can either let it estimate the Jacobian numerically or provide a function which provides the analytical Jacobian.