# Numerical scheme for coupled PDEs

I am trying to solve the three coupled PDEs;

$$\frac{\partial{Q}}{\partial{t}} = -RaPra^2\theta - Pra^2Q + Pr\frac{\partial^2{Q}}{\partial{z}^2}, \ \ \ \ \ \ \ \ \ (1)$$

$$\frac{\partial{\theta}}{\partial{t}} = w - a^2\theta + \frac{\partial^2{\theta}}{\partial{z}^2}, \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (2)$$

$$Q = -a^2w + \frac{\partial^2{w}}{\partial{z}^2}, \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (3)$$

where $$Ra, Pr, a$$ are constants. I want a method which is second order accurate in both time and space. I tried using the Crank-Nicolson scheme for the first two equations;

$$\frac{Q^{n+1} - Q^n}{\Delta{t}} = \frac{1}{2}(F^{n+1} + F^n)$$

$$\frac{\theta^{n+1} - \theta^n}{\Delta{t}} = \frac{1}{2}(G^{n+1} + G^n)$$,

and a centered space finite difference scheme for the third equaiton - $$F$$ and $$G$$ are the right hand sides of Eq.1 and Eq.2 respectively. My problem is that when using the Crank-Nicolson scheme I do not know $$\theta^{n+1}$$ or $$Q^{n+1}$$ and therefore $$w^{n+1}$$. So far, I have just used;

$$\frac{Q^{n+1} - Q^n}{\Delta{t}}= -RaPra^2\theta^n + \frac{1}{2}(H^{n+1} + H^n),$$

$$\frac{\theta^{n+1} - \theta^n}{\Delta{t}} = w^n + \frac{1}{2}(K^{n+1} + K^n)$$,

which isn't fully second order in time. I had an idea that was to solve these equations using the scheme above. Then use RK2 or a similar predictor-corrector method where the above scheme is my predictor. Does this make sense? What scheme can I use to solve these equations which will be accurate to second order in space and time?

I am similar with finite difference methods and Runge-Kutta methods so anything involving these would be the best.

In a standard method of lines approach we discretize in space first (using centered finite differences as you mentioned). We define new vector-valued time-dependent functions $$\boldsymbol{Q}, \boldsymbol{\theta}, \boldsymbol{w}$$ with $$\boldsymbol{Q}(t) = (Q_1(t),Q_2(t),\dots,Q_n(t))^{\top}$$, where $$Q_i(t) \simeq Q(z_i,t)$$, $$i = 1,2,\dots,n$$, and analogous for $$\boldsymbol{\theta}, \boldsymbol{w}$$.
For these new functions we now obtain the DAEs $$\begin{eqnarray} \boldsymbol{\dot{Q}} &=& -RaPra^2 \boldsymbol{\theta} - Pra^2 \boldsymbol{Q} + Pr (\boldsymbol{\underline{A}}_Q \boldsymbol{Q} + \boldsymbol{b}_Q),\\ \boldsymbol{\dot{\theta}} &=& \boldsymbol{w} - a^2 \boldsymbol{\theta} + \boldsymbol{\underline{A}}_{\theta} \boldsymbol{\theta} + \boldsymbol{b}_{\theta},\\ \boldsymbol{Q} &=& - a^2 \boldsymbol{w} + \boldsymbol{\underline{A}}_w \boldsymbol{w} + \boldsymbol{b}_w, \end{eqnarray}$$ with tridiagonal matrices $$\boldsymbol{\underline{A}}_Q, \boldsymbol{\underline{A}}_{\theta}, \boldsymbol{\underline{A}}_w$$ and with vectors $$\boldsymbol{b}_Q, \boldsymbol{b}_{\theta}, \boldsymbol{b}_w$$ which arise from the centered finite differences and from the boundary conditions on $$Q, \theta, w$$.
These are not ODEs because the time derivative of $$\boldsymbol{w}$$ is missing. However, because the third equation is linear, I would suggest to eliminate $$\boldsymbol{w} = (\boldsymbol{\underline{A}}_w - a^2 \boldsymbol{\underline{I}})^{-1}(\boldsymbol{Q} - \boldsymbol{b}_w)$$ using the third equation ($$\boldsymbol{\underline{I}}$$ denoting the identity matrix) and plug into the second equation to obtain an actual (linear) system of ODEs for $$\boldsymbol{Q}$$ and $$\boldsymbol{\theta}$$ only: $$\begin{eqnarray} \boldsymbol{\dot{Q}} &=& -RaPra^2 \boldsymbol{\theta} - Pra^2 \boldsymbol{Q} + Pr (\boldsymbol{\underline{A}}_Q \boldsymbol{Q} + \boldsymbol{b}_Q),\\ \boldsymbol{\dot{\theta}} &=& (\boldsymbol{\underline{A}}_w - a^2 \boldsymbol{\underline{I}})^{-1}(\boldsymbol{Q} - \boldsymbol{b}_w) - a^2 \boldsymbol{\theta} + \boldsymbol{\underline{A}}_{\theta} \boldsymbol{\theta} + \boldsymbol{b}_{\theta}. \end{eqnarray}$$ You can now use your favorite Runge-Kutta method to solve the system of ODEs for $$\boldsymbol{Q}, \boldsymbol{\theta}$$.