# Efficient Solution of Fokker-Planck Equations in Two Dimensions
## 5 Numerical Examples
### Duffing Oscillator  
The governing equation of the nonlinear Duffing oscillator is, in state form,  
\[g_{1}(x_{1},x_{2}) =x_{2}, \tag{31}\] \[g_{2}(x_{1},x_{2}) =-2\zeta\omega_{0}x_{2}-\omega_{0}^{2}x_{1}(\gamma+\varepsilon x_{ 1}^{2}), \tag{32}\]  
with \(\zeta=0.2\), \(\omega_{0}=1\), \(\gamma=-1\), \(\epsilon=0.1\), and \(D=0.4\). (Like the previous example, this system is taken from [6].) These parameters are those of a system with negative linear stiffness but positive cubic stiffness, resulting in a bimodal stationary pdf that is doubly symmetric about the origin of the phase plane and approximately Gaussian in \(x_{2}\). The nonlinearity thus controls the stationary response despite the small numerical value of \(\epsilon\). The initial conditions were taken to be \((\mu_{10},\mu_{20})=(0,10)\) and \((\sigma_{10}^{2},\sigma_{20}^{2})=(1/2,1/2)\), and the FE domain was \(-15\leq x_{1}\leq 15\), \(-20\leq x_{2}\leq 20\), with 100 elements used in each dimension. Guided by the projection coefficients of a composite vector \(\mathbf{u}\) and some experimentation, we determined that somewhat more than 4000 basis vectors were needed to capture the dynamics of this problem in reduced form. The diffusive boundary layer coefficients \(b_{1}=8\) and \(d_{1}=4.6\) and \(b_{2}=13\) and \(d_{2}=0.1\) were then chosen for the \(x_{1}\) and \(x_{2}\) directions, respectively. The discretized equation for the joint pdf of the Duffing oscillator was reduced onto a basis of 4500 eigenvectors, with and without introducing artificial boundary terms in the PDE. The nonstationary response in Figure 5a is shown for the point \((3.3,0)\), which is the approximate location in the phase plane of one of the two peaks of the bimodal stationary distribution. Both the full and reduced problems were integrated here by a standard Runge-Kutta routine. When artificial diffusive boundary terms are included in the formulation, the solution of the resulting reduced problem approximates well the solution of the full problem, after a transient error attributed to the inability of the reduced basis to represent the given initial conditions accurately. Without added boundary terms, however, the solution does not approach the correct stationary value, and begins to decrease from it as probability mass is lost over the edges of the FE domain. Figure 5b, a section of the stationary pdf along the \(x_{1}\) axis, shows spurious oscillations in the tails of the pdf. (The exact stationary solution shown there is from Caughey [3].) Both reduced problems underestimate the pdf in the valley between the peaks of the stationary distribution. While \begin{table}
\begin{tabular}{l r} \hline \hline
**Method** & **Time, \%** \\ \hline RK (9801) & 100.00 \\ RK (3000) & 15.90 \\ RK (2000) & 11.98 \\ RK (3000) and ABT & 17.17 \\ RK (2000) and ABT & 12.70 \\
1 st-order splitting (9801) & 17.36 \\
1 st-order splitting (3000) & 10.07 \\
1 st-order splitting (2000) & 9.69 \\
1 st-order splitting (3000) and ABT & 10.16 \\
1 st-order splitting (2000) and ABT & 9.81 \\ \hline \hline \end{tabular}
\end{table} 
Table 1: Integration times for various methods of solution for the linear 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. the inclusion of diffusive boundary terms reduces some of the extraneous ripples, it must be said that these results are not as clean as those obtained for the linear oscillator, despite the use of a much larger reduced basis (here, almost 46% of the size of the full problem). When the results of first- and second-order splitting are compared to those of a Runge-Kutta solution of the unsplit problem, second-order splitting is found to be more accurate, but both solutions from splitting are found to agree quite well with the RK solution during the entire response, including at stationarity.