We present an approximation method for geodesic circles on a spheroid. Our approximation curve is the intersection of two spheroids whose axes are parallel, and it interpolates four points of the geodesic circle. Our approximation method has two merits. One is that the approximation curve can be obtained algebraically, and the other is that the approximation error is very small. For example, our approximation of a circle of radius 1000 km on the Earth has error 1·13 cm or less. We analyze the error of our approximation using the Hausdorff distance and confirm it by a geodesic distance computation.