Building a solar system simulator

[Xilor]

TL;DR Summary

I'm trying to build a solar system simulator and have questions

Hiya
Hiya. I'm trying to build a solar system simulator for research purposes, and this does not seem to be all that easy. As in, my firsts tests have the moon shooting out of the solar system... I was wondering if you folks could offer some guidance.

So, since it's for research, accuracy is the most important aspect. Velocities of the objects in the solar system are especially important to me, and I'd like to be able to get as close as possible to the true values over some timespans. Ideally, I'd like to be able to reach a precision of the movements of masses within about 0.1 m/s of the true value over a simulation period of a year. I imagine this optimistic precision is orders of magnitude away from what I can actually manage, but it's the figure I'd like to keep in mind while building the simulator.
Because of the things I'd like to test, I basically need to use gravitational acceleration, I can't use Kepler's planetary motion or similar hacks.So, problem one. Using timesteps obviously introduces small errors that build on each other. How do I best prevent these inaccuracies? Should I just use the Euler method with small timesteps, or is there some better method I could use (in terms of precision / execution time)?

Problem two. What masses do I actually need for my desired accuracy level? Would I need the satellites of other planets? How about the Kuiper belt? (and how would I model that?) Dwarf planets? How about stars in 'close' proximity to the sun or the galaxy itself? For performance reasons, and for ease of development, I'd obviously like as few masses as necessary. Would it be ok to just add the masses of the moons of other planets to their mass?

Problem three. I'm going to need a good grond truth for positions and velocities for both the start and end time of my simulation for the masses I decide to use. Both to get everything started realistically and to check how accurate I am. Ideally, the ground truth is something I can access programatically so I can easily alter my simulation period. The positions/velocities will also have to be in x,y,z coordinates or need to be convertable in those. I don't neccessarily care if the 'ground truth' actually perfectly matches the true historical data. As long as the 'end truth' is a correct (within my desired accuracy) result of the 'start truth' given what we know, then it works for me.
The best library I found so far is this https://github.com/mgvez/planet-positions which is used in the jsOrrery solar system simulation, and I'm using it now, but I'm not sure if it's usable/reliable enough. For starters, it only has data for the Sun, planets and moon, so I couldn't add larger moons like Titan if I went with this. Does anyone know of a better way to obtain ground truth data?

Problem four. Is Newtonian gravitational acceleration actually enough for the desired precision? Are there any GR effects I need to consider for my desired accuracy? How about radiation pressure? SR velocity addition? Anything else I might've missed?

 

[Mister T]

Summary:: I'm trying to build a solar system simulator and have questions

Are you treating this as N-body problem?

So, problem one. Using timesteps obviously introduces small errors that build on each other. How do I best prevent these inaccuracies? Should I just use the Euler method with small timesteps, or is there some better method I could use (in terms of precision / execution time)?

Small time steps that are powers of two. So, for example, instead of a step size of 0.001 use $2^{-10}$. You could also use something like the Runge-Kutta method instead of the Euler method.

 

[DrStupid]

So, problem one. Using timesteps obviously introduces small errors that build on each other. How do I best prevent these inaccuracies? Should I just use the Euler method with small timesteps, or is there some better method I could use (in terms of precision / execution time)?

Euler is not suitable for your problem. I would suggest higher order Runge-Kutta or Runge-Kutta-Nyström methods with adaptive stepsize.

 

[Dale]

Are you insistent on coding the method yourself or would a library be acceptable?

 

[Xilor]

Are you treating this as N-body problem?

Yep

Small time steps that are powers of two. So, for example, instead of a step size of 0.001 use $2^{-10}$.

What would the advantage of a power of two step size be? Is that for variable precision? I was intending to use C#'s decimal type if it matters

You could also use something like the Runge-Kutta method instead of the Euler method.

Euler is not suitable for your problem. I would suggest higher order Runge-Kutta or Runge-Kutta-Nyström methods with adaptive stepsize.

Thanks, I'll try and see if I can implement something like that. Having looked at them briefly, I'm wondering whether the correct move would be to do these kind of half-steps and adjustments only on the single mass I'm dealing with, or whether I should be calculating the new positions of all masses at these half steps to account for their movement too. Do these Runge-Kutta methods account for changing environments?

Are you insistent on coding the method yourself or would a library be acceptable?

There's some things I'd like to be able to tweak, so I'd rather use my own code unless the library is highly adaptable. Libraries could be great sources of inspiration though, so if you know of any particularly good ones I'd love to hear of them

 

[DrStupid]

Having looked at them briefly, I'm wondering whether the correct move would be to do these kind of half-steps and adjustments only on the single mass I'm dealing with, or whether I should be calculating the new positions of all masses at these half steps to account for their movement too.

You need to performing each stage of the method for all bodies at once.

Do these Runge-Kutta methods account for changing environments?

Can you specify "changing environments"?

 

[Xilor]

You need to performing each stage of the method for all bodies at once.

Okay great, thanks!

Can you specify "changing environments"?

It was just an attempted clarification of the rest of the question. Like how the positions of the other masses influence things. It isn't just the position + velocity of a single mass that matters. If you'd treat the masses one at a time, the other masses would be 'environment'
But I see now that the positions and velocities of all the masses are basically a single high dimensional vector together so it makes sense that you'd want to it all at once.

 

[Dale]

Libraries could be great sources of inspiration though, so if you know of any particularly good ones I'd love to hear of them

For Python I would start with the jplephem package. At a minimum that would give you a “ground truth” to compare to. The JPL ephemeris is solid.

 

[Xilor]

For Python I would start with the jplephem package. At a minimum that would give you a “ground truth” to compare to. The JPL ephemeris is solid.

Ah! This data is awesome, thanks a ton! Using data from this package for the initial conditions even stopped my moon from flying off, to great rejoicing of my virtual moonians. I'll be checking out the different data packages to see which one fits my use case most

 

[Mister T]

What would the advantage of a power of two step size be?

The computer does its arithmetic in base 2, so you reduce rounding error. I'm not sure if that's relevant unless you're writing your own code.

If you really are interested in doing this as an N-body problem then that answers most of your questions. Every one of the objects you mentioned is a separate body. Otherwise you will have to make some approximations, such as how planets, moons, and other objects are to be entered into your equations.

 

[Filip Larsen]

Just a small comment that if you are aiming for high accuracy in your solver then you may want to consider using a symplectic integrator to minimize energy drift in your solutions i.e. the situation where truncation errors in the integrator "conspire" to make an accumulating non-physical change in system energy and/or other conserved quantities. Unfortunately, explicite Runge-Kutta methods are not symplectic so if you stay with those methods you need to find another way to correct for energy drift.

I am not aware of the current state-of-art for symplectic N-body solvers, but perhaps others here can point you to a good library, package or similar.

 

[Xilor]

The computer does its arithmetic in base 2, so you reduce rounding error. I'm not sure if that's relevant unless you're writing your own code.

If you really are interested in doing this as an N-body problem then that answers most of your questions. Every one of the objects you mentioned is a separate body. Otherwise you will have to make some approximations, such as how planets, moons, and other objects are to be entered into your equations.

I guess it makes sense to base 2, thanks for the tip!
Not really sure how thinking of it as an N-body problem makes any of the questions clearer, I already was from the start. It does also potentially make the asteroid/kuiper belts and the remainder of the galaxy possibly awkward. But I doubt I'll add these anyway.

Just a small comment that if you are aiming for high accuracy in your solver then you may want to consider using a symplectic integrator to minimize energy drift in your solutions i.e. the situation where truncation errors in the integrator "conspire" to make an accumulating non-physical change in system energy and/or other conserved quantities. Unfortunately, explicite Runge-Kutta methods are not symplectic so if you stay with those methods you need to find another way to correct for energy drift.

I am not aware of the current state-of-art for symplectic N-body solvers, but perhaps others here can point you to a good library, package or similar.

That does sound pretty important, even on the relatively short timescales I'm working with. I'll try and look into this when I'm reaching my accuracy limits.

So results so far: With the help of y'all, I managed to get a decently functional model working. Even using timesteps as high as 200 seconds over the course of a year, most larger objects are basically where they should be and at approximately the right velocities. If I don't include their moons, then Saturn, Jupiter, Uranus and Neptune actually all fall within my desired 0.1m/s accuracy!
The moon and mercury are the most problematic out of the main objects I'd like to deal with, which isn't entirely unexpected I guess. Mercury is off by as much as 3km/s on the 200 second timestep. It goes down to 105 m/s on a 8 second time, but that's still way off my goal. I'd like to try even lower timesteps, but it's going to take a while.

Strangely, Runge-Kutta has actually been decreasing my accuracy performance, while also taking several times as long. I may have made an error somewhere, but I've been hunting it for a while with no success.

Since they should have some minor impact on the orbits of the inner planets, I'd kind of like to add some of the most important asteroid masses, especially Ceres, Pallas and Vesta. But I can't seem to find a dataset that works with the jplephem library that contains them. I found one that's supposed to contain them (codes_300ast_20100725.bsp) but it's giving me errors.

Right now I'm mainly left wondering whether to include radiation pressure and whatever causes the perihelion shift of Mercury in GR. Does anyone know how to calculate the latter?

 

[Dale]

If I don't include their moons, then Saturn, Jupiter, Uranus and Neptune actually all fall within my desired 0.1m/s accuracy!

Excellent work!

Right now I'm mainly left wondering whether to include radiation pressure and whatever causes the perihelion shift of Mercury in GR. Does anyone know how to calculate the latter?

I would strongly doubt that this is the issue. To check if that is the issue then calculate your calculated precession of perihelion. If your Newtonian calculation is correct it should be 532 arcsec per Julian century. The actual observed precession is 574 arcsec per century. The GR correction is 43 arcsec per century. I suspect that your current calculation is nowhere near 532, so adding GR will not fix it.

 

[Xilor]

I would strongly doubt that this is the issue. To check if that is the issue then calculate your calculated precession of perihelion. If your Newtonian calculation is correct it should be 532 arcsec per Julian century. The actual observed precession is 574 arcsec per century. The GR correction is 43 arcsec per century. I suspect that your current calculation is nowhere near 532, so adding GR will not fix it.

Hmm, I'm definitely curious now. Guess I'll go and attempt to calculate that for a single orbit of Mercury and extrapolate (running a full Julian century on tiny timesteps would take ages). I'm guessing both the location of the perihelion and its precession are measured from the sun, and not from the solar systems center?

Meanwhile I tried to look for the GR adjustment I'd need and found this paper: https://arxiv.org/pdf/1803.01678.pdf which seems to suggest (formulas 14,15,16) that the adjusted formula is as simple as:
acceleration = (GM/r^2) * (1 + 3*(x X v)^2/(r^2c^2))
Which I'll probably end up trying anyway just for curiosities sake, even if it doesn't alter results much

 

[Xilor]

I would strongly doubt that this is the issue. To check if that is the issue then calculate your calculated precession of perihelion. If your Newtonian calculation is correct it should be 532 arcsec per Julian century. The actual observed precession is 574 arcsec per century. The GR correction is 43 arcsec per century. I suspect that your current calculation is nowhere near 532, so adding GR will not fix it.

So I just kept it running for a while on a 16 secs timestep, and it does seem to get closeish. On the first 70 orbits of Mercury the average seemed to slowly converge to the 500-550 arcsec area (there was quite a bit of variance per orbit). After that (until orbit 135 which is where I stopped) it started moving away towards 400ish arcseconds, I'm guessing because positions of some objects got messed up enough that it's not really accurate anymore.
With even shorter timesteps, which I'll use for the actual runs, I imagine it'll get quite a bit closer to the real value. So it does seem like a 43 arcsec per century difference could be noticeable. It's probably not the main contributor to my error right now, but if it results in noticeable accuracy gains I'll probably run with it

 

[Dale]

Wow, I am impressed! To get that level of accuracy on a home-brewed implementation this quickly is excellent.