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

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.

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.

• Does this help, scicomp.stackexchange.com/q/20710/3691 set $\theta=1/2$ for Crank-Nicolson. Commented Sep 13, 2015 at 13:54
• What coupled PDE? Commented Sep 13, 2015 at 14:10
• Welcome to SciComp.SE! Yes, please do add more details about the equation and about what precisely you have difficulties with -- is it the finite difference discretization, the Crank-Nicolson time-stepping, or the Newton iteration? Commented Sep 13, 2015 at 19:03
• Thanks all. I’m going to solve Poisson-Nernst-Planck system of equations. The Nernst-Planck equation is same as below: ∂ni/∂t=Ai*((∂^(2 ) ni)/(∂x^2 ))+Bi*(∂ni/∂x*(∂∅)/∂x)+Bi*(ni*(∂^2∅)/(∂x^2 )) .So the Poisson equation that couples by the previous equation is: ε (∂^2∅)/(∂x^2 )=F∑〖ni*zi〗 So how could I implement Poisson equation into Nernst-Planck and after that implement Newton-Raphson method for linearize it. Commented Sep 14, 2015 at 7:33
• This is a Differential-Algebraic Equation (DAE) since the Poisson equation has no time term. The simplest thing to do is to solve the constraints (Poisson) and then step the equations (NP). If you want something higher order, look at the book by Ascher and Petzold. Commented Sep 16, 2015 at 11:18

$$\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.