A New Approximation for the Perimeter of an Ellipse

Abstract

We consider the problem of approximating the perimeter of an ellipse, for which there is no known finite formula, in the context of high-precision performance. Ellipses are broadly used in many fields, like astronomy, manufacturing, medical imaging, and geophysics. They are applied on large and nanoscales, and while numerical integration can be used to obtain precision measurements, having a finite formula can be used for modeling. We propose an iterative symbolic regression approach, utilizing the pioneering work of Ramanujan’s second approximation introduced in 1914 and a known Padé approximation, leading to good results for both low and high eccentricities. Our proposed model is also compared with a very comprehensive historical collection of different approximations collated by Stanislav Sýkora. Compared with the best-known approximations in this centuries-old mathematical problem, our proposed model performs at both extremities while remaining consistent in mid-range eccentricities, whereas existing models excel only at one extremity.

1. Introduction

It can be argued that the idea of establishing “precision machine learning” would not be considered a new endeavor for many scientists, as it has indeed been the aim of mathematicians over many centuries to give the best possible approximations of objects when there is no closed-form formula. For instance, precisely calculating elliptical perimeters without complex integrals has been a long-standing problem in mathematics. The accuracy of these approximations can be critical in astronomy, manufacturing, medical imaging, and geophysics calculations, such as the orbital mechanics of celestial bodies, where small errors can lead to significant deviations in downstream calculations or decisions, affecting everything from scientific research to industrial processes.

Beyond the motivation for greater precision for the theoretical sake of future unknown applications, a note of particular importance is the wide range of scales on which perimeter approximations are necessary. Applications scale from the astronomical order in orbital bodies to atomical in elliptical lens design. Approximations at these scales are also subject to data accuracy, for instance, when visually eliciting ellipse characteristics from images with limited resolution. As hardware improves over time, the viability of applying different potential approaches could also increase, and reassessment of available methods is necessary to maintain the state of the art.

1.1. The Perimeter of the Ellipse

It is well known that while the area, A, of an ellipse with major and minor semi-axes a and b, respectively, is given by $A=\pi ab$; in contrast, there is no such finite formula for the perimeter S when $a\ne b$. An exact expression requires a formula to be written as an infinite series or the evaluation of an integral.

Following Sýkora [1,2], we are aware that many ellipse approximations have been categorized and benchmarked, namely into Keplerian, Padé, extreme, and algebraic categories. Sýkora uses the notation $S(a,b)$ for the perimeter and gives the following formula:

$$S(a,b)=4aE\left(e_c^2\right),$$

where

$$e_c={\displaystyle \frac{\sqrt{a^2-b^2}}{a}}$$

is the eccentricity of the ellipse, and

$$E\left(x\right)={\int }_0^{{\displaystyle \frac{\pi }{2}}}\sqrt{1-x{sin}^2\theta }d\theta$$

is the complete elliptic integral of the second kind.

For more than 400 years, researchers such as Kepler, Euler, and Peano, among many others, have proposed approximations [3], a subject that was even covered in letters in the journal Nature over a century ago [4,5,6].

1.2. Ramanujan’s Approximations

In 1914 [7], Srinivasa Ramanujan first proposed two approximations of interest, given in Equations (1) and (2):

$$S(a,b)\approx \pi \left(3(a+b)-\sqrt{(3a+b)(a+3b)}\right),$$

(1)

$$S(a,b)\approx m_R=\pi (a+b)\left(1+{\displaystyle \frac{3h}{10+\sqrt{4-3h}}}\right),$$

(2)

where

$$h={\left({\displaystyle \frac{a-b}{a+b}}\right)}^2.$$

(3)

This second approximation, $m_R(a,b)$, certainly stood the test of time, and although it requires the computation of a square root, it is admired for its simplicity and precision. Ramanujan presented it as

$$S(a,b)\approx \pi \left((a+b)+{\displaystyle \frac{3{(a-b)}^2}{10(a+b)+\sqrt{a^2+14ab+b^2}}}+ϵ\right),$$

(4)

where $ϵ$ is about $3a{e}_c^{20}/68719476736$. It is well known that Ramanujan’s approximation underestimates the perimeter, and a detailed error analysis is now known thanks to M.B. Villarino [8].

1.3. Symbolic Regression and Learning from Data

Symbolic regression (SR) is a particular type of mathematical modeling from data that is aimed at discovering mathematical expressions that best fit a given dataset. The emphasis is also on trying to reduce the complexity of the models. Instead of assuming a specific model form (like linear or polynomial), symbolic regression searches through the space of mathematical expressions to find one that explains the data. This is often carried out using evolutionary algorithms, such as genetic programming, which iteratively build upon equations by combining mathematical operators, constants, and variables.

Symbolic regression demonstrates strong performance in scenarios with limited data, contrasting with artificial neural networks (ANNs), which typically require additional strategies for few-shot or few-sample learning (FSL) to function effectively in such conditions [9,10]. Unlike the opaque models ANNs generate, SR produces explicit mathematical equations that domain experts can readily analyze, facilitating detailed inspection of variable relationships within the data.

A key advantage of SR lies in its capacity to inform further investigation and guide subsequent data collection, as the relationship between the variables can be analyzed in the formulae. This is particularly relevant in physical phenomena, where low-order, simple functions often govern underlying models. SR’s ability to train on minimal samples, derive concise, interpretable formulae, and maintain robustness in the presence of measurement noise positions it as an effective method for identifying underlying models. These characteristics can lead to fundamental insights and discoveries within the data.

In this work, we configure TuringBot software [11], which utilizes a novel simulated annealing metaheuristic rather than an evolutionary approach to uncover models. Our use of this software is detailed further in later sections.

1.4. Structure of the Paper

The process for utilizing symbolic regression to uncover perimeter approximations is described in Section 2, including the iterative method using upper and lower bounds. The series of approximations uncovered are described in Section 3, where it is also benchmarked against Sýkora’s collection of ellipse approximations. We consider the properties of the residuals in high precision from a circle to ellipses with very large eccentricity. In the final section, Section 5, we make our final remarks and suggest potential downstream research.

2. Materials and Methods

We used the 2.16.1 version of TuringBot in our experimentation. We note its commercial nature and evolving updates, which suggest that the application will be more robust than an unsupported and unevolving academic software package for SR. Additionally, the multithreading software leverages multiple CPU cores, enabling parallel evaluations of candidate models. This increases efficiency and allows the software to explore a broader section of the solution space over time by exploring more potential solutions simultaneously.

TuringBot displays candidate formulas and provides errors according to the selected objective function. The software evolves and refines these formulas over multiple iterations. Formulas are ranked based on a Pareto frontier, providing multiple solutions. The user can select those of better accuracy or simplicity, balancing the resulting models’ interpretability and effectiveness. Due to the software’s commercial and proprietary nature, we are unable to provide exact details of the novel simulated annealing approach that produces and develops the candidate formulas [11].

The dataset defined in

Table 1. The table shows the absolute relative error for each approximation of an ellipse for different eccentricities within the training samples, where a varies but $b=1$ is constant. The smallest errors are in bold within each row.

a	${\mathit{m}}_{\mathit{R}}$ Equation (2) Err	${\mathit{m}}_{\mathit{u}}$ Equation (5) Err	${\mathit{m}}_{\mathit{Ru}}$ Equation (8) Err	${\mathit{m}}_{\mathit{p}}$ Equation (9) Err	${\mathit{m}}_{\mathit{MC}}$ Equation (12) Err
1.05	1.707$\times {10}^{-21}$	2.793$\times {10}^{-07}$	1.707$\times {10}^{-21}$	$\mathbf{1}.\mathbf{743} \times {\mathbf{10}}^{-\mathbf{28}}$	$\mathbf{1}.\mathbf{743} \times {\mathbf{10}}^{-\mathbf{28}}$
1.15	6.296$\times {10}^{-17}$	2.284$\times {10}^{-06}$	6.296$\times {10}^{-17}$	$\mathbf{4}.\mathbf{317} \times {\mathbf{10}}^{-\mathbf{22}}$	$\mathbf{4}.\mathbf{317} \times {\mathbf{10}}^{-\mathbf{22}}$
1.25	6.679$\times {10}^{-15}$	5.784$\times {10}^{-06}$	6.679$\times {10}^{-15}$	$\mathbf{2}.\mathbf{959} \times {\mathbf{10}}^{-\mathbf{19}}$	$\mathbf{2}.\mathbf{959} \times {\mathbf{10}}^{-\mathbf{19}}$
1.35	1.268$\times {10}^{-13}$	1.037$\times {10}^{-05}$	1.268$\times {10}^{-13}$	$\mathbf{1}.\mathbf{825} \times {\mathbf{10}}^{-\mathbf{17}}$	$\mathbf{1}.\mathbf{825} \times {\mathbf{10}}^{-\mathbf{17}}$
1.05	1.707$\times {10}^{-21}$	2.793$\times {10}^{-07}$	1.707$\times {10}^{-21}$	$\mathbf{1}.\mathbf{743} \times {\mathbf{10}}^{-\mathbf{28}}$	$\mathbf{1}.\mathbf{743} \times {\mathbf{10}}^{-\mathbf{28}}$
1.15	6.296$\times {10}^{-17}$	2.284$\times {10}^{-06}$	6.296$\times {10}^{-17}$	$\mathbf{4}.\mathbf{317} \times {\mathbf{10}}^{-\mathbf{22}}$	$\mathbf{4}.\mathbf{317} \times {\mathbf{10}}^{-\mathbf{22}}$
1.25	6.679$\times {10}^{-15}$	5.784$\times {10}^{-06}$	6.679$\times {10}^{-15}$	$\mathbf{2}.\mathbf{959} \times {\mathbf{10}}^{-\mathbf{19}}$	$\mathbf{2}.\mathbf{959} \times {\mathbf{10}}^{-\mathbf{19}}$
1.35	1.268$\times {10}^{-13}$	1.037$\times {10}^{-05}$	1.268$\times {10}^{-13}$	$\mathbf{1}.\mathbf{825} \times {\mathbf{10}}^{-\mathbf{17}}$	$\mathbf{1}.\mathbf{825} \times {\mathbf{10}}^{-\mathbf{17}}$
1.45	1.049$\times {10}^{-12}$	1.574$\times {10}^{-05}$	1.049$\times {10}^{-12}$	3.516$\times {10}^{-16}$	$\mathbf{3}.\mathbf{516} \times {\mathbf{10}}^{-\mathbf{16}}$
1.55	5.328$\times {10}^{-12}$	2.166$\times {10}^{-05}$	5.328$\times {10}^{-12}$	3.422$\times {10}^{-15}$	$\mathbf{3}.\mathbf{422} \times {\mathbf{10}}^{-\mathbf{15}}$
1.65	1.966$\times {10}^{-11}$	2.794$\times {10}^{-05}$	1.966$\times {10}^{-11}$	2.130$\times {10}^{-14}$	$\mathbf{2}.\mathbf{130} \times {\mathbf{10}}^{-\mathbf{14}}$
1.75	5.799$\times {10}^{-11}$	3.445$\times {10}^{-05}$	5.796$\times {10}^{-11}$	9.684$\times {10}^{-14}$	$\mathbf{9}.\mathbf{683} \times {\mathbf{10}}^{-\mathbf{14}}$
1.85	1.450$\times {10}^{-10}$	4.109$\times {10}^{-05}$	1.448$\times {10}^{-10}$	3.493$\times {10}^{-13}$	$\mathbf{3}.\mathbf{492} \times {\mathbf{10}}^{-\mathbf{13}}$
1.95	3.193$\times {10}^{-10}$	4.778$\times {10}^{-05}$	3.186$\times {10}^{-10}$	1.056$\times {10}^{-12}$	$\mathbf{1}.\mathbf{055} \times {\mathbf{10}}^{-\mathbf{12}}$
2	4.560$\times {10}^{-10}$	5.112$\times {10}^{-05}$	4.546$\times {10}^{-10}$	1.738$\times {10}^{-12}$	$\mathbf{1}.\mathbf{737} \times {\mathbf{10}}^{-\mathbf{12}}$
3	3.298$\times {10}^{-08}$	1.122$\times {10}^{-04}$	3.107$\times {10}^{-08}$	6.960$\times {10}^{-10}$	$\mathbf{6}.\mathbf{687} \times {\mathbf{10}}^{-\mathbf{10}}$
4	2.499$\times {10}^{-07}$	1.584$\times {10}^{-04}$	2.009$\times {10}^{-07}$	1.180$\times {10}^{-08}$	$\mathbf{9}.\mathbf{571} \times {\mathbf{10}}^{-\mathbf{09}}$
5	8.514$\times {10}^{-07}$	1.923$\times {10}^{-04}$	5.339$\times {10}^{-07}$	6.511$\times {10}^{-08}$	$\mathbf{3}.\mathbf{926} \times {\mathbf{10}}^{-\mathbf{08}}$
6	1.965$\times {10}^{-06}$	2.175$\times {10}^{-04}$	8.884$\times {10}^{-07}$	2.081$\times {10}^{-07}$	$\mathbf{8}.\mathbf{627} \times {\mathbf{10}}^{-\mathbf{08}}$
7	3.629$\times {10}^{-06}$	2.363$\times {10}^{-04}$	1.080$\times {10}^{-06}$	4.867$\times {10}^{-07}$	$\mathbf{1}.\mathbf{325} \times {\mathbf{10}}^{-\mathbf{07}}$
8	5.821$\times {10}^{-06}$	2.506$\times {10}^{-04}$	9.776$\times {10}^{-07}$	9.345$\times {10}^{-07}$	$\mathbf{1}.\mathbf{600} \times {\mathbf{10}}^{-\mathbf{07}}$
9	8.487$\times {10}^{-06}$	2.614$\times {10}^{-04}$	5.271$\times {10}^{-07}$	1.570$\times {10}^{-06}$	$\mathbf{1}.\mathbf{572} \times {\mathbf{10}}^{-\mathbf{07}}$
10	1.156$\times {10}^{-05}$	2.695$\times {10}^{-04}$	2.684$\times {10}^{-07}$	2.398$\times {10}^{-06}$	$\mathbf{1}.\mathbf{181} \times {\mathbf{10}}^{-\mathbf{07}}$
11	1.497$\times {10}^{-05}$	2.756$\times {10}^{-04}$	1.371$\times {10}^{-06}$	3.415$\times {10}^{-06}$	$\mathbf{4}.\mathbf{146} \times {\mathbf{10}}^{-\mathbf{08}}$
12	1.865$\times {10}^{-05}$	2.800$\times {10}^{-04}$	2.728$\times {10}^{-06}$	4.610$\times {10}^{-06}$	$\mathbf{7}.\mathbf{110} \times {\mathbf{10}}^{-\mathbf{08}}$
13	2.254$\times {10}^{-05}$	2.831$\times {10}^{-04}$	4.281$\times {10}^{-06}$	5.968$\times {10}^{-06}$	$\mathbf{2}.\mathbf{153} \times {\mathbf{10}}^{-\mathbf{07}}$
14	2.660$\times {10}^{-05}$	2.852$\times {10}^{-04}$	5.977$\times {10}^{-06}$	7.471$\times {10}^{-06}$	$\mathbf{3}.\mathbf{853} \times {\mathbf{10}}^{-\mathbf{07}}$
15	3.078$\times {10}^{-05}$	2.864$\times {10}^{-04}$	7.768$\times {10}^{-06}$	9.104$\times {10}^{-06}$	$\mathbf{5}.\mathbf{744} \times {\mathbf{10}}^{-\mathbf{07}}$
16	3.504$\times {10}^{-05}$	2.869$\times {10}^{-04}$	9.615$\times {10}^{-06}$	1.085$\times {10}^{-05}$	$\mathbf{7}.\mathbf{758} \times {\mathbf{10}}^{-\mathbf{07}}$
17	3.935$\times {10}^{-05}$	2.869$\times {10}^{-04}$	1.148$\times {10}^{-05}$	1.269$\times {10}^{-05}$	$\mathbf{9}.\mathbf{830} \times {\mathbf{10}}^{-\mathbf{07}}$
18	4.370$\times {10}^{-05}$	2.864$\times {10}^{-04}$	1.335$\times {10}^{-05}$	1.461$\times {10}^{-05}$	$\mathbf{1}.\mathbf{190} \times {\mathbf{10}}^{-\mathbf{06}}$
19	4.805$\times {10}^{-05}$	2.855$\times {10}^{-04}$	1.519$\times {10}^{-05}$	1.660$\times {10}^{-05}$	$\mathbf{1}.\mathbf{393} \times {\mathbf{10}}^{-\mathbf{06}}$
20	5.239$\times {10}^{-05}$	2.843$\times {10}^{-04}$	1.699$\times {10}^{-05}$	1.864$\times {10}^{-05}$	$\mathbf{1}.\mathbf{586} \times {\mathbf{10}}^{-\mathbf{06}}$
21	5.671$\times {10}^{-05}$	2.829$\times {10}^{-04}$	1.874$\times {10}^{-05}$	2.072$\times {10}^{-05}$	$\mathbf{1}.\mathbf{767} \times {\mathbf{10}}^{-\mathbf{06}}$
22	6.100$\times {10}^{-05}$	2.812$\times {10}^{-04}$	2.044$\times {10}^{-05}$	2.284$\times {10}^{-05}$	$\mathbf{1}.\mathbf{934} \times {\mathbf{10}}^{-\mathbf{06}}$
23	6.524$\times {10}^{-05}$	2.794$\times {10}^{-04}$	2.206$\times {10}^{-05}$	2.499$\times {10}^{-05}$	$\mathbf{2}.\mathbf{085} \times {\mathbf{10}}^{-\mathbf{06}}$
24	6.943$\times {10}^{-05}$	2.775$\times {10}^{-04}$	2.362$\times {10}^{-05}$	2.715$\times {10}^{-05}$	$\mathbf{2}.\mathbf{220} \times {\mathbf{10}}^{-\mathbf{06}}$
25	7.356$\times {10}^{-05}$	2.754$\times {10}^{-04}$	2.511$\times {10}^{-05}$	2.932$\times {10}^{-05}$	$\mathbf{2}.\mathbf{337} \times {\mathbf{10}}^{-\mathbf{06}}$
26	7.763$\times {10}^{-05}$	2.733$\times {10}^{-04}$	2.652$\times {10}^{-05}$	3.150$\times {10}^{-05}$	$\mathbf{2}.\mathbf{436} \times {\mathbf{10}}^{-\mathbf{06}}$
27	8.164$\times {10}^{-05}$	2.710$\times {10}^{-04}$	2.787$\times {10}^{-05}$	3.367$\times {10}^{-05}$	$\mathbf{2}.\mathbf{518} \times {\mathbf{10}}^{-\mathbf{06}}$
28	8.558$\times {10}^{-05}$	2.688$\times {10}^{-04}$	2.914$\times {10}^{-05}$	3.585$\times {10}^{-05}$	$\mathbf{2}.\mathbf{584} \times {\mathbf{10}}^{-\mathbf{06}}$
29	8.945$\times {10}^{-05}$	2.665$\times {10}^{-04}$	3.035$\times {10}^{-05}$	3.801$\times {10}^{-05}$	$\mathbf{2}.\mathbf{634} \times {\mathbf{10}}^{-\mathbf{06}}$
30	9.325$\times {10}^{-05}$	2.641$\times {10}^{-04}$	3.149$\times {10}^{-05}$	4.016$\times {10}^{-05}$	$\mathbf{2}.\mathbf{669} \times {\mathbf{10}}^{-\mathbf{06}}$
40	1.274$\times {10}^{-04}$	2.407$\times {10}^{-04}$	3.976$\times {10}^{-05}$	6.060$\times {10}^{-05}$	$\mathbf{2}.\mathbf{398} \times {\mathbf{10}}^{-\mathbf{06}}$
50	1.554$\times {10}^{-04}$	2.194$\times {10}^{-04}$	4.395$\times {10}^{-05}$	7.852$\times {10}^{-05}$	$\mathbf{1}.\mathbf{560} \times {\mathbf{10}}^{-\mathbf{06}}$
60	1.783$\times {10}^{-04}$	2.009$\times {10}^{-04}$	4.578$\times {10}^{-05}$	9.393$\times {10}^{-05}$	$\mathbf{6}.\mathbf{352} \times {\mathbf{10}}^{-\mathbf{07}}$
70	1.974$\times {10}^{-04}$	1.851$\times {10}^{-04}$	4.626$\times {10}^{-05}$	1.072$\times {10}^{-04}$	$\mathbf{1}.\mathbf{905} \times {\mathbf{10}}^{-\mathbf{07}}$
80	2.134$\times {10}^{-04}$	1.714$\times {10}^{-04}$	4.595$\times {10}^{-05}$	1.186$\times {10}^{-04}$	$\mathbf{8}.\mathbf{676} \times {\mathbf{10}}^{-\mathbf{07}}$
90	2.271$\times {10}^{-04}$	1.596$\times {10}^{-04}$	4.521$\times {10}^{-05}$	1.285$\times {10}^{-04}$	$\mathbf{1}.\mathbf{399} \times {\mathbf{10}}^{-\mathbf{06}}$
100	2.390$\times {10}^{-04}$	1.493$\times {10}^{-04}$	4.421$\times {10}^{-05}$	1.372$\times {10}^{-04}$	$\mathbf{1}.\mathbf{805} \times {\mathbf{10}}^{-\mathbf{06}}$
500	3.572$\times {10}^{-04}$	4.225$\times {10}^{-05}$	1.780$\times {10}^{-05}$	2.295$\times {10}^{-04}$	$\mathbf{1}.\mathbf{409} \times {\mathbf{10}}^{-\mathbf{06}}$
1000	3.784$\times {10}^{-04}$	2.249$\times {10}^{-05}$	1.004$\times {10}^{-05}$	2.470$\times {10}^{-04}$	$\mathbf{6}.\mathbf{122} \times {\mathbf{10}}^{-\mathbf{07}}$
10,000	3.998$\times {10}^{-04}$	2.422$\times {10}^{-06}$	1.158$\times {10}^{-06}$	2.649$\times {10}^{-04}$	$\mathbf{1}.\mathbf{644} \times {\mathbf{10}}^{-\mathbf{08}}$
100,000	4.023$\times {10}^{-04}$	2.453$\times {10}^{-08}$	1.577$\times {10}^{-08}$	2.671$\times {10}^{-04}$	$\mathbf{5}.\mathbf{054} \times {\mathbf{10}}^{-\mathbf{15}}$
1,000,000,000	4.023$\times {10}^{-04}$	2.454$\times {10}^{-11}$	3.957$\times {10}^{-09}$	2.671$\times {10}^{-04}$	$\mathbf{3}.\mathbf{974} \times {\mathbf{10}}^{-\mathbf{15}}$

was used as input to the software to produce the symbolic models. The mean square error was the objective function used to minimize, and all samples were used during exploration. Only the basic operators were permitted in the exploration: multiplication, division, addition, subtraction, and exponentiation. Only integer coefficients were considered in the produced models.

Our results are based on a proposed method in which an iterative scheme of approximations is produced. We used symbolic regression and/or a known upper or lower bound at each step. Although the models were produced using the limited samples from Table 1, we evaluated them on uniformly sampled data in log space over the range $a\in [1,1\times {10}^9]$ and presented the results in plots with logarithmic scales on both axes.

The historical models curated by Sýkora, as well as the models produced by TuringBot, were subsequently re-implemented and re-evaluated in Mathematica at 200 points of decimal precision. When implementing the models from TuringBot at higher precision in Mathematica, the performance of some models declined and were excluded from consideration. We believe precision limitations within TuringBot led to this behavior, and it was something that affected only a few models that were generated in the last iteration of our approach.

Our Mathematica implementation validated the models from TuringBot and completed the approximations in Sýkora’s MATLAB implementation [1] by including S-class and algebraic equations in an independent implementation in an alternate language. We have made it available in the code repository accompanying this submission at https://github.com/impvau/ellipse (accessed 12 September 2024).

All Mathematica calculations were performed with a precision of 200 decimal places, including the “true” ellipse perimeter values computed with the EllipticE function. Using this EllipticE function allowed us to compute the exact perimeter values to the extent that we could compute the required integrals, which were subject to the numerical methods employed in this process. We confirmed that the computed values in Mathematica and MATLAB at 200 decimal places of accuracy were identical for all digits.

3. Results

This project started out of curiosity when one of the authors (PM) watched a YouTube video describing Equation (2). Arriving at Equation (8) using only the data for different values of a in Table 1 was relatively straightforward, so it is reasonable to give a more detailed rationale of the process involved since its simplicity has a certain beauty. Also, it may be potentially generalizable for other problems. Consequently, we first discuss the results of an approximation that can be obtained from Ramanujan’s second approximation. This defines the first iteration of our method and the first set of comparative results.

3.1. A New Asymptotic Bound Leads to a More Accurate Approximation for High Eccentricity

Ramanujan has not left any record on the rational process by which he proposed Equations (1) and (2) apart from mentioning that he used an “empirical” method (although his mastery in the use of continued fractions [12] may be behind it).

We propose here a new approximation; let $m_u$ be an upper bound defined as follows:

$$m_u=\pi (a+b)\left(1+{\displaystyle \frac{(44/\pi -11)h}{10+\sqrt{4-3h}}}\right),$$

(5)

The $m_u$ model was produced using the samples in Table 1. We note that most a values are below 100, while the last six points range between 100 and $1\times {10}^9$ to facilitate modeling the asymptotic behavior. Table 1 also shows approximation errors for key formulae mentioned in this section. Figure 1 shows the errors between $a=1$ and $a=80$, while Figure 2 shows the errors where both axes are in a log scale and the functions are re-evaluated using 10,000 uniform points in log space between $a=1$ and $a=1\times {10}^9$ with $b=1$.

Next, we take an iterative approach where we define the first iteration $h_0$ as

$$h_0={\left({\displaystyle \frac{2\pi a-2\pi b}{2\pi a+2\pi b}}\right)}^2=h,$$

(6)

arising from the crude lower and upper bounds $2\pi b$ and $2\pi a$, respectively.

Relying on the knowledge that has accumulated since 1914 to find an even better lower bound that can be used instead of Ramanujan’s second approximation, we now define $h_1$ as follows:

$$h_1={\left({\displaystyle \frac{m_R-m_u}{m_R+m_u}}\right)}^2.$$

(7)

Once again, using TuringBot software, we find the following approximation:

$$S(a,b)\approx m_{Ru}=m_R+\alpha m_u{h}^{\beta }{h}_1^{\gamma },\alpha =9938,\beta =7,\gamma =1$$

(8)

3.2. Use of a Highly Accurate Padé Approximation Formula for Low Eccentricities

Another approximation that is more accurate than Ramanujan’s second approximation and exceptionally good for low eccentricity values is also known (and it is included in Table 1 for comparison). In line with previous proposals by Bronshtein (1964) and Jacobsen (1985) (covered in [3]), a Padé approximation-based equation for the ‘normalized’ perimeter (i.e., the perimeter divided by $\pi (a+b)$) as a function of h is

$$S\approx m_p=\pi (a+b)\left({\displaystyle \frac{135168-85760h-5568h^2+3867h^3}{135168-119552h+22208h^2-345h^3}}\right).$$

(9)

Let us now analogously define $h_2$ as follows:

$$h_2={\left({\displaystyle \frac{m_p-m_u}{m_p+m_u}}\right)}^2.$$

(10)

Once again, using symbolic regression and the TuringBot package, we have obtained a new approximation given by the following formula:

$$S\approx m_i=m_p{(m_u/m_p)}^{{\left({\displaystyle \frac{h_1^2}{615}}\right)}^{\left({\displaystyle \frac{h_1/h_2}{a}}\right)}}.$$

(11)

Utilizing $m_i$, we obtain the following improved approximation, valid for any ellipse (except for a circle):

$$S\approx m_{MC}=m_i\left(1+{\displaystyle \frac{h^{18}}{a^2-a^3}}\right).$$

(12)

Once again, the integer coefficients for the formula for $m_{MC}$ and its functional form were found via symbolic regression using the same software cited before with the values of a given in Table 1 and $b=1$, as well as the corresponding exact solutions obtained via multiple precision calculations involving the integral definition of the perimeter. As a consequence, the values obtained for $m_{MC}$ are also displayed in Table 1. It is clear now that $m_{MC}$ maintains good asymptotic behavior for low a values, while the relative error for an extremely high eccentricity of $a={10}^9$ is of the order of ${10}^{-13}$.

In this case, we searched for formulas that, in addition to the basic arithmetic functions of multiplication, division, addition, and subtraction, can also raise to a power. We used as variables the set $\{h,h_1,h_2,m_u,m_R,m_p\}$, while we continued to only use integers as coefficients in the formulas.

The range of a values can be thought to be “measured in units of b” and covers a wide range of eccentricities in orbits around the Sun, including the comet Ikeya–Seki (C/1965 S1-A), as per https://theskylive.com/ikeyaseki-info (accessed on 12 September 2024), which has an eccentricity of 0.999915 and an a in units of b of approximately 76.7.

We note for the periodic comet called 103P/Hartley 2 (https://www.wikiwand.com/en/articles/103P/Hartley (accessed on 12 September 2024)), with an eccentricity of 0.694 and an a in units of b of approximately 1.39, $m_{MC}$ has a relative error of the order of ${10}^{-16}$, thus maintaining the good approximation behavior of Equation (9), but for highly eccentric ellipses, clearly the introduction of our simple $m_u$ (Equation (5)) approximation formula gave $m_{MC}$ a better asymptotic behavior. This motivated us to look more widely to see if other proposed approximations also do well for very eccentric ellipses.

4. Comparative Results with Sýkora’s Curated Collection of Approximations

Sýkora has put together an extensive collection of ellipse perimeter approximations and provides a classification into groups that share properties. Some approximations may fall into multiple categories (for instance, both Keplerian and Padé) and are explicitly classified as such in those circumstances.

We define in this section a few selected examples in each category of Sýkora’s collection, as well as all approximations from the S-class. A comprehensive list is provided in the code repository accompanying this submission, and later in Figure 3, Figure 4 and Figure 5, the nature of the residuals for each class can be seen. Note that Sýkora’s analysis considers the ratio of the semi-axes $y=b/a$.

4.1. Keplerian Equations

Keplerian approximations have been defined as approximations “whose relative error curves attain absolute maximum at $y=0$ (degenerate ellipses) and decay towards zero at $y=1$ (circles)”. Examples are given as follows:

$$k_1=2\pi \sqrt{ab},$$

$$k_2=2\pi {\displaystyle \frac{{(a+b)}^2}{{(\sqrt{a}+\sqrt{b})}^2}},$$

$$k_3=\pi (a+b),$$

$$\dots$$

4.2. Keplerian Padè Equations

Padé approximations are classified “by fractions $n/m$, where n and m are the degrees of the polynomials in h appearing, respectively, in the nominator and denominator of the rational function”. Some examples are

$$k_a=\pi (a+b){\displaystyle \frac{16+3h}{16-h}},$$

$$k_b=\pi (a+b){\displaystyle \frac{64+16h}{64-h^2}},$$

$$k_c=\pi (a+b){\displaystyle \frac{64-3h^2}{64-16h}},$$

$$\dots$$

4.3. Exact Extremes (No Crossing) Equations

“Exact extreme” approximations are exact for both $y=1$ and $y=0$. Sýkora divides them into those with “crossing” or “no crossing”, where the latter means “$S^{\prime }(a,b)$ remains …$S(a,b)$” ($S^{\prime }$ being the approximation to S). A few examples are given as follows:

$$e_1={\displaystyle \frac{\pi (a-b)}{arctan\left(\frac{a-b}{a+b}\right)}},$$

$$e_2=4(a+b)-{\displaystyle \frac{(8-2\pi )ab}{0.410117(a+b)+(1-2\times 0.410117)\left(\frac{\sqrt{(a+74b)(74a+b)}}{1+74}\right)}},$$

$$e_3=2\sqrt{{\pi }^{2}ab+4{(a-b)}^2},$$

$$\dots$$

4.4. Combined Padè Equations with Exact Extremes (No Crossing)

Examples of this class follow:

$$e_{ad1}={\displaystyle \frac{\pi }{4}}\left({\displaystyle \frac{81}{64}}\right)-1,e_{ad2}={\displaystyle \frac{\pi }{4}}\left({\displaystyle \frac{19}{15}}\right)-1,e_{ap}={\displaystyle \frac{e_{ad1}}{e_{ad1}-e_{ad2}}},$$

$$e_a=\pi (a+b)\left(e_{ap}{\displaystyle \frac{16+3h}{16-h}}+(1-e_{ap}){\left(1+{\displaystyle \frac{h}{8}}\right)}^2\right),$$

$$e_{bd1}={\displaystyle \frac{\pi }{4}}\left({\displaystyle \frac{80}{63}}\right)-1,e_{bd2}={\displaystyle \frac{\pi }{4}}\left({\displaystyle \frac{61}{48}}\right)-1,e_{bp}={\displaystyle \frac{e_{bd1}}{e_{bd1}-e_{bd2}}},$$

$$e_b=\pi (a+b)\left(e_{bp}{\displaystyle \frac{64-3h^2}{64-16h}}+(1-e_{bp}){\displaystyle \frac{64+16h}{64-h^2}}\right),$$

$$\dots$$

4.5. Exact Extremes and Crossing Equations

These are exact models where $S^{\prime }(a,b)$ does not always remain either above or below $S(a,b)$.

$$c_1=\pi \sqrt{2\left(a^2+b^2\right)}\left({\displaystyle \frac{sin\left(c_{1t}\right)}{c_{1t}}}\right),c_{1t}={\displaystyle \frac{\pi }{4}}\left({\displaystyle \frac{a-b}{b}}\right),$$

$$c_2=4a+2(\pi -2)a{\left({\displaystyle \frac{b}{a}}\right)}^{1.456},$$

$$c_3=4\left({\displaystyle \frac{\pi ab+{(a-b)}^2}{a+b}}\right)-{\displaystyle \frac{89}{146}}{\left({\displaystyle \frac{b\sqrt{a}-a\sqrt{b}}{a+b}}\right)}^2,$$

$$\dots$$

4.6. Algebraic Equations

The following Sýkora “algebraic” approximations are those “which use only the four basic operations of the real numbers algebra”.

$$a_1={\displaystyle \frac{4(a^2+b^2+(\pi -2)ab)}{a+b}}$$

$$a_2={\displaystyle \frac{4(a^3+b^3+a_{2t1}ab(a+b))}{a^2+b^2+2a_{2s1}ab}},a_{2t1}=2.49808365277126,a_{2s1}=1.22694921875000,$$

$$a_{3t1}=6.16881239339582,a_{3t2}=4.06617730084445,a_{3s1}=6.15241658169936,$$

$$a_3={\displaystyle \frac{4(a^4+b^4+a_{3t1}ab(a^2+b^2)+2a_{3t2}{a}^2{b}^2)}{a^3+b^3+a_{3s1}ab(a+b)}},$$

$$\dots$$

4.7. All S-Class Equations

The S-class approximations contain approximations identified after the prior classes’ classification and are limited to “a single square root evaluation”. These are the most performant approximations from Sýkora’s approximations regarding the maximum absolute relative error.

$$s_{0q}=4.16102118885517,s_{0u}={\displaystyle \frac{4-\pi }{4-2\sqrt{2+s_{0q}}}},s_{0p}=1-s_{0u},$$

$$s_0=s_{0p}(a+b)+s_{0u}\sqrt{a^2+b^2+s_{0q}ab},$$

$$s_{1q}=92.28480788617108,s_{1v}=0.04522028227769,s_{1p2}=0.99983439391729,$$

$$s_{1u}=1+s_{1v}-s_{1p2},s_{1p1}={\displaystyle \frac{\pi }{2}}(2+s_{1v}\sqrt{2+s_{1q}})-(2s_{1p2}+2s_{1u}\sqrt{2+s_{1q}}),$$

$$s_1={\displaystyle \frac{s_{1p2}(a^2+b^2)+s_{1p1}ab+s_{1u}(a+b)\sqrt{a^2+b^2+s_{1q}ab}}{(a+b)+s_{1v}\sqrt{a^2+b^2+s_{1q}ab}}}$$

$$s_{2q1}=13.66022044408346,s_{2q2}=37.30921886231118,s_{2v}=-1.03788930003090,$$

$$s_{2t}=5.51954143485218,s_{2p2}=1.24957869093182,s_{2u}=1+s_{2v}-s_{2p2},$$

$$s_{2p1}={\displaystyle \frac{\pi }{4}}(2+s_{2t}+s_{2v}\sqrt{2+2s_{2q1}+s_{2q2}})-(s_{2p2}+s_{2u}\sqrt{2+2s_{2q1}+s_{2q2}}),$$

$$s_{2R}=\sqrt{(a^4+b^4)+s_{2q1}ab(a^2+b^2)+s_{2q2}{a}^2{b}^2},$$

$$s_2={\displaystyle \frac{s_{2p2}(a^3+b^3)+s_{2p1}ab(a+b)+s_{2u}(a+b)s_{2R}}{(a^2+b^2)+s_{2t}ab+s_{2v}{s}_{2R}}}$$

$$s_{1rq}=74.01745125408363,s_{1rv}=0.05027328304233,s_{1ru}=s_{1rv},$$

$$s_{1rp1}={\displaystyle \frac{\pi }{2}}(2+s_{1rv}\sqrt{2+s_{1rq}})-(2+2s_{1ru}\sqrt{2+s_{1rq}}),$$

$$s_{1r}={\displaystyle \frac{(a^2+b^2)+s_{1rp1}ab+s_{1ru}(a+b)\sqrt{a^2+b^2+s_{1rq}ab}}{(a+b)+s_{1rv}\sqrt{a^2+b^2+s_{1rq}ab}}}$$

$$s_{ac1}=\pi -3,s_{ac2}=\pi ,s_{ac3}=0.5,s_{ac4}={\displaystyle \frac{1+\pi }{2}},s_{ac5}=4,$$

$$s_{ak}=1-{\displaystyle \frac{s_{ac1}ab}{(a^2+b^2)+s_{ac2}\sqrt{s_{ac3}{\left(ab\right)}^2+ab\sqrt{ab\left(s_{ac4}(a^2+b^2)+s_{ac5}ab\right)}}}},$$

$$s_a=4{\displaystyle \frac{\pi ab+s_{ak}{(a-b)}^2}{a+b}}$$

$$s_{aoc1}=0.14220038049945,s_{aoc2}=3.30596250119242,s_{aoc3}=0.00135657637724,$$

$$s_{aoc4}=2.00637978782056,s_{aoc5}=5.3933761426286,$$

$$s_{aok}=1-{\displaystyle \frac{s_{aoc1}ab}{(a^2+b^2)+s_{aoc2}\sqrt{s_{aoc3}{\left(ab\right)}^2+ab\sqrt{ab\left(s_{aoc4}(a^2+b^2)+s_{aoc5}ab\right)}}}},$$

$$s_{ao}=4{\displaystyle \frac{\pi ab+s_{aok}{(a-b)}^2}{a+b}}$$

$$s_{ard1}=0.14220096377128,s_{ard2}=3.93490847789660,s_{ard3}=2.691437204515743,$$

$$s_{ark}=1-{\displaystyle \frac{s_{ard1}ab}{(a^2+b^2)+s_{ard2}\sqrt{\sqrt{{\left(ab\right)}^3(a^2+b^2+s_{ard3}ab)}}}},$$

$$s_{ar}=4{\displaystyle \frac{\pi ab+s_{ark}{(a-b)}^2}{a+b}}$$

4.8. A Detailed Comparison of Results with the Best Performers in the S-Class of Sýkora’s Collection

Towards the end of our benchmarking study, we looked at other alternative proposals described by Sýkora. We have continued to group the existing formulas and compare them within the same range. Sýkora compares approaches by measuring the maximum relative absolute error in parts per unit (generally million or billion). These functions and categories are defined in the code repository accompanying this submission.

We have taken our best model produced in the prior sections and show the comparison to each category in Figure 3, Figure 4 and Figure 5, as well as Sýkora’s best performers (the S-class) in Figure 6.

Sýkora’s analysis focuses on the maximum point of errors; however, we focus on consistency at the extremities and in the regions of critical practical performance ($a\le 80$) for astronomical importance. The categories in Figure 3, Figure 4 and Figure 5 show that the approach excels at approximating the extremities of a. We observe that we generally out-compete all categories, except for $e_4$ out-competing at very high a values, and notably, $c_6$ is the only method offering a close approximation at both extremities. In most circumstances, $m_{MC}$ out-competes across the entire spectrum of a, except for the region before $a=1000$, where many methods outperform $m_{MC}$.

Interestingly, in the most competitive S-class comparison in Figure 6, we exclude $S_0$, $S_1$, and $S_2$, which fail to achieve good performance when evaluated at high precision (see code repository accompanying this submission). The unoptimized Ahmedi $s_a$ is perhaps the closest to our approximation $m_{MC}$ due to its performance at high a values, although it is not as competitive at low a values.

5. Discussion

Le Verrier’s 1846 work predicting Neptune due to perturbations in Uranus’s orbit is a famous example of the importance of mathematical modeling in scientific discovery [13]. Key orbital parameters, such as the semi-major axis and eccentricity, were central to his calculations, demonstrating how additional gravitational forces affect planetary paths. Considering that later work in [14] suggests that comparable accuracy can be achieved with less numerical calculation and given that we can easily perform computations directly from integral representations or series expansions with modern computers, we might ask ourselves then, “Why is this quest for high-precision approximations still relevant”?

In our case, we are motivated by the development of new methods for “precision machine learning” [15] since high-precision approximations are often crucial to scientific discovery. We are also motivated by the quest of finding approximations with simple models since sometimes we have functions that are defined via integral representations or series expansions; having a simple model can help if these computations need to be carried out for several values of parameters many times during simulations.

Given these motivations, the problem of approximating the ellipse’s perimeter seems to be an ideal test bed of new ideas and a challenging one with centuries of approximations given, and we hope our work encourages others to keep moving in the same direction.

Overall, this case study indicates that an iterative approach to creating high-performing approximations using symbolic regression can (at least sometimes) be achieved, given existing data, by recursively computing the upper and lower bounds of the function of interest interspersed with the computation of extra variables as in Equations (3) and (7) to be incrementally added to generate new and more precise bounding models. This procedure can be understood as a constructive heuristic approach to building a “network” of increasingly better approximations, but, in conjunction with many other connectionist approaches, it requires symbolic regression at different steps to build intermediate models, each with an increasing approximation.

In addition, we also found that minimization of the mean relative error may be useful to produce approximations that work well for a large range of parameters. The functional forms of Equations (8) and (9) suggests that continued fraction regression [16] using dynamic depth strategies [17] can be used to obtain more concise and precise approximations iteratively.

As a side note, the search process did not result in a simpler and/or equally fitting approximation than Equation (8), using symbolic regression, even when we used other ellipse characteristics such as the eccentricity or the semi-latus rectum. We note that nonlinear optimization may also be helpful to identify other values of $\alpha$, $\beta$, and $\gamma$ that can even improve the approximation for particular ranges of a measured in units of b that have specific applications. For example, in fitting planetary orbiting objects, a training set that uses the most typical values of a could possibly lead to other coefficient values that would better approximate for the particular range of interest.