I'm trying to solve the following system of coupled differential equations, the two-temperature model for $e$ = electrons and $l$ = lattice.

$$ \rho_{e}C_{p,e}\frac{\partial T_{e}}{\partial t} = k_{e}\frac{\partial^{2}T_{e}}{\partial x^{2}} - G(T_{e}-T_{l}) + S(x,t) $$ $$ \rho_{l}C_{p,l}\frac{\partial T_{l}}{\partial t} = k_{l}\frac{\partial^{2}T_{l}}{\partial x^{2}} + G(T_{e}-T_{l}) $$

I was thinking of re-casting them into a form suitable for the Crank-Nicolson method, which would be the following if $T_{e}=u_{i}^{n}$ and $T_{l}=v_{i}^{n}$ where $i$ is the space node and $n$ is the time node.

$$ A_{o}u_{i+1}^{n+1} + B_{o}u_{i}^{n+1} + A_{o}u_{i-1}^{n+1} = A_{1}u_{i+1}^{n} + B_{1}u_{i}^{n} + A_{1}u_{i-1}^{n+1} + C_{1}v_{i}^{n} + D_{1}S_{i}^{n} $$ $$ A_{2}v_{i+1}^{n+1} + B_{2}v_{i}^{n+1} + A_{2}v_{i-1}^{n+1} = A_{3}v_{i+1}^{n} + B_{3}u_{i}^{n} + A_{3}v_{i-1}^{n+1} + C_{3}v_{i}^{n} $$

where $A_{o}=-\frac{dt}{\rho_{e}C_{p,e}}\frac{k_{e}}{2dx^{2}}, \ A_{1} = -A_{o}, $ $B_{o}=1+\frac{dt}{\rho_{e}C_{p,e}}\frac{k_{e}}{dx^{2}}, \ B_{1} = 1-\frac{dt}{\rho_{e}C_{p,e}}\frac{k_{e}}{dx^{2}} + G, \ C_{1} = \frac{dt}{\rho_{e}C_{p,e}}G, \ D_{1} = \frac{dt}{\rho_{e}C_{p,e}}$

$A_{2}=-\frac{dt}{\rho_{l}C_{p,l}}\frac{k_{l}}{2dx^{2}}, \ A_{3}=-A_{2}$ $B_{2} = 1+\frac{dt}{\rho_{l}C_{p,l}}\frac{k_{l}}{2dx^{2}}, \ B_{3}=1-\frac{dt}{rho_{e}Cp_{e}}\frac{k_{l}}{2dx^{2}}+G, \ C_{3} = -\frac{dt}{\rho_{l}C_{p,l}}G$

Analytically, this works out but I'm drawing a blank on how to implement it (I want to use Python). I can do each of the equation systems by themselves easily enough, but I was trying to solve both in the same matrix system. I've seen certain attempts at this sort of problem by converting the two unknowns at the updated time into a vector $\vec{W}_{i}^{n+1}=(u_{i}^{n+1},v_{i}^{n+1})^{T}$ but due to the coefficients being different for each (apart from the coupling constant G) I wasn't sure how to go about it.

  • 1
    $\begingroup$ I suggest taking a look at the FiPy PDE solver (written in Python) rather than trying to write your own solver. $\endgroup$ Commented Jul 5 at 18:52
  • 3
    $\begingroup$ Why do you want to put both in the same system matrix? What benefit does this have? To be sure, it is easy enough to do, I just don't think you gain anything. If you have $Mu^{k+1} = b(u^k, v^k)$ and $Nv^{k+1} = c(u^k,v^k)$ then: $$\begin{bmatrix} b(u^k,v^k)\\ c(u^k, v^k)\end{bmatrix} = \begin{bmatrix} M & 0 \\ 0 & N\end{bmatrix}\begin{bmatrix} u^{k+1} \\ v^{k+1}\end{bmatrix}.$$ For an iterative solver you anyways can avoid constructing the matrix explicitly, and for a direct one why would you grow your problem size? $\endgroup$
    – lightxbulb
    Commented Jul 7 at 7:27

1 Answer 1


Good afternoon,

If the coefficients are correctly found (I have not checked), you can present the system as @lightxbulb says. Being $N$ the total number of space nodes, let me expand the answer a little bit for completeness, therefore

$u^{n+1} = \left[u_1^{n+1},\,u_2^{n+1},\,...,\,u_N\right]^\top$

$v^{n+1} = \left[v_1^{n+1},\,v_2^{n+1},\,...,\,v_N\right]^\top$

$b(u^n,v^n) = \left[b_1,\,b_2,\,...,b_N\right]^\top$

where $b_i,\forall i=2,...,N-1$ can be computed as Assuming $A_1u_{i-1}^{n+1}$ to be a typo in your question and be $A_1u_{i-1}^{n}$ instead, which makes more sense to me

$b_i = A_{1}u_{i+1}^{n} + B_{1}u_{i}^{n} + A_{1}u_{i-1}^{n} + C_{1}v_{i}^{n} + D_{1}S_{i}^{n}$

and $b_1$ and $b_N$ will depend on the boundary conditions of the problem.

$c(u^n,v^n) = \left[c_1,\,c_2,\,...,c_N\right]^\top$

where $c_i,\forall i=2,...,N-1$ can be computed as

$c_i = A_{3}v_{i+1}^{n} + B_{3}u_{i}^{n} + A_{3}v_{i-1}^{n+1} + C_{3}v_{i}^{n}$

and $c_1$ and $c_N$ will depend on the boundary conditions of the problem.

Finally, we have that

$M = \begin{bmatrix} B_0&A_0&0&\ldots & 0 & 0& 0\\ A_0 &B_0&A_0&\ldots& 0 & 0 & 0\\ \vdots &\vdots&\vdots&\ddots & \vdots & \vdots & \vdots\\ 0 &0&0&\ldots& A_0 &B_0&A_0\\ 0 &0&0&\ldots& 0 &A_0&B_0\\ \end{bmatrix} $

$N = \begin{bmatrix} B_2&A_2&0&\ldots & 0 & 0& 0\\ A_2 &B_2&A_2&\ldots& 0 & 0 & 0\\ \vdots &\vdots&\vdots&\ddots & \vdots & \vdots & \vdots\\ 0 &0&0&\ldots& A_2 &B_2&A_2\\ 0 &0&0&\ldots& 0 &A_2&B_2\\ \end{bmatrix} $

As @lightxbulb says, it becomes a decoupled system maybe you don't want to threat the system matrix as a whole. You can use a huge variety of methods for solving this linear system to obtain the vector of unkowns you mentioned

$\vec{W}^{n+1} = \begin{bmatrix} u^{n+1}\\v^{n+1} \end{bmatrix}$

I hope it helps,


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.