# Numerically solving a system of stiff nonlinear PDEs

I am attempting to numerically solve the following:

\begin{align} \frac {\partial y_1}{\partial t} &= i(y_2y_3 - y_2^*y_3^*) - y_1 \tag{1}\\ \frac {\partial y_2}{\partial t} &= y_1^*y_3 - y_2 \tag{2}\\ \frac {\partial y_3}{\partial t} &= y_4 \tag{3} \\ \frac {\partial y_4}{\partial t} &= y_4+ \frac {\partial^2 y_3}{\partial z^2} + \frac {\partial y_3}{\partial z} + \frac {\partial^2 y_2}{\partial z^2} + \frac {\partial y^*_2}{\partial z} + y_1^*y_4 + y_3[i(y_2y_3 - y_2^*y_3^*) - y_1] - iy_1y_3 -y_2 \tag{4} \end{align}

The initial/boundary conditions are:

\begin{align} y_1(z,0) = 10^{-5}; y_2(z,0) = 10^4; y_3(z,0) = 0; y_4 (z,0) = 0; \end{align} \begin{align} y_1(0,t) = 10^{-5}e^{-t}; y_2(0,t) = 10^{4}e^{-t}; y_3(0,t) = 0; \frac {\partial y_3}{\partial z} (0,t) &= 0 \end{align}

With a current approach using method of lines, I implemented a central difference scheme for the spatial derivatives and Runge-Kutta (4th order) for the time derivatives. This worked but only for particular conditions within a narrow time range (after which the output blew up).

I have attempted using backward Euler instead of RK4 for each of the equations, although since it involves complex conjugates, I end up needing to split the real and imaginary parts which then yields a system of 8 equations. Splitting the real and imaginary parts is no issue itself, however, since these equations are nonlinear, I cannot isolate all the desired variables and insert them sequentially into the last equations.

Do you have any suggestions on how to tackle this problem? Any better-suited explicit methods I should try for this problem? Advice on formulating/implementing an implicit method to solve the above equations? I am using Python and have looked into existing packages, although have had no luck with a couple (e.g. FiPy).

• Have you considered trying to use an explicit adaptive time stepping scheme? Something like Runge-Kutta-Cash-Karp? – spektr Aug 29 '16 at 17:15
• Have you looked at LSODE in ODEPACK? That's what it's for, and you can call it from just about any language. – Mike Dunlavey Aug 29 '16 at 17:25
• @choward I am in the process of implementing that right now. Still working on it and was just curios if people had other suggestions. – Mathews24 Aug 29 '16 at 17:59
• @MikeDunlavey I have heard about LSODE and went through its documentation, and was just having trouble trying to implement it with my set of DEs with terms such as spatial derivatives. By making these real-valued first-order derivatives, I would have 10 nonlinear DEs in this system and it just didn't seems like a good approach for large time- and spatial-domains. I'm still in the process of trying new things, but do you think this is a good solution? – Mathews24 Aug 29 '16 at 19:12
• Consider to approximate the first derivatives in the 4th equation not by central difference, but one-sided one, I would guess the forward difference might help, if your blow up means some oscillations grow up. You decrease the accuracy, but you may increase stability. I do not have experiences with your type of PDEs, it is not usual to have "convection-diffusion" term of $y_3$ (and $y_2^*$ and $y_2$?!) in the equation for $y_4$. If e.g. so called Peclet number is important also in your case, your space discretization step, say $\Delta z$, shall be less than 1/2 when using central difference. – Peter Frolkovič Aug 30 '16 at 8:24