In this work, we explore fluctuations during phase transitions of uniaxial and biaxial liquid crystals using a phenomenological free energy functional. We rely on a continuum-level description of the liquid crystal ordering with a tensorial parameter and a temperature dependent Landau polynomial expansion of the tensor's invariants. The free energy functional, over a three-dimensional periodic domain, is integrated with a Gaussian quadrature and minimized with a theoretically informed Monte Carlo method. We reconstruct analytical phase diagrams, following Landau and Doi's notations, to verify that the free energy relaxation reaches the global minimum. Importantly, our relaxation method is able to follow the thermodynamic behavior provided by other non-phenomenological approaches; we predict the first order character of the isotropic-nematic transition, and we identify the uniaxial-biaxial transition as second order. Finally, we use a finite-size scaling method, using the nematic susceptibility, to calculate the transition temperatures for 4-Cyano-4'-pentylbiphenyl (5CB) and N-(4-methoxybenzylidene)-4-butylaniline (MBBA). Our results show good agreement with experimental values, thereby validating our minimization method. Our approach is an alternative towards the relaxation of temperature dependent continuum-level free energy functionals, in any geometry, and can incorporate complicated elastic and surface energy densities.