# Efficient Solution of Fokker-Planck Equations in Two Dimensions
## 4 Operator Splitting  
If the convection and diffusion terms in the original equation are separated, perhaps in the expectation that their physics may develop on different time scales, the terms in the element \(\mathbf{k}\) matrix can be grouped and that matrix can be divided into  
\[\mathbf{k}_{1}=\iint_{\Omega_{x}}D\frac{\partial\mathbf{N}}{\partial x_{2}} \frac{\partial\mathbf{N}^{\mathrm{T}}}{\partial x_{2}}\,\mathrm{d}x_{1}\, \mathrm{d}x_{2}\] (22a) and \[\mathbf{k}_{2}=-\iint_{\Omega_{x}}\left(g_{1}\frac{\partial\mathbf{N}}{\partial x _{1}}+g_{2}\frac{\partial\mathbf{N}}{\partial x_{2}}\right)\mathbf{N}^{ \mathrm{T}}\,\mathrm{d}x_{1}\,\mathrm{d}x_{2}. \tag{22b}\]Obviously, this splitting satisfies \(\mathbf{k}=\mathbf{k}_{1}+\mathbf{k}_{2}\). Then, global matrices \(\mathbf{K}_{1}\) and \(\mathbf{K}_{2}\) can be assembled as the previous \(\mathbf{K}\) matrix was and transformed to obtain \(\tilde{\mathbf{K}}_{1}=\mathbf{S}^{\mathsf{T}}\mathbf{K}_{1}\mathbf{S}\) and \(\tilde{\mathbf{K}}_{2}=\mathbf{S}^{\mathsf{T}}\mathbf{K}_{2}\mathbf{S}\). For the system of ordinary differential equations of the form \(\dot{\boldsymbol{\eta}}+\tilde{\mathbf{K}}\boldsymbol{\eta}=\mathbf{0}\) obtained using \(\mathbf{S}= \begin{bmatrix}\boldsymbol{\phi}_{1}&\dots&\boldsymbol{\phi}_{n} \end{bmatrix} \) (so that \(\hat{\mathbf{M}}=\mathbf{I}\)), we begin with the simple first-order splitting method and  
\[\dot{\boldsymbol{\eta}}+(\tilde{\mathbf{K}}_{1}+\tilde{\mathbf{K}}_{2}) \boldsymbol{\eta}=\mathbf{0}, \tag{23}\]  
and separate the original equation into two parts. Over a time step \(\Delta_{t}\), we solve the two equations  
\[\dot{\mathbf{q}}+\tilde{\mathbf{K}}_{1}\mathbf{q}=\mathbf{0},\quad\dot{ \mathbf{p}}+\tilde{\mathbf{K}}_{2}\mathbf{p}=\mathbf{0} \tag{24}\]  
separately, using the solution of the first as the initial conditions for the second, and then combine the two separate solutions to form a solution to the original Equation (23) [35]. Using the state transition matrices associated with Equation (24), we can advance the complete solution from time step \(n\) to step \(n+1\) by computing  
\[\mathbf{q}_{n+1}=\mathrm{e}^{-\mathbf{K}_{1}\Delta_{t}}\mathbf{p}_{n},\quad \mathbf{p}_{n+1}=\mathrm{e}^{-\mathbf{K}_{2}\Delta_{t}}\mathbf{q}_{n+1}. \tag{25}\]  
Combining these equations for a fixed time step \(\Delta_{t}\), the result  
\[\boldsymbol{\eta}_{n+1}=\mathrm{e}^{-\tilde{\mathbf{K}}_{2}\Delta_{t}}\mathrm{ e}^{-\tilde{\mathbf{K}}_{1}\Delta_{t}}\boldsymbol{\eta}_{n} \tag{26}\]  
is obtained, which may be used to efficiently advance the solution over any number of steps by simply multiplying the previous value by a constant matrix. By comparison to the Taylor series for the exact solution of a canonical equation, it can be shown that this splitting method is first-order accurate (\(\mathcal{O}(\Delta_{t})\)) even when the matrices \(\mathbf{K}_{1}\) and \(\mathbf{K}_{2}\) do not commute under multiplication [49]. This scheme is both simple and stable once the matrix exponentials are computed and multiplied to form the single-state transition matrix needed to advance the solution over a fixed step. It is of course possible to choose a time step so large that the exponentiation (carried out in floating-point arithmetic) fails, but our experience has been that using time steps small enough to resolve the dynamics of interest in a system usually avoids any such difficulties. Symmetric splitting, as introduced by [PERSON] [35], can achieve second-order accuracy. Here, \(\mathbf{K}\) is split into  
\[\mathbf{K}=\left(\tfrac{1}{2}\mathbf{K}_{1}+\mathbf{K}_{2}\right)+\tfrac{1}{ 2}\mathbf{K}_{1}. \tag{27}\]  
Examining the Taylor series shows that  
\[\mathrm{e}^{\frac{1}{2}\mathbf{K}_{1}\Delta_{t}}\mathrm{e}^{\mathbf{K}_{2}} \mathrm{e}^{\frac{1}{2}\mathbf{K}_{1}\Delta_{t}}-\mathrm{e}^{(\mathbf{K}_{1}+ \mathbf{K}_{2})\Delta_{t}}=\mathbf{C}\Delta_{t}^{3}+\mathcal{O}(\Delta_{t}^{4}), \tag{28}\]  
where  
\[\mathbf{C}= \begin{bmatrix}[\mathbf{K}_{1},\mathbf{K}_{2}],\mathbf{K}_{1} \end{bmatrix} +2