# Integrate the law of universal gravitation to get position [closed]

I'm trying to write Physics simulation software which simulates gravitational and electric forces between $n$ bodies.

I know how constant acceleration is integrated to get distance, but I have no clue on how to do the same for a varying force (and acceleration). Let alone multiple forces.

Can someone guide me through this?

• This is going to be a differential equation. For an arbitrary number of n bodies, I'd recommend not so much solving it exactly, but letting the computer solve it approximately. You just have to write down your equations in terms of $x$, $\dot x$, $\ddot x$, $t$, and whatever else you have. Then you set up a numerical integrator that can handle systems of differential equations. Then you click run. Then you wait 6 years. Badda-bing badda-boom, you're done. Easy as that – Jim Jan 29 '15 at 22:16
• I'm actually writing the code. No "run" button. The question is probably not clear enough. I'm struggling with the physics (and maths), not with the computer part. – rrrrrrrrrrrrrrrr Jan 29 '15 at 22:29
• Would Computational Science be a better home for this question? – Qmechanic Jan 29 '15 at 23:02
• Probably not, since I'm really asking for the physics and maths here. I wouldn't be able to solve this with pen and paper either. Maybe I should rewrite the question. Sorry, I'm really new to StackExchange. – rrrrrrrrrrrrrrrr Jan 29 '15 at 23:03
• Guiding you through this is too much to ask here, chapters of textbooks are written about doing this problem. However, the basic idea can be found in this answer of mine. – Kyle Kanos Jan 30 '15 at 1:13

If you want the conceptual answer of what the equations would be for the exact answer, rather than practical numerical rules that approximate the results of the integral, the basic answer is that you first break down each force into components along each coordinate axis ($x$, $y$ and $z$ if you're using Cartesian coordinates), then you add the components of all the forces along any given axis to get the net force on that axis, say $F_x(t)$, as a function of time. Then you divide by the mass of the object being accelerated to get the acceleration along that axis as a function of time, $a_x(t)$, and velocity along that axis as a function of time will be the integral of that: $v_x (t) = \int_{t_0}^t a_x(t^{\prime}) \, dt^{\prime}$ (here $t^{\prime}$ is a dummy variable used inside the integral to distinguish it from the upper limit of integration $t$, $t_0$ is whatever time you want to start from, and also note that you would use the initial velocity $v_0$ at $t_0$ to determine the value of the constant of integration). Likewise, position along that axis as a function of time is $x(t) = \int_{t_0}^{t} v_x(t^{\prime}) \, dt^{\prime}$ (using the initial position $x_0$ at $t_0$ to determine the constant of integration). Then you do the same thing to derive $y(t)$ from $F_y(t)$ and $z(t)$ from $F_z(t)$. How difficult this is to solve exactly depends on the particular function $F_x(t)$, $F_y(t)$, $F_z(t)$ that you start with.

Also note that in many contexts you don't actually know the functions for forces as a function of time to begin with, but only get them by solving one or more differential equations telling you how the instantaneous force at any given moment relates to each physical object's position (and possibly velocity) at that moment. Often a numerical technique will be to start with such differential equations and then use some type of finite difference method to approximate a solution, see this pdf or this one for some discussion.

I'm trying to write Physics simulation software which simulates gravitational and electric forces between n bodies.

You can have a look at my code at http://yukterez.ist.org/4kp2, sorry if I won't post the whole 1 meter long code in Latex (:

Anyway, if you have multiple forces you just add them up as vectors. For instance here the movement of particle 1 with mass m1 and charge q1 on the x axis if attracted by the particles 2 und 3 with masses m2 und m3 and the charges q2 and q3:

$a1(t) = x1''(t) =$

$$\frac{G \text{m2} (\text{x2}(t)-\text{x1}(t))}{\sqrt{\left((\text{x2}(t)-\text{x1}(t))^2+(\text{y2}(t)-\text{y1}(t))^2+(\text{z2}(t)-\text{z1}(t))^2\right)^3}}+$$

$$\frac{G \text{m3} (\text{x3}(t)-\text{x1}(t))}{\sqrt{\left((\text{x3}(t)-\text{x1}(t))^2+(\text{y3}(t)-\text{y1}(t))^2+(\text{z3}(t)-\text{z1}(t))^2\right)^3}}+$$

$$-\frac{\text{q1} \text{q2} (\text{x2}(t)-\text{x1}(t))}{((4 \pi \text{\varepsilon 0}) \text{m1}) \sqrt{\left((\text{x2}(t)-\text{x1}(t))^2+(\text{y2}(t)-\text{y1}(t))^2+(\text{z2}(t)-\text{z1}(t))^2\right)^3}}+$$

$$-\frac{\text{q1} \text{q3} (\text{x3}(t)-\text{x1}(t))}{((4 \pi \text{\varepsilon 0}) \text{m1}) \sqrt{\left((\text{x3}(t)-\text{x1}(t))^2+(\text{y3}(t)-\text{y1}(t))^2+(\text{z3}(t)-\text{z1}(t))^2\right)^3}}$$

Then you repeat the same procedure for y and z axis, and for the other two particles. If you have more forces or particles just add them up.

Post Script:

But what I'm really asking is for the position as a function of time, not the acceleration.

Then just solve for the 0th derivative x(t) instead for the second one x''(t). Here is an example how you solve for acceleration, velocity and position after a given time (I choose 5 seconds, and for simplicity I reduce to 1 dimension x and 2 bodies):

Here the Objects are placed 1000 m ahead of each other on the x axis, one with mass 10000 kg and charge +1 Coulomb, the other with 5000 kg and -1 Coulomb. Initial velocity for example 10 m/s and -20 m/s, and solving for the properties of particle 1 after 5 seconds: So the acceleration after 5 sec is 1.363 m/s², the velocity 15.45 m/s and the position 62.69 m. Understand now?

If you want to solve for particle 2 replace the x1 by x2 in the last 3 evaluations. If there are y and z axes follow the described pattern; to combine them into a plot or simulation I choose this method: Link

Post post script:

I work with Python. How do I solve it numerically?

I don't have experience with Python, but since you say

I'm struggling with the physics (and maths), not with the computer part.

I'm sure you'll manage to translate that because the Syntax is self explaining.

If you want to upgrade the math even farther for solar system simulations in this PDF, page 7 you can find the corrections for asteroid belt, solar quadrupole moment and relativity. However, you always get the terms for acceleration as second derivative which has to be solved differentially.

• Thank you. That will be extremely useful by itself. I would really appreciate if you'd talk about the calculus behind that code, as my main problem is that I can't figure the integrations out by myself. I know I can add up forces, but I can't really grasp the part where the force become the position as a function of time. Which is the whole problem. :) – rrrrrrrrrrrrrrrr Jan 29 '15 at 22:58
• I really appreciate your edit. But what I'm really asking is for the position as a function of time, not the acceleration. That, I can figure out. :) – rrrrrrrrrrrrrrrr Jan 29 '15 at 23:07
• You can not integrate a n-Body simulation, you have to solve it numerically. Which program do you work with? If Mathematica use NDSolve instead of NIntegrate. – Gendergaga Jan 29 '15 at 23:08
• I work with Python. How do I solve it numerically? I just take a very short time, assume that the acceleration is constant though that time and calculate the displacement? – rrrrrrrrrrrrrrrr Jan 29 '15 at 23:27
• Right, if you solve it differentially then the "very short time" becomes infinitesimal. If extra precision is needed on Mathematica I would add the settings MaxSteps -> Infinity, Method -> {"ImplicitRungeKutta", "DifferenceOrder" -> 20}, InterpolationOrder -> All – Gendergaga Jan 30 '15 at 1:50