# Efficient Solution of Fokker-Planck Equations in Two Dimensions
## 2 Problem Formulation  
Consider an SDOF system (an oscillator) governed by the stochastic differential equations\[\dot{X}_{1} =g_{1}(\mathbf{X}), \tag{1a}\] \[\dot{X}_{2} =g_{2}(\mathbf{X})+w(t), \tag{1b}\]  
where the state vector is \(\mathbf{X}= \begin{bmatrix}X&\dot{X}\end{bmatrix} ^{\mathrm{T}}= \begin{bmatrix} X_{1}&X_{2}\end{bmatrix} ^{\mathrm{T}}\) and the overdot denotes differentiation with respect to time. The excitation \(w(t)\) is Gaussian white noise, which is fully defined by its first two moments,  
\[E[w(t)]=0\qquad\text{and}\qquad E[w(t)w(t+t^{\prime})]=2D\delta(t^{\prime}), \tag{2}\]  
where \(\delta(\cdot)\) is the Dirac delta function and \(D/\pi\) is the magnitude of the constant two-sided spectral density function of the excitation. Equation (1) describes a Markov process, and its behavior is completely determined by the transition probability density function \(p(\mathbf{x},t|\mathbf{x}_{0})\), which satisfies the (deterministic) Fokker-Planck equation  
\[\frac{\partial p}{\partial t}=D\frac{\partial^{2}p}{\partial x_{2}^{2}}-\frac {\partial}{\partial x_{1}}[g_{1}(\mathbf{x})p]-\frac{\partial}{\partial x_{2 }}[g_{2}(\mathbf{x})p] \tag{3}\]  
and the initial condition \(\lim_{t\to 0}p(\mathbf{x},t|\mathbf{x}_{0})=\delta(\mathbf{x}-\mathbf{x}_{0})\). Homogeneous boundary conditions are imposed where \(p\to 0\) as \(|x_{1}|\), \(|x_{2}|\to\infty\). The joint probability density function of the response is related to the transition probability density by  
\[f(\mathbf{x},t)=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}p(\mathbf{x},t| \mathbf{x}_{0})f(\mathbf{x}_{0})\,\mathrm{d}\mathbf{x}_{0}, \tag{4}\]  
where for convenience the initial condition may be taken in the form of a binormal distribution. Substituting this function into the FP equation and integrating term-wise gives  
\[\frac{\partial f}{\partial t}=D\frac{\partial^{2}f}{\partial x_{2}^{2}}-\frac {\partial}{\partial x_{1}}(g_{1}f)-\frac{\partial}{\partial x_{2}}(g_{2}f). \tag{5}\]  
The desired joint probability density function \(f(\mathbf{x},t)\) satisfies this PDE subject to the initial conditions on \(f\) and the boundary conditions (which are the same as those on \(p\)). A Galerkin finite element method is used where within a single element with domain \(\Omega_{e}\), the probability density function is interpolated according to  
\[f(\mathbf{x},t)=\sum_{r=1}^{n}N_{r}(\mathbf{x})f_{r}^{e}(t), \tag{6}\]  
where the \(N_{r}(\mathbf{x})\), \(r=1,2,\ldots,n\), are bilinear trial functions of class \(C^{0}\), \(n\) is the number of nodes in a single element, and \(f_{r}^{e}\) are the nodal values of the density. Substituting Equation (6) into the weak form of Equation (5) leads to  
\[\sum_{s=1}^{n}m_{rs}f_{s}^{e}+\sum_{s=1}^{n}k_{rs}f_{s}^{e}=0,\qquad r=1,\ldots,n, \tag{7}\]  
where the element matrices are given by  
\[m_{rs}=\iint_{\Omega_{e}}N_{r}N_{s}\,\mathrm{d}x_{1}\,\mathrm{d}x_{2}\] (8a) and \[k_{rs}=\iint_{\Omega_{e}}\left[D\frac{\partial N_{r}}{\partial x_{2}}\frac{ \partial N_{s}}{\partial x_{2}}\right.\left.-N_{s}\left(g_{1}\frac{\partial N _{r}}{\partial x_{1}}+g_{2}\frac{\partial N_{r}}{\partial x_{2}}\right) \right]\,\mathrm{d}x_{1}\,\mathrm{d}x_{2}. \tag{8b}\]Assembling the global coefficient matrices by standard procedures produces  
\[\mathbf{M}\mathbf{f}(t)+\mathbf{K}\mathbf{f}(t)=\mathbf{0}, \tag{9}\]  
where \(\mathbf{f}(t)\) is the global vector of nodal values of the joint probability density function and \(\mathbf{M}\) and \(\mathbf{K}\) are constant, square matrices. We note some properties of the coefficient matrices that are relevant to the discussion to follow. Both \(\mathbf{M}\) and \(\mathbf{K}\) are real and sparse.