# Numerically solving a partial differential equation

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.

• Your finite difference equation is explicit with respect to the time step, so why do you even need to work in terms of matrices? – Chester Miller Oct 28 '18 at 0:15
• @ChesterMiller I trying to understand how to implement for the explicit from so that I can implement the implicit form later. – Natasha Oct 28 '18 at 2:03
• Do you have access to the IMSL library of subroutines? – Chester Miller Oct 28 '18 at 2:27
• @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. – Natasha Oct 28 '18 at 3:00
• I'm voting to close this question as off topic because it is primarily a mathematics question, even if there could be physics application. – Aaron Stevens Oct 28 '18 at 3:25

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?