Abstract: The International GNSS Service (IGS) precise ephemeris provides satellite coordinates every 15 minutes. In order to improve the positioning accuracy, it is often necessary to obtain the satellite position at any time. Interpolation and fitting of IGS precise ephemeris is a common method to obtain satellite coordinates of continuous epochs. The improved Legendre polynomial algorithm is used to fit the satellite orbit coordinates, and compared with the conventional algorithm. The results show that the conventional algorithm can maintain high accuracy only when the fitting order is low. When the fitting period is 6 h, LU decomposition method and singular value decomposition (SVD) method can maintain high accuracy in solving singular matrix, while when the fitting period is 12 h, SVD decomposition method decomposes the matrix \boldsymbolB with low condition number to obtain polynomial coefficient \boldsymbolC, so as to avoid the error caused by ill conditioned matrix, so it can still maintain high accuracy. In high-order fitting, SVD decomposition method is superior to LU decomposition method and conventional algorithm in both accuracy and stability.