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} $$


$$\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.

  • $\begingroup$ Your finite difference equation is explicit with respect to the time step, so why do you even need to work in terms of matrices? $\endgroup$
    – Chester Miller
    Oct 28, 2018 at 0:15
  • $\begingroup$ @ChesterMiller I trying to understand how to implement for the explicit from so that I can implement the implicit form later. $\endgroup$
    – Natasha
    Oct 28, 2018 at 2:03
  • $\begingroup$ Do you have access to the IMSL library of subroutines? $\endgroup$
    – Chester Miller
    Oct 28, 2018 at 2:27
  • $\begingroup$ @ChesterMiller Sorry, I haven't had a chance to use IMSL before. I work on MATLAB.I'm trying to formulate the model on my own to get a better understanding on how the system works. $\endgroup$
    – Natasha
    Oct 28, 2018 at 3:00
  • 1
    $\begingroup$ I'm voting to close this question as off topic because it is primarily a mathematics question, even if there could be physics application. $\endgroup$
    – Aaron Stevens
    Oct 28, 2018 at 3:25

1 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?

  • $\begingroup$ Could you please have a look at the updated post? $\endgroup$
    – Natasha
    Oct 28, 2018 at 13:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.