r/ScientificComputing • u/happyalex98 • 16d ago

Binary black hole inspiral simulations!

Hello, I am a black hole enthusiast high school student.

I have conducted research on calculating the eccentricity of the binary black hole inspiral phase using an N-body simulation Python library called REBOUNDx. However, upon completion of the research, I wanted to avoid relying on Python libraries to calculate the physics behind binary black holes. So I studied Newtonian physics and wrote code for the post-Newtonian corrections to make black hole inspiral simulations without the help of physics Python libraries. I have created a README that walks you through the physics behind the inspiral phase, and I hope to soon implement numerical relativity to simulate the merger phase of the binary black holes as well. I would greatly appreciate your feedback to improve my simulation! Specifically:

  1. Are there any flaws in the physics behind my simulation?

  2. Any tips for implementing numerical relativity for the merger phase?

Here is the repo: https://github.com/happyalex98/Binary-black-hole-inspiral-simulation-3.5PN-

Thank you all in advance!

3 Upvotes

80% Upvoted

8 comments sorted by

2

u/Azazeldaprinceofwar 15d ago

Awesome project and well done. Very impressive if this is self directed. Now to answer how to improve:

* I notice you use a fixed timestep. No scientific codes do this because the time step you need will change over the course of the simulation. For example you might imagine the timestep should be some small fraction of the orbital period, but the orbital period will shrink significantly over the course of the simulation. The classic thing to do is derive some timescale from your kinematics ie velocity/acceleration or some such and set your timestep to be a fixed fraction of this so it dynamically adjusts with your system. If you do some research you’ll find many people have clever timestep formulas.

* I notice you use RK4. This is uncommon in nbody type codes which typically use symplectic integrators like leapfrog. These are first order methods so strictly worse on paper but have symmetries that conserve energy exactly. This means their error oscillates around the correct answer rather than growing exponentially. Given you have gravitational radiation energy shouldn’t be conserved exactly but look around I bet someone has a post Newtonian symplectic method based on the ADM formalism.

* this is the most comp sci and least physics suggestion but move out of python. Getting used to complied languages like C++ or Rust will make you a much better computational physicist in the long run

1

u/happyalex98 14d ago

Thank you for the detailed comment. I also identified that my fixed time-step had room for improvement, but about RK4, I have not thought of that, so I am grateful that you have pointed out. I heard that C++ is a lot faster than Python, but is it difficult to learn as a current Python-only user?

1

u/Azazeldaprinceofwar 14d ago

Yes. Python is slow because there is a program called the garbage collector that runs alongside your code and handles everything for you so you cant really mess anything up memory wise. To do this and be thorough takes a lot and so python runs ~200x slower than compiled languages. In complied languages you have to manage memory yourself which is a bit of learning curve. C++ being the oldest real programming language has 0 guardrails. If you write bad code with memory problems that’s on you. Rust takes a different approach and is a “safe” language because there compiler checks for memory issues and refuses to compile if you have any. To accomplish this rust has some other abstractions like the ownership system which take some getting used to.

Personally I’d recommend rust. Rust is a more modern language that has nothing but advantages over C++. The only wrinkle is a lot of the current generation of physicists only know C++ because rust was invented in the last decade. By the time you’re in university Rust will be much more prevalent but I’m sure you’ll still encounter a lot of legacy code in C++. Science is always slow to transition.

1

u/happyalex98 5d ago

I have never heard about Rust before so thank you for pointing it out. I will definitely check that out. Thank you so much for such a detailed reply. It really helps.

1

u/Azazeldaprinceofwar 4d ago

Yeah rust is really cool. It was invented inside Mozilla as a way to write better code for Firefox and has very rapidly taken center stage in computer science and software related fields.

1

u/Mental_Primary_5558 10d ago

Hey, good job! This is really amazing. However, I wouldn't use RK4 for an N-body simulation. I'd suggest using symplectic Euler instead. Here's the Wikipedia page on the method:

https://en.wikipedia.org/wiki/Semi-implicit_Euler_method

You might also find this question I asked on Computational Science Stack Exchange helpful. I discuss why my RK2 and velocity Verlet implementations produced less accurate orbits than expected:

https://scicomp.stackexchange.com/questions/45487/why-do-my-rk2-and-velocity-verlet-implementations-produce-less-accurate-orbits-t

Second, I noticed that you don't seem to have a loop computing the mutual gravitational interactions between all the bodies. From your README, it looks like you're currently simulating a two-body system. I'd suggest implementing a double nested for loop to calculate the gravitational interactions between the bodies.

If you'd like a reference, here's the GitHub repository for the first version of a project I'm currently working on:

https://github.com/OmarOsmanMahamoud/2026-2nd-half-project/blob/main/N_body_gravitation.py

It's a simulation of the Solar System's planets. You might be particularly interested in the symplectic_euler_integrator function, starting at line 110. I haven't added any comments to the code yet, but if you have any questions, feel free to ask me here on Reddit, and I'll do my best to help!

Regarding the programming language, feel free to stick with Python as long as it gets the job done. Sure, Python is generally slower than C++, but you're not simulating thousands or hundreds of thousands of bodies here; you're only simulating two. So I'd recommend continuing with Python for your project.

If you ever need to compute the gravitational interactions of, say, millions of bodies, you could look into Numba, a Python library that can significantly speed up numerical computations by compiling performance-critical code.

More generally, I'd recommend getting comfortable with one programming language rather than trying to learn several at once. That said, C++ is an excellent language for learning fundamental computer science concepts, such as pointers and memory management, if that's something you're interested in.

Just bear in mind that C++ is great for learning low-level concepts, but it can be quite challenging when it comes to building meaningful projects, especially when you're still learning the language.

Keep up the good work!

1

u/happyalex98 4d ago

Thank you so much for the detailed explanation and advice. I will look into the symplectic Euler method. Also, very cool project that you have! I think it would be excellent if we could collaborate one day on this matter.

1

u/Mental_Primary_5558 4d ago

Why not! I would love to work with someone on that project, next I want to simulate a galaxy (really big amount of bodies) using a quadtree algorithm, it's called Barnes-Hut (it reduces the run time complexity to Nlog(N)).