# Numerical methods for coupled stiff PDEs

I'm dealing with a set of nonlinear coupled PDEs that have the form:

\begin{align} \frac {\partial y_1}{\partial t} &= y_2y_3 - y_1 \tag{1}\\ \frac {\partial y_2}{\partial t} &= y_1y_3 - y_2 \tag{2}\\ \frac {\partial^2 y_3}{\partial t^2} - \frac {\partial y_3}{\partial t} &= \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} - \frac {\partial^2 y_2}{\partial t^2} + \frac {\partial y_2}{\partial t} \tag{3} \end{align}

with initial/boundary conditions:

\begin{align} y_1(z,0) &= 10^{-5}\\ y_2(z,0) &= 10^4\\ y_3(z,0) &= 0\\ \frac {\partial y_3}{\partial t} (z,0) &= 0\\ 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}

An initial approach has been to use a finite difference scheme. One being where I simply use a method of lines and treat it as a set of ODEs. But with this setup, I am having trouble formulating the above equations in such a way to actually use a stiff solver such as Gear's method since I lack an explicit expression for $$\frac {\partial y_3}{\partial t}$$

By taking the time derivative of $(2)$, I can "simplify" $(3)$ to:

$$\frac {\partial^2 y_3}{\partial t^2} + (y_1-1) \frac {\partial y_3}{\partial t} = \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} + \frac {\partial y_2}{\partial t} - y_3 \frac {\partial y_1}{\partial t}$$

but that still doesn't help too much.

Any suggestions on methods/approaches to solve the above equations?

• Don't you also need an initial condition for $\frac{\partial y_2}{\partial t}$ since there is a second derivative of $y_2$ in equation 3? – Bill Greene Aug 17 '16 at 19:59
• @BillGreene $\frac{\partial y_2}{\partial t}$ is set by equation [2] to $y_1y_3-y_2$, which is the IC. – Kirill Aug 17 '16 at 20:06
• What is the spatial domain? What are your boundary conditions? If they are simple enough, there may be better options. – David Ketcheson Aug 18 '16 at 11:46
• @DavidKetcheson At this moment, I am planning to test multiple scales. Therefore $z$ may range from $10^3 - 10^{15}$, and $t$ from $10^{-3} - 10^3$ – Mathews24 Aug 18 '16 at 15:14
• @DavidKetcheson The boundary conditions are included in the main post now. – Mathews24 Aug 18 '16 at 17:16

The first step is to add one more equation defined as:

$$\frac{\partial y_3}{\partial t} = v$$

Then if you substitute equations 1 and 2 and the new equation into your modified equation 3 you get:

$$\frac{\partial v}{\partial t}=(1-y_1)v + \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_1y_3-y_2 -y_3(y_2y_3-y_1)$$

• Earlier, I did this with RK4 but it did not work. I then switched to Gear's method but was having trouble actually finding all the necessary terms since to find the updated $\frac {\partial v}{\partial t}$ requires me finding the updated $v$ —but this term is implicitly defined by the earlier term. Any recommendations? (Perhaps, I am missing something obvious; an example of finding the necessary terms would be helpful.) – Mathews24 Aug 17 '16 at 21:07
• Also, besides a central difference scheme, is there a better way to handle the spatial derivatives? – Mathews24 Aug 17 '16 at 21:07
• Having the v term on the RHS is no different from from having the $y_i$ terms on the RHS. Obviously, equation [4] and the modified equation [3] are coupled in v but this should present no problems. As part of debugging this, I would simply try backward Euler as an integrator; it should work perfectly fine. Similarly, central difference for the spatial derivatives should also work perfectly fine. – Bill Greene Aug 17 '16 at 21:32
• Thank you for the response. I'll try it out and get back with the results. – Mathews24 Aug 17 '16 at 21:47
• I suggested trying backward Euler on all four equations because it is very simple and therefore easier to debug. When you get that working, you might want to change to a more sophisticated integrator. – Bill Greene Aug 18 '16 at 11:29