4
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Does this help, scicomp.stackexchange.com/q/20710/3691 set $\theta=1/2$ for Crank-Nicolson. $\endgroup$ – boyfarrell Sep 13 '15 at 13:54
  • 2
    $\begingroup$ What coupled PDE? $\endgroup$ – Bill Barth Sep 13 '15 at 14:10
  • 3
    $\begingroup$ 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? $\endgroup$ – Christian Clason Sep 13 '15 at 19:03
  • $\begingroup$ 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. $\endgroup$ – Milad Rezaee Sep 14 '15 at 7:33
  • $\begingroup$ 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. $\endgroup$ – Matt Knepley Sep 16 '15 at 11:18
2
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.