Numerically solving a partial differential equation

Question

I am trying to numerically solve the following PDE,

$$\frac{\partial u^A}{\partial t} = c_1\frac{\partial^2 u^A}{\partial^2x} \,,$$

where $c_1$ is a constant. The above can be discretized using the numerical approximations $$ \begin{align} \frac{\partial u^A}{\partial t} & ~\Rightarrow~ \frac{u_{i}^{t+1} - u_i^t}{\Delta t } \\[5px] \frac{\partial^2 u^A}{\partial^2x} &~\Rightarrow~ \frac{u_{i+1}^t -2u_i^t+ u_{i-1}^t}{\Delta x^2 } \end{align} $$

Gives,

$$\frac{u_{i}^{t+1} - u_i^t}{\Delta t} = c_1\frac{u_{i+1}^t -2u_i^t+ u_{i-1}^t}{\Delta x^2 }$$

The above can be rearranged,

$$ \begin{align} u_{i}^{t+1} &= u_{i}^t + \frac{c_1 \Delta t}{\Delta x}(u_{i+1}^t -2u_i^t+ u_{i-1}^t) \\[5px] u_{i+1}^{t+1} &= u_{i+1}^t + \frac{c_1 \Delta t}{\Delta x}(u_{i+2}^t -2u_{i+1}^t+ u_{i}^t) \end{align} $$ and so on.

In matrix notation, $$U^{t+1} = [u_{i}^{t+1},u_{i+1}^{t+1}, \, \dots ,u_{i+n}^{t+1}]$$

Therefore, $$ U_A^{t+1} = U_A^{t} + c_1* \text{tridiagonal matrix} * U_A^t $$

The above can be solved to obtain $U_A^{t+1}$.

Now when there is another species B, the equation will be $U_B^{t+1} = U_B^{t} + c_2* \text{tridiagonal matrix} * U_B^t$

I am solving for the time evolution of function $u$ of species A and B as two separate matrix equations. Could someone suggest if there is a way to formulate this as a single matrix equation? I'm trying to create a diagonal matrix of the constants c and then combine the matrices.

Any suggestions?

Edit: I'm trying an alternate approach. If I were to use the method of lines the equations would be

$$\frac{du^A}{dt} = c_1\frac{u_{i+1}^{t^A} -2u_i^{t^A}+ u_{i-1}^{t^A}}{\Delta x^2 }$$ $$\frac{du^B}{dt} = c_2\frac{u_{i+1}^{t^B} -2u_i^{t^B}+ u_{i-1}^{t^B}}{\Delta x^2 }$$

In the matrix form, $[\frac{du^A}{dt} \frac{du^B}{dt}]^T = \frac{1}{\Delta x^2}*diagonalmatrix*[u_{i+1}^{t^A} -2u_i^{t^A}+ u_{i-1}^{t^A} ; u_{i+1}^{t^B} -2u_i^{t^B}+ u_{i-1}^{t^B} ]$ ,

where diagonalmatrix contains the constants $c_1$ and $c_2$ along its diagonal.

Would this be the right way to proceed? I have the initial condition and boundary conditions.

Answer 1

Assuming that the time-step taken between the two species is the same, you should be able to just simply combine the two species into a single vector: $$ U_\text{tot}=\left[\begin{array}{c}u_0^A \\ u_1^A\\ \vdots \\u_{n+1}^A \\ u_0^B \\ u_1^B \\ \vdots \\ u_{n+1}^B \end{array}\right] $$ And turn the two individual matrix operators into the single matrix, $$ \mathsf{X}_\text{tot}=\left(\begin{array}{cc}\mathsf{X}_A & \mathbf 0 \\ \mathbf 0 & \mathsf{X}_B\end{array}\right) $$ where $\mathsf{X}_A =c_1 * \text{tridiagonal matrix}$ and $\mathbf 0$ a matrix of 0's. Thus you still have a tridiagonal problem to solve: $$U^{n+1}_\text{tot}=\mathsf{X}_\text{tot}U^n_\text{tot}$$

On a programming aspect, there probably won't be much of a difference in the timing between doing them separately and together like the above. Maybe a collective few microseconds due to the duplicate calls to the same update function with different vectors and matrices, but you'll probably lose those gains with populating & de-populating $U_\text{tot}$ and $\mathsf {X}_\text{tot}$. Perhaps the combined method would allow you to add in some cross terms (in the operator), so that may be one advantage to running with this method?