# How to implement Newton method in solving 1D PDE system? (ie. Poisson eq, continuity eq, drift-diffusion eq.)

I want to solve PDE system, which consists of Poisson equation, continuity equations for electron and hole with drift-diffusion equation numerically, by using method called Newton's method.

This method is explained in several literatures, including

(But I'm not sure I understand this method at all)

$$\begin{gathered} \frac{\partial }{{\partial x}} \left( {\varepsilon \frac{\partial }{{\partial x}}V} \right) = - q\left( {p - n + N_D^ + - N_A^ - } \right) \hfill \\ \frac{\partial }{{\partial x}}{J_n} = q\left( {\frac{{\partial n}}{{\partial t}} + \frac{{n - {n_0}}}{{{\tau _n}}} - {G_{op}}} \right) \hfill \\ \frac{\partial }{{\partial x}}{J_p} = - q\left( {\frac{{\partial p}}{{\partial t}} + \frac{{p - {p_0}}}{{{\tau _p}}} - {G_{op}}} \right) \hfill \\ {J_n} = - qn{\mu _n}\frac{\partial }{{\partial x}}V + q{D_n}\frac{\partial }{{\partial x}}n \hfill \\ {J_p} = - qp{\mu _p}\frac{\partial }{{\partial x}}V - q{D_p}\frac{\partial }{{\partial x}}p \hfill \\ \end{gathered}$$

Arranging these, I got $$\begin{gathered} 0 = \frac{\partial }{{\partial x}} \cdot \left( {\varepsilon \frac{\partial }{{\partial x}}V} \right) + q\left( {p - n + N_D^ + - N_A^ - } \right) \hfill \\ 0 = \frac{{\partial n}}{{\partial t}} + {\mu _n}\frac{{\partial n}}{{\partial x}}\frac{{\partial V}}{{\partial x}} + {\mu _n}n\frac{{{\partial ^2}V}}{{\partial {x^2}}} - {D_n}\frac{{{\partial ^2}n}}{{\partial {x^2}}} + \frac{{n - {n_0}}}{{{\tau _n}}} - {G_{op}} \hfill \\ 0 = \frac{{\partial p}}{{\partial t}} - {\mu _p}\frac{{\partial p}}{{\partial x}}\frac{{\partial V}}{{\partial x}} - {\mu _p}p\frac{{{\partial ^2}V}}{{\partial {x^2}}} + {D_p}\frac{{{\partial ^2}p}}{{\partial {x^2}}} + \frac{{p - {p_0}}}{{{\tau _p}}} - {G_{op}} \hfill \\ \end{gathered}$$

Now, what I understand from general Newton's method for solving equation system is following; I have system of equation \begin{aligned} f_1(x,y,z)=0 \\ f_2(x,y,z)=0 \\ f_3(x,y,z)=0 \end{aligned} \qquad\rightarrow\qquad \mathbf{F}(x,y,z)=0 and initial value (starting value) $$\mathbf{X}_0=(x_0,y_0,z_0)$$ Then I create Jacobian Matrix $$\mathbf{J}(x,y,z)= \left[ \begin{array}{ccc} \partial_xf_1 & \partial_yf_1 & \partial_zf_1\\ \partial_xf_2 & \partial_yf_2 & \partial_zf_2\\ \partial_xf_3 & \partial_yf_3 & \partial_zf_3 \end{array} \right]$$ The $\mathbf{X}$ that give $\mathbf{F}=0$ is $$\mathbf{X}_{n+1}=\mathbf{X}_n-\mathbf{J}^{-1}(\mathbf{X}_n)\mathbf{F}(\mathbf{X}_n) \qquad (n=0,1,2...)$$

I got this idea, and Now I want to apply this into the poisson, continuity, drift-diffusion eq. system. First, I have my system of equation, with independent variables $V, n, p$.

\begin{aligned} 0 =& f_1(V,n,p)= \frac{\partial }{{\partial x}} \left( {\varepsilon \frac{\partial }{{\partial x}}V} \right) + q\left( {p - n + N_D^ + - N_A^ - } \right) \\ 0 =& f_2(V,n,p)= \frac{{\partial n}}{{\partial t}} + {\mu _n}\frac{{\partial n}}{{\partial x}}\frac{{\partial V}}{{\partial x}} + {\mu _n}n\frac{{{\partial ^2}V}}{{\partial {x^2}}} - {D_n}\frac{{{\partial ^2}n}}{{\partial {x^2}}} + \frac{{n - {n_0}}}{{{\tau _n}}} - {G_{op}} \\ 0 =& f_3(V,n,p)= \frac{{\partial p}}{{\partial t}} - {\mu _p}\frac{{\partial p}}{{\partial x}}\frac{{\partial V}}{{\partial x}} - {\mu _p}p\frac{{{\partial ^2}V}}{{\partial {x^2}}} + {D_p}\frac{{{\partial ^2}p}}{{\partial {x^2}}} + \frac{{p - {p_0}}}{{{\tau _p}}} - {G_{op}} \end{aligned} I assumed that the initial condition might be V(x)=\left\{ \begin{aligned} 0 \qquad & (x_1<x<0)\\ k_BT\ln\left(\frac{N_AN_D}{n_i^2}\right)\qquad &(0<x<x_2)\\ \end{aligned}\right. \qquad n(x)=\left\{ \begin{aligned} N_D \qquad & (x_1<x<0)\\ \frac{n_i^2}{N_A} \qquad &(0<x<x_2)\\ \end{aligned}\right. \qquad p(x)=\left\{ \begin{aligned} \frac{n_i^2}{N_D} \qquad & (x_1<x<0)\\ N_A \qquad &(0<x<x_2)\\ \end{aligned}\right.

I got stuck when I was creating Jacobian $\mathbf{J}_{11}$ $$\partial_Vf_1=\frac{\partial }{\partial V}\left[ \frac{\partial }{{\partial x}} \left( {\varepsilon \frac{\partial }{{\partial x}}V} \right) + q\left( {p - n + N_D^ + - N_A^ - } \right) \right] \stackrel{?}{=} \frac{\partial }{{\partial x}} \left( {\varepsilon \frac{\partial }{{\partial x}}} \right)$$ Apparently, I got an Jacobian element which is just operator! Am I proceeding correctly? If I may go further, I have my Jacobian that looks like $$\mathbf{J}(V,n,p)= \left[ \begin{array}{ccc} \frac{\partial }{{\partial x}} \left( {\varepsilon \frac{\partial }{{\partial x}}} \right) & -q & q \\ \mu_n\frac{\partial^2 n}{\partial x^2}+mu_n n\frac{\partial^2 }{\partial x^2} & \frac{{\partial }}{{\partial t}} + {\mu _n}\frac{{\partial^2 V}}{{\partial x^2}} + {\mu _n}\frac{{{\partial ^2}V}}{{\partial {x^2}}} - {D_n}\frac{{{\partial ^2}}}{{\partial {x^2}}} + \frac{1}{{{\tau _n}}} & 0 \\ -\mu_p\frac{\partial^2 p}{\partial x^2}-mu_p p\frac{\partial^2 }{\partial x^2} & 0 & \frac{{\partial }}{{\partial t}} - {\mu _p}\frac{{\partial^2 V}}{{\partial x^2}} - {\mu _p}\frac{{{\partial ^2}V}}{{\partial {x^2}}} + {D_p}\frac{{{\partial ^2}}}{{\partial {x^2}}} + \frac{1}{{{\tau _p}}} \end{array} \right]$$

$$= \left[ \begin{array}{ccc} \frac{\partial }{{\partial x}} \left( {\varepsilon \frac{\partial }{{\partial x}}} \right) & -q & q \\ \mu_n\frac{\partial^2 n}{\partial x^2}+\mu_n n\frac{\partial^2 }{\partial x^2} & \frac{{\partial }}{{\partial t}} + 2{\mu _n}\frac{{\partial^2 V}}{{\partial x^2}} - {D_n}\frac{{{\partial ^2}}}{{\partial {x^2}}} + \frac{1}{{{\tau _n}}} & 0 \\ -\mu_p\frac{\partial^2 p}{\partial x^2}-\mu_p p\frac{\partial^2 }{\partial x^2} & 0 & \frac{{\partial }}{{\partial t}} - 2{\mu _p}\frac{{{\partial ^2}V}}{{\partial {x^2}}} + {D_p}\frac{{{\partial ^2}}}{{\partial {x^2}}} + \frac{1}{{{\tau _p}}} \end{array} \right]$$ I don't have any idea how to apply inverse matrix of this "operator" matrix

Is this correct result?

• You're going to need to do more than linearize system of PDEs, you're also going to need to discretize the derivatives in space and time in order to find a computational solution. There are lots of different ways to do that. – Bill Barth Mar 13 '16 at 22:52
• Can you show me how to do it actually? – user65452 Mar 13 '16 at 23:06
• Nope. There are whole books on this. – Bill Barth Mar 14 '16 at 2:37
• Then can you recommend me some books on this?, I mean, rather in details. – user65452 Mar 14 '16 at 4:04
• Actually, you need much more than only Newton to solve this. You need methods for discretization in space (e.g. FEM, FDM, FVM, spectral, etc.) and time (e.g. Euler, Crank-Nicholson, etc.). – Bill Barth Mar 15 '16 at 15:49