# Stabilization of solution to one-dimensional system of PDE

I am trying to solve numerically next PDE system: $$\frac{\partial c}{\partial t}=\epsilon\frac{\partial}{\partial x}(\frac{\partial c}{\partial x}+\rho\frac{\partial \varphi}{\partial x}+\frac{vc}{1-vc}\frac{\partial c}{\partial x})$$ $$\frac{\partial \rho}{\partial t}=\epsilon\frac{\partial}{\partial x}(\frac{\partial \rho}{\partial x}+c\frac{\partial \varphi}{\partial x}+\frac{v\rho}{1-vc}\frac{\partial c}{\partial x})$$ $$-\epsilon^2\frac{\partial \varphi^2}{\partial x^2}=\rho$$ with completely blocking boundary conditions at x= ±1; $v$ and $\epsilon$ are constants. However, at some point of calculation the oscillations near the boundaries occur due to the too big gradients near the boundaries. I have found some information that there are several techniques for handling numerical instabilities. (Petrov-Galerkin (SUPG) or Galerkin least-squares (GLS) method). Unfortunately all of them usually described for the case of simple Convection-Diusion equation. Can anybody help me how to adjust some of those methods (or maybe another method) for my case? Also I need to add, that problem is not in the fact that $(1-vc)$ may tended to zero in some points.

• How do you solve the non-linearity? – Dr_Sam Jun 24 '13 at 7:31
• I don't quite understand your question. I can say only, that I'm using special software for solving PDE, not my own C++ code. – oleksandrv Jun 24 '13 at 11:57
• You have a non-linear system of PDEs, so I wanted to know how you linearize it to be able to solve using finite elements. – Dr_Sam Jun 24 '13 at 12:19
• Newton's method. – oleksandrv Jun 24 '13 at 13:44

• Maybe if you plotted the troublesome solution for us, gave the values of $\nu$ and $\epsilon$, and gave the mathematical form of the boundary conditions, we could help out a little more. – Bill Barth Jun 23 '13 at 23:19
• Boundary conditions can be written in the next form: $$\frac{\partial c}{\partial x}+\rho\frac{\partial \varphi}{\partial x}+\frac{vc}{1-vc}\frac{\partial c}{\partial x}=0$$ $$\frac{\partial \rho}{\partial x}+c\frac{\partial \varphi}{\partial x}+\frac{v\rho}{1-vc}\frac{\partial c}{\partial x}=0$$ $$\varphi\pm\delta\epsilon\frac{\partial \varphi}{\partial x}=\pm v$$ where $\delta$ is a constant which measures the effective thickness of the surface insulating layer (according to the physics model). – oleksandrv Jun 24 '13 at 9:27
• To be honest I'm using these equations in a little bit another form (which is a little harder for analysis, for my opinion) $$\frac{\partial cn}{\partial t}=\frac{\partial }{\partial x}(A*cn \frac{\partial }{\partial x}(-B*\phi+C*(\log{cn}-\log{(1-D*(cp+cn))}))$$ $$\frac{\partial cp}{\partial t}=\frac{\partial }{\partial x}(A*cp \frac{\partial }{\partial x}(B*\phi+C*(\log{cp}-\log{(1-D*(cp+cn))}))$$ $$\frac{\partial \phi^2}{\partial x^2}=D(cp-cn)$$ However, this system can be easely converted to previous form with simple change of variables. – oleksandrv Jun 24 '13 at 9:27
• Boundary conditions: $$\mathbf{n}\cdot((A*cn \frac{\partial }{\partial x}(-B*\phi+C*(\log{cn}-\log{(1-D*(cp+cn))})))=0$$ $$\mathbf{n}\cdot((A*cp \frac{\partial }{\partial x}(B*\phi+C*(\log{cp}-\log{(1-D*(cp+cn))})))=0$$ $$\mathbf{n}\cdot(80*\frac{\partial \phi}{\partial x})=0$$ Values of parameters: A=7.26E-13 B=96485 C=2436 D=2.21E-4 Initial conditions: $cn$=$cp$=250; $\phi$=3 (smoothly increasing from 0 to 3) Oscillations: imgur.com/OY85XOE – oleksandrv Jun 24 '13 at 9:32