diff --git a/content/resources/mind-bytes/suncrash/featured.jpg b/content/resources/mind-bytes/suncrash/featured.jpg new file mode 100644 index 00000000000..87e62e0b4f2 Binary files /dev/null and b/content/resources/mind-bytes/suncrash/featured.jpg differ diff --git a/content/resources/mind-bytes/suncrash/index.md b/content/resources/mind-bytes/suncrash/index.md new file mode 100644 index 00000000000..fdb470dd006 --- /dev/null +++ b/content/resources/mind-bytes/suncrash/index.md @@ -0,0 +1,273 @@ +--- +# Documentation: https://sourcethemes.com/academic/docs/managing-content/ + +title: "The day the Earth stood still" +subtitle: "The advantages an ease of just-in-time compilation in Python" +summary: "" +authors: ["reto"] +tags: ["python", "physics"] +categories: ["Physics", "Coding"] +date: "2023-09-06" +lastmod: "2023-09-06" +featured: false +draft: false + +# Featured image +# To use, add an image named `featured.jpg/png` to your page's folder. +# Focal points: Smart, Center, TopLeft, Top, TopRight, Left, Right, BottomLeft, Bottom, BottomRight. +image: + caption: "" + focal_point: "" + preview_only: false + +# Projects (optional). +# Associate this post with one or more of your projects. +# Simply enter your project's folder or file name without extension. +# E.g. `projects = ["internal-project"]` references `content/project/deep-learning/index.md`. +# Otherwise, set `projects = []`. +projects: [] +--- + +Assume the following thought experiment: +The Earth circling the Sun at a cozy distance of 1 AU (around 150 million km) +suddenly stands still and doesn't move anymore. +The only force left on the Earth is now the Sun's gravitational pull. +The question we want to answer here is: How long will it take for the Earth to crash into the Sun? + + + +## An analytical solution + +The Sun with mass \$m\_s\$ and the Earth with mass \$m\_e\$ are located at distance \$r\$ from each other. +In this scenario (an leaving all vector arrows away since everything takes place in one dimension here), +the gravitational force between the two bodies can be written as + +\begin{equation} + F = G \frac{m_s m_e}{r^{2}}. +\end{equation} + +The force and acceleration the Earth is feeling is + +\begin{equation} + F = m_e a +\end{equation} + +and setting both equations equal we get + +\begin{equation} + a = \ddot{r} = G \frac{m_s}{r^2}. +\end{equation} + +Note that the mass of the Earth dropped out of the equation. +We now have a differential equation that we could solve in order to determine the free fall time. +The solution when solving for the time as a function of position can be seen, +e.g., on [Wikipedia](https://en.wikipedia.org/wiki/Free_fall#Inverse-square_law_gravitational_field). + + +## Kepler's third law + +A common way to solve this problem in astronomy is +by using +[Kepler's third law](https://en.wikipedia.org/wiki/Kepler%27s_laws_of_planetary_motion) +of planetary motions, +which states that the square of a planet's orbital period +is proportional to the cube of the length of the semi-major axis of the planet's orbit. +Let us assume that the Earth's orbit is a circle +(which is close to reality). +The semi-major axis is therefore simply twice the previously defined distance \$r\$. +Furthermore, the orbital period of the Earth is 1 year (or 365.24 days). + +The free fall time \$t\$ can be calculated by imagining that the Earth travels around the Sun +in a degenerate elliptic orbit, +such that it travels radially inward towards it and then on the same line back out again. +Using Kepler's third law allows us to write the following equation: + +\begin{equation} + \frac{(2t)^2}{r^3} = \frac{p_e^2}{(2r)^3} +\end{equation} + +Here, \$p\_e\$ is the orbital period of the Earth, +i.e., 1 year. +Solving this equation for \$t\$ results in + +\begin{equation} + t = \sqrt{\frac{p_e^2}{32}} = 64.52\\,\\mathrm{d}. +\end{equation} + + +## A numerical solution + +If we do not know about Kepler's third law and do not want to solve the differential equation, +we can still solve this problem numerically. +The difficulty we need to keep in mind is that the acceleration of Earth is not constant +but depends on the distance \$r\$. +We can split the distance between the Sun and the Earth into small segments \$dr\$. +Within such a segment, we can assume that the acceleration is constant. +This allows us to apply the equation for distance traveled under constant acceleration: + +\begin{equation} + dr = \frac{a}{2} dt^2 + v_0 dt +\end{equation} + +Here, \$v\_0\$ is the initial velocity of the Earth at the beginning of the segment under consideration +and is zero for the first step (see the title of this post). +The time it takes to traverse segment \$dr\$ is \$dt\$. +Rewriting above equation and solving for the time \$dt\$ results in + +\begin{align} + \frac{a}{2} dt^2 + v_0 dt - dr &= 0 \\\\ + dt_{1,2} &= \frac{-v_0 \pm \sqrt{v_0^2 + 2 a dr}}{a}. +\end{align} + +Clearly, only the positive solution is of interest here +since it is the only one that results in a positive time \$t\$. + +To calculate the time it takes for the Earth to crash into the Sun +we can now simply sum over all segments \$dr_i\$, +calculate the time \$dt_i\$ it takes to cross segment \$i\$, +and then calculate the total time \$t\$ as: + +\begin{equation} + t = \sum dt_i. +\end{equation} + +### An implementation in Python + +The following Python code implements this algorithm: + +```python +import time + +import numpy as np + +# Constants +GRAV = 6.67408e-11 # m^3 kg^-1 s^-2 +M_SUN = 1.98847e30 # kg +R_START = 1.495979e11 # m + +# Step size +dr = 1e3 # m + +def total_time(dr): + # Initial conditions + t_total = 0 # s: total time + v_curr = 0 # m/s: initial velocity + r_curr = R_START # m: initial distance + + # Calculate time until Earth crashes into the Sun + while r_curr > 0: + # Calculate acceleration + a_curr = GRAV * M_SUN / r_curr**2 + + # Calculate time it takes to travel dr + t_curr = (-v_curr + np.sqrt(v_curr**2 + 2 * a_curr * dr)) / a_curr + + # Update total time + t_total += t_curr + + # Update velocity + v_curr += a_curr * t_curr + + # Update distance + r_curr -= dr + + return t_total + +# Print result + +tic = time.time() +t_total = total_time(dr) +toc = time.time() +print(f"Time until Earth crashes into the Sun: {t_total / 3600 / 24:.2f} days") +print(f"Computation time: {toc - tic:.2f} seconds") +``` + +Here we wrote the actual calculation into a little function, +which makes it easier to time it. +Running this routine on my computer results in the following output: + +``` +Time until Earth crashes into the Sun: 64.57 days +Computation time: 136.57 seconds +``` + +While we get the correct answer, +the computation time is over two minutes. +Let's see if we can speed this up a bit. + + +### A faster implementation in Python using [`numba`](https://numba.pydata.org/) + +The +[`numba`](https://numba.pydata.org/) +package allows us to take the function `total_time(dr)` +and compile it to machine code before running it. +Here, just-in-time compilation is used, +which means that the compilation happens at runtime when the function is used for the first time. + +The only thing we need to change is to import the package +and add a decorator to the function: + +```python +... + +from numba import njit + +... + +@njit +def total_time(dr): + ... +``` + +Three dots in above example indicate that the rest of the code is unchanged. +Running the code again results in the following output: + +``` +Time until Earth crashes into the Sun: 64.57 days +Computation time: 1.84 seconds +``` + +The computation time is now only 1.84 seconds, +which is a speedup of a factor of 74! + +We can ask the question how long it actually took to compile the function. +To find out, we can simply run the calculation twice and only measure the time for the second run: +The output from this is: + +``` +Time until Earth crashes into the Sun: 64.57 days +Computation time: 1.57 seconds +``` + +Just-in-time compilation thus took around 250 ms. + +### A comparison with Rust + +Just as a comparison, here is the output when implementing the same algorithm in Rust +and building the code with the `--release` flag: + + +``` +Time until Earth crashes into the Sun: 64.56 days +Computation time: 1.58s. +``` + +The computation time is basically identical to the +[`numba`](https://numba.pydata.org/) +implementation in Python. + +This clearly shows the power of +[`numba`](https://numba.pydata.org/) +and just-in-time compilation. + +### Source code + +The final python and Rust source code can be found +[here on GitHub](https://gist.github.com/trappitsch/35d04ef43a8a53a2bd9ea513032da61a). + + +### Resources + +- [Website for `numba`](https://numba.pydata.org/) +- [Documentation for `numba`](https://numba.readthedocs.io/en/stable/user/index.html) \ No newline at end of file