# Efficient Solution of Fokker-Planck Equations in Two Dimensions
## 6 Results  
We have investigated two techniques to improve the speed of solution of Fokker-Planck equations as discretized using finite element analysis in the phase plane. In the first of these, reduced bases were constructed from selected eigenvectors of the coefficient matrices \(\mathbf{M}\) and \(\mathbf{K}\) appearing in the resulting first-order differential equation. Basis vectors obtained by solving a regular eigenproblem in \(\mathbf{M}\) were attractive because they were real-valued, could be used to diagonalize \(\mathbf{M}\), and were independent of the parameters of the underlying dynamical system and the intensity of the excitation. The size of the basis required in order for the reduced problem to capture well the dynamics of the full problem was typically one-fifth to one-half that of the original system. Considering that the number of elements of the reduced coefficient matrices goes roughly as the square of this ratio, this represents an attractive savings in storage and computation; however, it does not always enable the drastic improvement in computational efficiency sometimes associated with a reduced-order model. In addition, even when a large basis was used, the reduced problem was found to be susceptible to spurious oscillations in the tails of the computed pdf, and to loss of probability mass over the edges of the computational domain. These phenomena were effectively addressed by the introduction into the original PDE of artificial diffusive terms near the boundaries of the domain. While greater gains were achieved in the linear example, in the more representative, nonlinear Duffing oscillator example, the typical reduction in run time, compared to integration by a conventional adaptive Runge-Kutta algorithm, was approximately 80%. This is certainly significant; however, we must reiterate that the reduced bases were selected with the benefit of the composite solution vector \(\mathbf{u}\), and that even a modest reduction in the size of the basis used was found to require the introduction of artificial diffusive terms to stabilize the solution near the boundaries of the \begin{table}
\begin{tabular}{l r} \hline \hline
**Method** & **Time, \%** \\ \hline RK (9801) & 100.00 \\ RK (4500) & 19.84 \\ RK (4500) and ABT & 19.76 \\
1 st-order splitting (9801) & 1.57 \\
2 nd-order splitting (9801) & 1.61 \\
1 st-order splitting (4500) & 1.00 \\
1 st-order splitting (4500) and ABT & 1.03 \\ \hline \hline \end{tabular}
\end{table} 
Table 2: Integration times for various methods of solution for the Duffing oscillator (\(\tau\) = 0:0.01:4). The size of the basis used in each method is shown in parentheses. ABT denotes the use of artificial boundary terms. Run times are normalized with respect to the Runge–Kutta solution of the full FE problem. computational domain. In view of these limitations, we cannot enthusiastically recommend the use of this reduction. The second numerical method applied herein to the same discretized PDE was operator splitting, where the singular coefficient matrix \(\mathbf{K}\) was written as the sum of two nonsingular matrices. This allowed the rapid computation of the joint response pdf with a uniform time step using the constant state transition matrix given by the product of the exponentials of the two component matrices. Both first- and second-order-accurate splitting schemes were considered. It was observed that, although the cost of second-order methods was not much greater, first-order splitting was usually sufficiently accurate while being simpler to program, and that separating the FE coefficient matrix \(\mathbf{K}\) according to the physical source of the terms in the original PDE (i.e., convection or diffusion) was effective in the linear and Duffing oscillator examples studied. The resultant saving in computational effort was comparable to that obtained with basis reduction for the linear oscillator example, a little more than 80%, but savings of over 98% were observed when operator splitting was used to compute the joint pdf of a Duffing oscillator.