# Coupled Partial Differential Equations

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.

• I suggest taking a look at the FiPy PDE solver (written in Python) rather than trying to write your own solver. Commented Jul 5 at 18:52
• 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? Commented Jul 7 at 7:27

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,