Euler's Solutions

In the 1760s, Leonhard Euler made the first serious attempt to solve the three-body problem after Newton's formulation. Euler introduced a crucial simplification — the restricted three-body problem, in which the mass of one body is negligible compared to the other two. This approximation proved extremely fruitful for celestial mechanics.

Restricted Three-Body Problem

Euler assumed that one of the three bodies has negligible mass (m₃ → 0) and does not affect the motion of the two massive bodies. The two massive bodies orbit their common center of mass following Kepler's laws, while the small body moves in their gravitational field. This simplification reduces the problem to the motion of a single body in a time-varying gravitational field.

Euler's Method

Euler also developed the simplest numerical integration method for ordinary differential equations (ODEs), now known as Euler's method. For an ODE system dy/dt = f(t, y), the method approximates the solution by moving forward along the tangent line given by the derivative f(t_n, y_n) at each point, using a fixed step size h.

yn+1=yn+h⋅f(tn,yn)\mathbf{y}_{n+1} = \mathbf{y}_n + h \cdot \mathbf{f}(t_n, \mathbf{y}_n)

Significance for the Three-Body Problem

Euler's work on the restricted three-body problem became the foundation for lunar theory and asteroid orbit calculations. His numerical integration method, though rarely used today in its pure form, remains the conceptual basis for more advanced methods. This simulation uses Euler's method as a tribute to his pioneering work.

Numerical Example

First Euler step (Δt = 0.016). At t = 0: body 1 at (10, 0), v = (0, 14.14), m = 100; body 2 at (−10, 0), v = (0, −14.14), m = 100; body 3 at (0, 120), v = (−11.49, 0), m = 0.1. G = 80 (for simulation convenience; real G = 6.674 × 10⁻¹¹).

Gravitational force between massive bodies 1 and 2:

F12=Gm1m2r122=80⋅100⋅100202=2000F_{12} = G \frac{m_1 m_2}{r_{12}^2} = 80 \cdot \frac{100 \cdot 100}{20^2} = 2000

Body 1: acceleration from F₁₂, Euler step:

a1=F12m1=2000100=20a1=(−20,  0)\begin{aligned} a_1 &= \frac{F_{12}}{m_1} = \frac{2000}{100} = 20 \\ \mathbf{a}_1 &= (-20,\; 0) \end{aligned}
v1(Δt)=v1(0)+a1Δt=(0,  14.14)+(−20,  0)⋅0.016=(−0.32,  14.14)\begin{aligned} \mathbf{v}_1(\Delta t) &= \mathbf{v}_1(0) + \mathbf{a}_1 \Delta t \\ &= (0,\; 14.14) + (-20,\; 0) \cdot 0.016 \\ &= (-0.32,\; 14.14) \end{aligned}
r1(Δt)=r1(0)+v1(0)Δt=(10,  0)+(0,  14.14)⋅0.016=(10,  0.226)\begin{aligned} \mathbf{r}_1(\Delta t) &= \mathbf{r}_1(0) + \mathbf{v}_1(0) \Delta t \\ &= (10,\; 0) + (0,\; 14.14) \cdot 0.016 \\ &= (10,\; 0.226) \end{aligned}

Body 2: acceleration from F₁₂ (symmetric, opposite direction):

a2=F12m2=2000100=20a2=(20,  0)\begin{aligned} a_2 &= \frac{F_{12}}{m_2} = \frac{2000}{100} = 20 \\ \mathbf{a}_2 &= (20,\; 0) \end{aligned}
v2(Δt)=v2(0)+a2Δt=(0,  −14.14)+(20,  0)⋅0.016=(0.32,  −14.14)\begin{aligned} \mathbf{v}_2(\Delta t) &= \mathbf{v}_2(0) + \mathbf{a}_2 \Delta t \\ &= (0,\; -14.14) + (20,\; 0) \cdot 0.016 \\ &= (0.32,\; -14.14) \end{aligned}
r2(Δt)=r2(0)+v2(0)Δt=(−10,  0)+(0,  −14.14)⋅0.016=(−10,  −0.226)\begin{aligned} \mathbf{r}_2(\Delta t) &= \mathbf{r}_2(0) + \mathbf{v}_2(0) \Delta t \\ &= (-10,\; 0) + (0,\; -14.14) \cdot 0.016 \\ &= (-10,\; -0.226) \end{aligned}

Body 3: net acceleration from bodies 1 and 2 (x-components cancel):

∣a31∣=∣a32∣=Gm1r2=80⋅100102+1202=0.552a3x=0.552⋅1014500−0.552⋅1014500=0a3y=−0.552⋅12014500−0.552⋅12014500=−1.100a3=(0,  −1.100)\begin{aligned} |\mathbf{a}_{31}| &= |\mathbf{a}_{32}| = \frac{G m_1}{r^2} = \frac{80 \cdot 100}{10^2 + 120^2} = 0.552 \\ a_{3x} &= 0.552 \cdot \frac{10}{\sqrt{14500}} - 0.552 \cdot \frac{10}{\sqrt{14500}} = 0 \\ a_{3y} &= -0.552 \cdot \frac{120}{\sqrt{14500}} - 0.552 \cdot \frac{120}{\sqrt{14500}} = -1.100 \\ \mathbf{a}_3 &= (0,\; -1.100) \end{aligned}
v3(Δt)=v3(0)+a3Δt=(−11.49,  0)+(0,  −1.100)⋅0.016=(−11.49,  −0.0176)\begin{aligned} \mathbf{v}_3(\Delta t) &= \mathbf{v}_3(0) + \mathbf{a}_3 \Delta t \\ &= (-11.49,\; 0) + (0,\; -1.100) \cdot 0.016 \\ &= (-11.49,\; -0.0176) \end{aligned}
r3(Δt)=r3(0)+v3(0)Δt=(0,  120)+(−11.49,  0)⋅0.016=(−0.184,  120)\begin{aligned} \mathbf{r}_3(\Delta t) &= \mathbf{r}_3(0) + \mathbf{v}_3(0) \Delta t \\ &= (0,\; 120) + (-11.49,\; 0) \cdot 0.016 \\ &= (-0.184,\; 120) \end{aligned}

Simulation