Functional analysis algorithms as developed by Kantorovich provide a rigorous theory for the convergence of iterative methods, specifically Newton's method, to a wide class of nonlinear functional equations on Banach spaces. The numerical application of the Newton-Kantorovich theorem, quasilinearization, has been demonstrated by Bellman. Here this theorem has also been applied symbolically in order to explore the feasibility of obtaining accurate analytic solutions by automatic formula manipulation. These solution methods are applied to satellite perturbations in polar coordinates or in orbital elements for the relativistic, oblateness, and atmospheric perturbations. The accuracy of these calculations is a priori set by error bounds; e.g., the well-known relativistic pericenter precession and a correction to it are found in one iteration. By virtue of the small deviations of perturbed planetary orbits from a Keplerian orbit, accurate initial approximations can be selected that yield very accurate first iterates, and the inverse Frechet derivatives of the appropriate nonlinear operators are linear and easily constructed. The application to differential equations of satellite geodesy provide a potentially efficient method of obtaining analytical expressions for the dependence of satellite orbits, e.g., on geoid and atmospheric parameters.