Flexoelectricity represents the electric polarization response to strain gradient and is governed by a boundary value problem (BVP) of fourth-order partial differential equations (PDEs). Due to the complicities of the BVP, we are frequently prevented from searching for analytical formulas of solutions and have to make a detour to numerical approximations. The well-posedness of the BVP has not been rigorously confirmed yet, although it serves as the cornerstone of numerical methods for flexoelectricity. We address this issue by establishing a theorem indicating a unique solution to the electric Gibbs variational principle (EGVP) that is a weak form of the BVP. The theorem leads us to bear out the solvability and error-estimation formula for some Galerkin methods, e.g., the finite element method (FEM) and isogeometric analysis (IGA). The reliability of the numerical theory is further verified by numerical experiments of benchmark problem, where the convergence rate given by FEM results agrees well with that predicted by error-estimation formulas. This study confirms the well-posedness of the flexoelectric theory and further establishes computational flexoelectricity on rigorous mathematical foundations instead of numerical tests.