Crank-Nicolson for coupled PDE's

Question

$\newcommand{\T}{T}$ $\newcommand{\partiald}[2]{\frac{\partial #1}{\partial #2}}$ $\newcommand{\partialdd}[2]{\frac{\partial^2 #1}{\partial #2^2}}$

I am trying to solve a set of coupled PDE's with the Crank-Nicolson method. So far I have used it to solve a single PDE, the 1D diffusion problem in the Wikipedia article I have linked. This involved turning the set of equations into a matrix equation

$$ \begin{bmatrix} A & 1 & 0 & \cdots & & 0 \\ 1 & A & 1 & \cdots & & \\ 0 & 1 & A & \ddots & & \\ \vdots & \vdots & \ddots & \ddots \\ & & & & A & 1 \\ 0 & & & & 1 & A \end{bmatrix} \begin{bmatrix} \T^{n+1}_1\\ \T^{n+1}_2\\ \T^{n+1}_3\\ \vdots \\ \T^{n+1}_{N-2}\\ \T^{n+1}_{N-1} \end{bmatrix} = \begin{bmatrix} d_1^n - \T_0^{n+1}\\ d_2^n \\ d_3^n \\ \vdots \\ d_{N-2}^n\\ d_{N-1}^n - \T_{N}^{n+1} \end{bmatrix} $$ where $A = -\frac{1+2r}{r}$ and $$d_i^n = - \left[\T_{i+1}^{n} + \left(\frac{1-2r}{r}\right) \T_i^{n} + \T_{i-1}^{n}\right]$$ I specify $T_0^{n+1}$ and $T_N^{n+1}$ as boundary conditions, and provide an initial set of values for $T_i^0$. From there it is simply a matter of iteratively solving the set of linear equations for the $T_i^{n+1}$.

Now I want to apply this method to solve the set of coupled equations \begin{align*} \partiald{v}{t} &= - \partiald{\rho}{x} \\ \partiald{\rho}{t} &= - \partiald{v}{x} \end{align*}

Applying the Crank-Nicolson algorithm I find the following relation \begin{align*} v_i^{n+1} + r\left(\rho_{i+1}^{n+1} - \rho_{i-1}^{n+1}\right) &= v_i^n - r\left(\rho_{i+1}^n - \rho_{i-1}^n\right) \\ \rho_i^{n+1} + r\left(v_{i+1}^{n+1} - v_{i-1}^{n+1}\right) &= \rho_i^n - r\left(v_{i+1}^n - v_{i-1}^n\right) \end{align*} where $r \equiv \Delta t/ 4\Delta x$. But from here I am stuck, I don't know how to turn this into a matrix equation due to having two variables. I am assuming I want to turn this into a vector equation, like the way I did in this answer, but I'm not sure how in this case I would then convert it into a matrix equation which can be iteratively solved.

Edit:

I realize that I can convert these into two separte equations \begin{align*} \partialdd{\rho}{x} &= \partialdd{\rho}{t}\\ \partialdd{v}{x} &= \partialdd{v}{t} \end{align*} which I think I can solve using the previous method, though I think I will need to use an explicit algorithm for the first time step. But my question above still stands as to how to solve the set of coupled equations.

Answer 1

$\renewcommand{\d}{\vec{d}} \renewcommand{\S}{\vec{S}} \renewcommand{\v}{v} \renewcommand{\r}{\rho}$

I solved the problem by doing the following: starting from

\begin{align*} \v_i^{n+1} + r\left(\r_{i+1}^{n+1} - \r_{i-1}^{n+1}\right) &= \v_i^n - r\left(\r_{i+1}^n - \r_{i-1}^n\right) \\ \r_i^{n+1} + r\left(\v_{i+1}^{n+1} - \v_{i-1}^{n+1}\right) &= \r_i^n - r\left(\v_{i+1}^n - \v_{i-1}^n\right) \end{align*}

Define a vector $\vec{S}_i^n = (\tilde{v}_i^n , \r_i^n)^{\mathrm T}$. Then we can write these two equations as $$ \S_i^{n+1} + \begin{bmatrix} 0 & r \\ r & 0 \end{bmatrix} \left(\S_{i+1}^{n+1} - \S_{i-1}^{n+1}\right) = \S_i^n - \begin{bmatrix} 0 & r \\ r & 0 \end{bmatrix} \left(\S_{i+1}^n - \S_{i-1}^n\right) $$ call this matrix $A$. Then we write our set of equations as $$ \begin{bmatrix} 1 & A & & & 0 \\ -A & 1 & \ddots & & \\ & \ddots & \ddots & & \\ & & & 1 & A \\ 0 & & & -A & 1 \end{bmatrix} \begin{bmatrix} \S_1^{n+1} \\ \S_2^{n+1} \\ \vdots \\ \S_{N-2}^{n+1} \\ \S_{N-1}^{n+1} \end{bmatrix} = \begin{bmatrix} 1 & -A & & \! & 0 \\ A & 1 & \ddots & \! & \\ & \ddots & \ddots & \! & \\ & & & \! 1 & -A \\ 0 & & & \! A & 1 \\ \end{bmatrix} \begin{bmatrix} \S_1^{n} \\ \S_2^{n} \\ \vdots \\ \S_{N-2}^{n} \\ \S_{N-1}^{n} \\ \end{bmatrix} + \begin{bmatrix} A\left(\S_0^{n} + \S_0^{n+1}\right) \\ 0 \\ \vdots \\ 0 \\ -A\left(\S_N^{n} + \S_N^{n+1}\right) \\ \end{bmatrix} $$