# Moment-Based Stochastic Analysis of a Bistable Energy Harvester with Matching Network
## 6 Nonlinear System Analysis
\right] \tag{34}\]  
Substituting (30) and (34) into (28) and (31), respectively, we obtain a system of nonlinear differential algebraic equations (DAEs) that can be solved numerically to determine the admissible values for \(\boldsymbol{\mu}\) and \(\overline{\boldsymbol{\sigma}}\), from which the covariance matrix \(\boldsymbol{\sigma}\) is readily derived. Figure 4 shows the root mean square output voltage for the nonlinear bistable energy harvester, as a function of the matching network parameters \(L_{S}\) and \(C_{P}\). The elastic constant model parameters are set to the values \(k_{1}=5.4046\) kN/m and \(k_{3}=5.4046\) TN/m. The latter is chosen to be very large to achieve a reasonably small value for the displacement \(x\). The other parameters are the same as those of Table 1. The output voltage achieves the maximum value of \(v_{o(\text{rms})}^{\text{max}}=20.59\) V for \(L_{S}^{\text{(opt)}}=11.6\) H and \(C_{P}^{\text{(opt)}}=8.63\) nF. Interestingly, the optimal inductance value is a little lower for the nonlinear system with respect to the linear case. Note that since the noise intensity is very small, the figure resembles that of the linear system. However, since the equilibrium position is not the vertical resting state, the optimal values of the parameters are different. Figure 5 shows the root mean square output voltage as a function of the noise intensity. Red squares refer to the energy harvester with matched load with optimal values of the circuit components given above. Blue circles refer to the energy harvester with a simple resistive load. The output voltage was obtained from the numerical integration of the SDEs system using the stochastic Runge-Kutta method. For each value of the noise intensity, we averaged over 20 simulations, each one with a different realization of the noise process, with time length \(\Delta T=10^{4}\) s and time integration step \(\delta t\approx 50\) us. Root mean square values were calculated averaging over all simulations, after removing the initial transients. We notice that although the matching network optimization is carried out for a small noise intensity, the matched load setup offers higher output voltage, and therefore higher output average power and better power efficiency, also for a relatively high noise intensity. Figure 6 shows the displacement \(x\) versus time for two values of the noise intensity. On the left, for \(\varepsilon=10^{-3}\), the beam-inertial mass system remains confined in the right potential well, vibrating around the equilibrium position with small amplitude. On the right, for \(\varepsilon=5\cdot 10^{-3}\), the beam-inertial mass system jumps irregularly back and forth from one well to the other. In all our simulations, the initial condition was set to \(\mathbf{X}(0)=0\), with the  
To validate the accuracy of the method based on Gaussian moment closure, we compare theoretical predictions with results from numerical simulations. Figure 7 shows the relative error between theoretical predictions obtained using the Gaussian moment closure method and numerical simulations changing the noise intensity, for both the root mean square output voltage (red squares) and the power efficiency (blue circles). The relative error is evaluated as  
\[\Delta x=\left|\frac{x_{\text{theoretical}}-x_{\text{numerical}}}{x_{\text{numerical }}}\right| \tag{35}\]  
As expected, the relative error is very small for weak noise, and it increases along with the noise intensity. The error becomes rapidly significant when the beam begins to jump between the two potential wells because the corresponding probability density function is no longer a multivariate Gaussian.