To allow everyone to search for content regarding "Celestial Mechanics" and to improve my own expression and calculation skills, BoJone has conceived the theme "Equations and the Universe." This series mainly focuses on using mathematics to discuss astronomical problems in a relatively in-depth manner. In fact, I have always felt that science cannot be perfectly described without formulas (though it shouldn’t be purely formulas either). I remember books like Stephen Hawking’s A Brief History of Time and The Universe in a Nutshell strive to express viewpoints with as few mathematical formulas as possible. This mode is excellent for the general public, but it is difficult for friends who wish to conduct in-depth research. Therefore, I advocate: The universe is calculated!
Every word of this theme is typed by BoJone, including content cited from Introduction to Celestial Mechanics as well as BoJone’s personal insights. Due to the length of the content and limited time, BoJone plans to write and publish it in several installments, making it as accessible as possible so that friends with a basic foundation in calculus can understand it. Here, the first part is released. Due to the rush, there may be some oversights; corrections are welcome!
The two-body problem assumes there are only two celestial bodies, ignoring interference from other bodies, and studies how they move under the influence of universal gravitation (Newtonian classical mechanics). This is the simplest type of problem in celestial mechanics and is currently the only type that can be solved exactly. If the effects of three or more bodies are considered, only approximate analysis or numerical algorithms can be used.
Suppose there is only the Sun S and a certain planet P. Establish an arbitrary spatial coordinate system O-xyz. Let \boldsymbol{r}_s be the position vector of the Sun and \boldsymbol{r}_p be the position vector of the planet. \boldsymbol{r} is the position vector of the planet relative to the Sun. Let M and m represent the masses of the Sun and the planet, respectively. The gravitational force exerted on the Sun by the planet is: \boldsymbol{F}_s = G\frac{Mm}{r^2}\frac{\boldsymbol{r}}{|\boldsymbol{r}|} = G\frac{Mm}{r^3}\boldsymbol{r} Simultaneously, the planet is also subjected to the Sun’s gravity: \boldsymbol{F}_p = G\frac{Mm}{r^2}\frac{-\boldsymbol{r}}{|\boldsymbol{r}|} = -G\frac{Mm}{r^3}\boldsymbol{r} According to Newton’s Second Law, both the Sun and the planet move under the action of these forces: \boldsymbol{F} = \frac{d(m\boldsymbol{v})}{dt} = m\frac{d^2 \boldsymbol{r}}{dt^2} Thus, we have: \begin{aligned} \boldsymbol{F}_s = G\frac{Mm}{r^3}\boldsymbol{r} = M\frac{d^2 \boldsymbol{r}_s}{dt^2} \\ \boldsymbol{F}_p = -G\frac{Mm}{r^3}\boldsymbol{r} = m\frac{d^2 \boldsymbol{r}_p}{dt^2} \end{aligned} Since \boldsymbol{r} = \boldsymbol{r}_p - \boldsymbol{r}_s, we have: \frac{d^2 \boldsymbol{r}}{dt^2} = \frac{d^2}{dt^2}(\boldsymbol{r}_p - \boldsymbol{r}_s) = -\frac{\mu}{r^3}\boldsymbol{r} where \mu = G(M+m). If we establish a Cartesian coordinate system x-y-z with the Sun at the origin, we obtain the system of differential equations for the two-body problem: \begin{aligned} \frac{d^2 x}{dt^2} = -\frac{\mu}{r^3}x \\ \frac{d^2 y}{dt^2} = -\frac{\mu}{r^3}y \\ \frac{d^2 z}{dt^2} = -\frac{\mu}{r^3}z \end{aligned} In these equations, r^2 = x^2 + y^2 + z^2. Now let us look at the solution to this system. From theoretical mechanics, we know that the motion of a planet is under the action of a central force. Motion under a central force is always planar motion. We can, of course, confirm this from the differential equations above. From the last two differential equations, we have: y\ddot z - z\ddot y = 0 \Rightarrow \frac{d}{dt}(y\dot z - z\dot y) = 0 Integrating both sides: y\dot z - z\dot y = A \tag{1} where A is an integration constant. Similarly, we can obtain: z\dot x - x\dot z = B \tag{2} x\dot y - y\dot x = C \tag{3} Equations (1), (2), and (3) are known as the "Angular Momentum Integrals."
Multiplying equations (1), (2), and (3) by x, y, z respectively and adding them together, we get: Ax + By + Cz = 0. This is the equation of a plane passing through the origin (the Sun), which indicates that the planet and the Sun always remain in the same plane! Therefore, we only need to consider the equations in the O-xy plane: \ddot x = -\frac{\mu x}{r^3}, \quad \ddot y = -\frac{\mu y}{r^3} \tag{4} The angular momentum integral reduces to one: x\dot y - y\dot x = h \tag{5}
Now we discuss this using polar coordinates: let x = r\cos\theta, y = r\sin\theta. Transforming (4) yields: \ddot{r} - r\dot{\theta}^2 = -\frac{\mu}{r^2} \tag{6} Transforming (5) yields: r^2 \dot{\theta} = h \tag{7} Equation (7) is called the "Area Integral," where h is twice the areal velocity swept by the radius vector. Since the areal velocity is constant, the area swept in equal time intervals is naturally equal. Thus, we have proven "Kepler’s Second Law."
Note: A brief description of the transformation process:
Let x = r\cos\theta, y = r\sin\theta …(00)
Then x^2 + y^2 = r^2. Substituting into the original equations:
\ddot{x} = -\frac{\mu \cos\theta}{r^2}, \ddot{y} = -\frac{\mu \sin\theta}{r^2} …(0)
Differentiating (00) with respect to t:
\dot{x} = \dot{r}\cos\theta - \dot{\theta}r\sin\theta …(01)
\dot{y} = \dot{r}\sin\theta + \dot{\theta}r\cos\theta …(02)
Differentiating (01) and (02) again with respect to t:
\ddot{x} = \ddot{r}\cos\theta - 2\dot{r}\dot{\theta}\sin\theta - \ddot{\theta}r\sin\theta - \dot{\theta}^2r\cos\theta …(03)
\ddot{y} = \ddot{r}\sin\theta + 2\dot{r}\dot{\theta}\cos\theta + \ddot{\theta}r\cos\theta - \dot{\theta}^2r\sin\theta …(04)
Calculating (03)\cos\theta + (04)\sin\theta gives:
\ddot{x}\cos\theta + \ddot{y}\sin\theta = \ddot{r} - \dot{\theta}^2r …(05)
From (0), we also have:
\ddot{x}\cos\theta + \ddot{y}\sin\theta = -\mu/r^2 …(06)
Comparing (05) and (06), we get: \ddot{r} - \dot{\theta}^2 \cdot r = -\mu/r^2As for transforming equation (5) into (7), it does not require as many techniques; simply substitute the derivatives and simplify.
Now we solve the system of equations formed by (6) and (7). If we wish to find the type of curve for the planetary orbit, we need to find the relationship between r and \theta. If we wish to calculate the planet’s position at a certain time, we must find the relationship between r or \theta and time t.
First, we seek the type of orbital curve. Here I directly excerpt content from Introduction to Celestial Mechanics. Let u = 1/r, then (7) becomes: \dot{\theta} = hu^2, and we have: \begin{aligned} \dot{r} = \frac{dr}{d\theta}\frac{d\theta}{dt} = \frac{d(1/u)}{d\theta}\dot{\theta} = -1/u^2 \cdot \frac{du}{d\theta} \cdot hu^2 = -h\frac{du}{d\theta} \\ \ddot{r} = -h\frac{d}{d\theta}\left(\frac{du}{d\theta}\right)\dot{\theta} = -h^2 u^2\frac{d^2 u}{d\theta^2} \end{aligned} Substituting into (6) gives: \frac{d^2 u}{d\theta^2} + u = \frac{\mu}{h^2} This is a second-order linear differential equation. Its general solution is: u = \frac{\mu}{h^2}[1 + e \cos(\theta - \omega)] \tag{8} where e and \omega are constants to be determined. Reverting to r, it becomes: r = \frac{h^2/\mu}{1 + e \cos(\theta - \omega)} \tag{9}
From analytic geometry, we know that this is a conic section with the origin at one focus (ellipse, hyperbola, or parabola; see Wikipedia: Conic Section). Thus, we have proven "Kepler’s First Law," and in a more general sense (including not only ellipses but also hyperbolas and parabolas). From this, we can derive h^2 = \mu a(1 - e^2), where a is the semi-major axis of the conic section and e is the eccentricity. When \theta = \omega, r is at its minimum, meaning the planet is at perihelion; thus \omega is the argument of perihelion.
Note: How to solve this differential equation?
For details on solving linear differential equations, refer to Wikipedia: Linear Differential Equation.According to the solution method, we can find that the general solution of \frac{d^2 u}{d\theta^2} + u = 0 is y = C_1 \cos\theta + C_2 \sin\theta. Adding the particular solution \frac{\mu}{h^2} from the right side gives the general solution of the differential equation. But how do we convert it to the form of (8)? It is not difficult to find that C_1 \cos\theta + C_2 \sin\theta + \frac{\mu}{h^2} is equivalent to (8), because: e \cos(\theta - \omega) = e\cos\theta \cos\omega + e \sin\theta \sin\omega That is, e\cos\omega = C_1 h^2/\mu and e\sin\omega = C_2 h^2/\mu.
At this point, we have completed most of the work. There is only one small part left to finish solving the two-body problem: finding r or \theta as a function of time t to calculate the planet’s position at any given moment.
Note: The relationship between t and r
From (7) we get \dot{\theta} = h/r^2. Substituting this into (6) gives: \ddot{r} - h^2/r^3 = -\frac{\mu}{r^2} \tag{10} This is a second-order differential equation. Its solution is easy to find, but the integral is quite complex: \dot{r}\frac{d\dot{r}}{dr} = h^2/r^3 - \frac{\mu}{r^2} \dot{r}d\dot{r} = (h^2/r^3 - \frac{\mu}{r^2})dr. Integrating both sides: \dot{r}^2 = 2\mu/r - h^2/r^2 + K_1 \tag{11} \begin{aligned} \Rightarrow dt/dr = \frac{r}{\sqrt{K_1 r^2 + 2\mu r - h^2}} \\ t = \int \frac{r}{\sqrt{K_1 r^2 + 2\mu r - h^2}} dr \end{aligned} This integral is straightforward (refer to the Integral Table on this site), but the disadvantage is that the final result is too complex in form!
Introduction to Celestial Mechanics finds the relationship between \theta and t (i.e., Kepler’s Equation)... (To be continued)
When reposting, please include the original address of this article: https://kexue.fm/archives/549
For more detailed reposting matters, please refer to: Scientific Space FAQ