This paper reports an iteration-driven method to numerically study the thermal nonlinearity in lithium niobate (LN) based MEMS resonators. In comparison to the state of the art, this technique adopts an approximation-free algorithm and thus more accurately captures the complex nonlinear dynamics that often evades the description by Duffing equation. For the first time, the nonlinearity of LN-based laterally vibrating resonators is theoretically investigated and experimentally validated. The admittance response of both S0 and SH0 mode devices was simulated and measured in this work by forward and backward sweeping the excitation frequency at different power levels. Excellent agreement between simulations and measurements has been achieved.