# 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

https://nanohub.org/resources/1565/download/ddmodel_solution_details_word.pdf

(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

## 2 Answers

Here is a really good book on implementing these equations:

Analysis and Simulation of Semiconductor Devices by S. Selberherr

Before solving the full set of drift diffusion equations, an initial guess would be solving only the Poisson equation at equilibrium. This is also discussed in the book.

Here's an extended discussion about how to do it for a single, nonlinear PDE: https://dealii.org/developer/doxygen/deal.II/step_15.html . You may also want to watch lectures 31.5 and following on nonlinear PDEs here: http://www.math.tamu.edu/~bangerth/videos.html .

(Disclaimer: I co-wrote that program and these are my lectures.)