How can the orbital motion of planets around the sun be calculated taking into account the force exerted by other planets?

Question

In The Feynman Lectures on Physics (Vol. 1, Chapter 9), Feynman explained how to calculate the planetary motion by just using the Newton's laws of motion and the Euler's method as a numerical solution to the differential equations. He first explained how to calculate the simple scenario of one planet orbiting another massive body (such as Earth around the sun), supposing that the massive body is infinitely heavy and is therefore not moving. This is easy since the force that the planet feels is always directed towards the fixed location massive body, and the only body in motion is the planet. As a matter of fact, I was able to model the motion of the Earth around the sun to a good approximation.

Then, Feynman explained how to calculate the motion of any planet taking into account the force exerted by more than one body in space:

Now let us see how we can calculate the motion of Neptune, Jupiter, Uranus, or any other planet. If we have a great many planets, and let the sun move too, can we do the same thing? Of course we can. We calculate the force on a particular planet, let us say planet number i, which has a position xi,yi,zi (i=1 may represent the sun, i=2 Mercury, i=3 Venus, and so on). We must know the positions of all the planets. The force acting on one is due to all the other bodies which are located, let us say, at positions xj,yj,zj. Therefore the equations are

Where r(ij) is the distance between the body 'i' in question and the body 'j'. He then said the following:

Also, ∑ means a sum over all values of j—all other bodies—except, of course, for j=i. Thus all we have to do is to make more columns, lots more columns. We need nine columns for the motions of Jupiter, nine for the motions of Saturn, and so on. Then when we have all initial positions and velocities we can calculate all the accelerations from Eq. (9.18) by first calculating all the distances, using Eq. (9.19).

The thing is that I am confused about how would this be done. This is not as simple as the case of the Earth around the sun since in here, every other body is moving! Thus, consider for example that I wanted to calculate the motion of Uranus taking into account the forces exerted by Neptune, Jupiter, and the sun. For simplicity we could still suppose that the sun does not move. But, Uranus, Neptune, and Jupiter would certainly be moving. Therefore, the data required to model this motion would be the initial position and initial velocities of Uranus, Neptune, Jupiter, and the sun (which if we suppose that it doesn't move, the velocity would be 0). Then, what would be the next step? This is what confuses me.

If we calculate the initial acceleration (at time t) that Uranus feels due to all the other forces, we can calculate its velocity and position at one time interval (∆t) later (at t+∆t). But at this time t+∆t, now Jupiter and Neptune have both changed its position, so to calculate the position and velocity of Uranus another time interval later (at t+2∆t), the new position of Jupiter and Neptune (at t+∆t) would have to be taken into account, but how can their positions be calculated?

What I think is that I should apply the 2nd Law of Newton for each planet, and in that way, it would be possible to apply Euler's method to solve the differential equation of each of the bodies in question, but I am confused about how this should be done. I will try to describe my strategy to solve this problem: Taking into account the positions of all the bodies involved, the initial acceleration for each planet 'i' (a_i(t)) can be calculated. Therefore, using Euler's method or another numerical method for solving differential equations, the velocities and positions at a time t+∆t can be calculated for each planet 'i' (x_i(t+∆t) and v_i(t+∆t)). Now, with all these new positions, the acceleration at the time t+∆t can be calculated for each planet 'i' (a_i(t+∆t)), and the cycle can be repeated 'n' times to calculate the position, velocity, and acceleration of each planet 'i' at at a time t+n∆t.

Is my strategy correct? If someone can explain to me how this problem can be solved in an organized manner I would really appreciate it, since my thoughts about how to solve this aren't really organized right now.

Answer 1

Use a computer.

Actually such calculations are routinely done by "numerical integration", which is just to say, exactly as you describe: Step by step.

Your thoughts are as organized as they can be, there is nothing more to it. Start with initial positions and velocities at time $t$. From those, calculate the force on each body, and thus its velocity and speed and the next time step $t+\Delta t$. Rinse and repeat. A computer does that without complaining.

When you have a result, you may want to confirm that your result does not depend on the step size that you used. That can be quite simple: re-do the same calculation but with smaller steps. If the result is the same (to within the precision that you care about), then you know that your numerical integration is accurate.

BTW, there is no analytic solution to this problem. Numerical methods are the only way to solving this.