# Efficient Solution of Fokker-Planck Equations in Two Dimensions
## 5 Numerical Examples
### Linear Oscillator  
The FP equation for a linear system has an analytical solution [1; 2], and it is frequently used as a benchmark problem for testing numerical methods [4; 6; 8]. We consider the well-known SDOF linear oscillator (LO)  
\[g_{1}(x_{1},x_{2}) =x_{2}, \tag{30a}\] \[g_{2}(x_{1},x_{2}) =-2\zeta\omega_{0}x_{2}-\omega_{0}^{2}x_{1}, \tag{30b}\]  
with the parameters \(\zeta=0.05\) and \(\omega_{0}=1\), the excitation \(D=0.1\), and the initial conditions \((\mu_{10},\mu_{20})=(5,5)\) and \((\sigma_{10}^{2},\sigma_{20}^{2})=(1/9,1/9)\). The domain for computation is square, \(-10\leq x_{1},x_{2}\leq 10\). The FE mesh is uniform, with an element size of 0.2 in both the \(x_{1}\) and \(x_{2}\) directions. We begin with the FEA represented by Equations (8) and (9), applied to the state Equation (30). Using a composite solution vector \(\mathbf{u}\) composed of snapshots taken from a reference Runge-Kutta solution of the full FE equations, the coefficients resulting from the projection of this vector onto all the eigenvectors of \(\mathbf{M}\) are plotted in Figure 1, where they have been sorted from large to small. The norm of the error resulting when \(\mathbf{u}\) is approximated with an increasing number of eigenvectors follows a similar trend. When artificial diffusive terms were used in the solution of this problem on reduced bases, parameters were chosen after a limited amount of experimentation. In the \(x_{1}\) direction, the parameters \(b_{1}=2.4\) and \(d_{1}=1.6\) were selected, such that the added terms are nonzero only for \(|x_{1}|>7.6\) and rise to a maximum value of 1.6 on the boundary of \(\Omega\). Likewise, in the \(x_{2}\) direction, with \(b_{2}=2.5\) and \(d_{2}=1.3\), the boundary terms enter where \(|x_{2}|>7.5\), with coefficients that increase to 1.3 on \(\partial\Omega\). While an effort was made to minimize the effect of these terms on the original PDE (beyond stabilizing clearly erroneous solutions), the results obtained were not particularly sensitive to the values chosen. Figure 2 compares solutions computed using a reduced basis of 2000 eigenvectors, with and without artificial diffusive boundary terms, to the exact solution. The basis vectors were selected according to the projection coefficient magnitudes in Figure 1. Figure 2a shows the evolution with time of the nonstationary pdf at the origin of the phase plane. Both reduced solutions exhibit initial large errors because of the distortion of the initial conditions caused by the transformation onto the smaller basis. Both then track the exact solution until approximately \(\tau=6\), when the reduced solution computed without added boundary terms begins to visibly exceed the reference curve. This trend continues, and would eventually result in divergence rather than a stationary solution. In contrast, the numerical solution computed on the same basis but with diffusive terms introduced  
Figure 1: Coefficients resulting from the projection of the composite solution vector \(\mathbf{u}\) onto the eigenvectors of \(\mathbf{M}\) for the linear oscillator. near the boundary continues to agree with the exact analytical result, exceeding it only slightly as \(\tau\) approaches 20. A cross section of the stationary pdf along the displacement axis of the phase plane is plotted in Figure 2b, where it may be seen that the artificial diffusive boundary terms both improve the accuracy of the peak value and reduce spurious oscillations in the tails of the distribution. However, both reduced solutions do oscillate rather than going smoothly to zero in the tails, and consequently, they take on physically meaningless negative values at some points. Figures 3 and 4 reveal the effect of basis size, depicting the same responses discussed above, this time computed using a basis of 1000, 1500, or 2000 eigenvectors. All of the reduced problems were solved using artificial diffusive terms to improve the solutions.