We study the solutions of an equation of the form Lu=f, where L is a pseudo-differential operator defined for functions on the unit sphere embedded in a Euclidean space, f is a given function, and u is the desired solution. We give conditions under which the solution exists, and deduce local smoothness properties of u given corresponding local smoothness properties of f, measured by local Besov spaces. We study the global and local approximation properties of the spectral solutions, describe a method to obtain approximate solutions using values of f at points on the sphere and polynomial operators, and describe the global and local rates of approximation provided by our polynomial operators.