# General Formula For Gravitational Vector Field

## Recommended Posts

As if I could make title of this thread more daunting enough, the length of this post should scare off 90% of you I've buried two questions in this post, which I've highlighted bold to make easier to see.

I dont know why I've been thinking about it recently, but I was wondering about (Newtonian) gravity and the equations to calculate the movement of objects. I have a naive idea of how to do it:

- For every object in space, sum up all of the gravitational vector forces acting on it, "press" the force against your object to plot its new coordinates.

Basically, that does work, but its not very efficient. If you were trying to create a computer model of gravitation using that approach, it would bring your computer to its knees, because the amount of processing overhead increases at a rate of O( (x^x) * n) for each successive movement, where x is the number of planets you're trying to model and n number of actual moves you want to make. (At just 10 objects, you're computer has to go through 10 000 000 000 loop iterations for each step).

I've thought of a better way to model gravitation:

- Calculate the gravitational vector field as a series of time-based parametric equations, plot the coordinates with respect to time.

That works a lot more efficiently, you only have to calculate the vector field once, and then loop through all of your planets and plot their new coordinates. So the processing overhead is closer to O(xn), because you need to loop through your planets twice (once to get the values necessary to make your vector field, and once more to move your planets).

I understand the process of calculating a gravitational vector field, its just a very simple gradient of your force vector, where the force of gravity is given as:

F = Force of Gravity =

G * mass1 * mass2

---------------------

(mass1.x+mass2.x)^
+ (mass1.y+mass2.y)^2
+ (mass1.z+mass2.z)^2

Assuming mass1 is located at the origin,

Vector r = x2[i]i[/i] + y2[i]j[/i] + z2[i]k[/i]
Unit vector u = r / ||r||

G * mass1 * mass2
-----------------   *   u
||r||^2

Thats works nice, but I guess physics isnt my strong area, because that equation only works for 2 objects. What if I wanted to model 12 objects? I cant, because theres no such thing as a distance between 3 or more non-coplanar points. At this point, my calculus knowledge comes to a hault, because the simple vector field above is a little too simple, it doesnt work for more than 2 planets. How exactly do we generalize the gravitational vector field for more than 2 planets at any specific time t?

Ok, this is just a followup to the previous question, but I just noticed that its still not a very good model. It tells you the instantaneous vector field at any given time t, but the vector field itself is constantly shifting (because the masses are moving). To make the model above work in a visual demonstration, you need to create a loop that recalculates the vector field after each time you move your planets around, and that could be very resource intensive.

So, in order to have a better model, we need to account for all of the masses in the vector field moving, so the vector field itself is actually in fluid motion. If we do this, then we only need to calculate our vector field once, then we can run our model indefinitely with O(x+n) (linear) processing overhead. To model a fluid vector field, my guess is that we actually need a gradient of a gradient, which is actually very easy to do mathematically, but I'm clueless as to how you'd even set up this kind of equation. How do we generalize a gravitational vector field at all points in time?

##### Share on other sites

Firstly multibody equations are AMAZINGLY complicated when they are time dependent.

For non-time dependent gravity at anypoint is just a sum of all the bodies that you are taking into account.

$\bold {g} (\bold {r})_{total} = \sum_i -G \frac {m_i} {r_i^2}\bold {\hat r}$

To add a time dependence into this you have to make all of the r's time dependent, and they will ALL depending on every other objects gravity. So it becomes very complex, very quickly.

##### Share on other sites

A correction for my formula (thanks Atheist)

$\bold {g} (\bold {r})_{total} = \sum_i -G \frac {m_i} {(r-r_i)^2}\bold {\hat (r-r_i)}$

or:

$\vec g (\vec r) = -G \sum_i m_i \frac{ (\vec r - \vec r_i)}{\| \vec r - \vec r_i \| ^3}$

Which means the same just written in a differnt way... ri is the vector to the object being considered for each sum point, r is the vector to the point in space you are considering the field in.

And it's not the equations that are so complicated it's analytically solving them. Atheist says you could solve them numerically using O(x2). But I'm about to go to sleep so I need to think about this...

##### Share on other sites

How do we generalize a gravitational vector field at all points in time

The short answer is, you don't. The two-body problem can be solved analytically, but not when you have 3 or more. You can assume that the other bodies don't perturb the existing field, or exclude some bodies if you set limits on how much error you're willing to introduce into your answer, or you can do the really messy full-blown math you've started and Klaynos has described (with an apparent assist from Atheist) and solve them numerically

##### Share on other sites

And this is why accurately modelling a solar system with multiple bodies, over a large period of time, requires a complex program running on a supercomputer.

##### Share on other sites

If I was to write a computer program to simulate graverty of spherical masses (like planets) over time I would do it like this. This is probbebly a massive oversimplifacation but it works and is really fast. This is in 2D but it could be expanded into 3 dimentions with a bit more code.

Iv'e tried to write it in a semi code where most people should be able to understand it, however you still need to be able to understand more than '>', less than '<', Variables and Arrays.

grav(0) = 256
dist_max = 40
np = the number of planets
[color="DarkGreen"]Some code I havn't included to start the masses in starting positions
and optionally with starting velocitys, Its not there because it would
be diffrent depending on what you where trying to simulate[/color]
Loop (where n goes from 1 to dist_max)
grav(n) = grav(n - 1) * 0.5
End of Loop
Loop (where n goes from 1 to np)
x(n) = x(n) + fx(n)
y(n) = y(n) + fy(n)
[color="Navy"]Loop (where nn goes from 1 to np)
If (if n is not equal to nn then)
dist_x = x(n) - x(nn)
If dist_x < 0 then dist_x = -dist_x
dist_y = y(n) - y(nn)
If dist_y < 0 then dist_y = -dist_y
dist = dist_x + dist_y
dist = the integer of dist
dist2 = 0
If dist < dist_max then dist2 = grav(dist)
fx(n) = fx(n) + ((x(n) - x(nn)) * (dist2 * 0.0001))
fy(n) = fy(n) + ((y(n) - y(nn)) * (dist2 * 0.0001))
End of If
End of Loop[/color]
[color="DarkGreen"]Draw (x_old(n), y_old(n)) in the background colour
Draw (x(n), y(n)) in another colour[/color]
x_old(n) = x(n)
y_old(n) = y(n)
End of Loop

It could be easely modified to be points of negative graverty, Include friction etc...

obviously there is no friction in space.

Could this be converted to one LaTeX equation?

##### Share on other sites

- For every object in space, sum up all of the gravitational vector forces acting on it, "press" the force against your object to plot its new coordinates.

Basically, that does work, but its not very efficient. If you were trying to create a computer model of gravitation using that approach, it would bring your computer to its knees, because the amount of processing overhead increases at a rate of O( (x^x) * n) for each successive movement, where x is the number of planets you're trying to model and n number of actual moves you want to make. (At just 10 objects, you're computer has to go through 10 000 000 000 loop iterations for each step).

Hi there,

I know this doesn't answer your question, but actually you can make this method work very efficiently. All it takes is a little extra effort when optimizing your code. I've attached some very simple example programs to prove my point. I don't know how many loops/step they spend, but it's definitely much less than the number you mention above.

The first example has five hundred bodies that all interact, and it should run nicely on any newer pc.

If you're simulating a system with a lot of bodies with neglible mass you can make it even faster:

The second example has five thousand bodies, of which 10 have masses high enough to take into consideration.

If anyone's interested in seeing some code, Id be glad to post it and explain how it works. These examples are simple and only use 2 dimenions, but making it run in three dimensions will only make it a tad slower. I've used the same code to simulate some of the planets in our own system, and it returns some pretty accurate data. In other words there is no loss of precision due to high performance.

The programs are both approx. 100 lines, and I've written them in FreeBasic.

Best regards,

Michael

Gravity_example_1.zip

Gravity_example_2.zip

##### Share on other sites

Is it simerlar to my program? If not I would love to know how it works.

##### Share on other sites

The funny thing about simulating such large systems is that the solutions are chaotic. The qualitative behavior of the systems (e.g. spectral properties of positions of function of time, shape of orbits) can be reproduced very well with such simulation software, even by using simple first order approximations for the next time step. The precise quantitative behavior as function of time (exact position of all objects at a certain point in time) is amazingly hard to reproduce. Luckily, when one wants to study the behavior of such systems, then the qualitative behavior is sufficient.

The fact that precise quantitative behavior is very hard to reproduce tells something about the real system. If such a system exists in reality, then the tiniest perturbation by means of an external (not modelled) force also screws up the precise quantitative behavior, while hardly affecting the qualitative behavior.

What I told here is not relevant for simulating relatively simple systems with a few tens of objects (e.g. our solar system with the planets and the heaviest asteroids and several moons), but it becomes important for simulating systems of 1000's or even 1000000's of objects. I did experiments with this (I have written a program, using DASSL and later DASPK) and have found that this behavior in fact is quite common for many natural phenomena.

##### Share on other sites

Is it simerlar to my program? If not I would love to know how it works.

Hi alan2here,

I'm afraid there's quite a bit of difference between you code snippet and my programs. Actually it does'nt look like your program is going to work at all - did you test it?

Anyway, I'm simply too tired righ now to explain how my programs works, but here's the important part of the code for you to look at. I'll get back to the explanation soon, but It'd be good practice for you to try and figure it out yourself

  For i = Lbound(Body) to Ubound(Body)

For i2 = i+1 To Ubound(Body)

Dist_X = Body(i).X-Body(i2).X
Dist_Y = Body(i).Y-Body(i2).Y
Dist_sqared = Dist_X^2 + Dist_Y^2
Dist = sqr(Dist_Sqared)

Force = Body(i2).Mass/Dist_sqared
Body(i).Xvec -= (Dist_X/Dist)*Force
Body(i).Yvec -= (Dist_Y/Dist)*Force

Force = Body(i).Mass/Dist_sqared
Body(i2).Xvec += (Dist_X/Dist)*Force
Body(i2).Yvec += (Dist_Y/Dist)*Force

Next

Body(i).X += Body(i).XVec
Body(i).Y += Body(i).YVec

Next

Edit:

If you want to tamper with the whole piece of code, I'll post it for you. I'd recommend downloading the FreeBasic and FBide bundle here. It's a great and really fast basic compiler and it's very easy to use. Anyway I had lots of fun with it, and I never really bothered learning C++ since freebasic is pretty much just as fast and versatile.

Best regards,

Michael

##### Share on other sites

Sorry I forgot to say that the code starting after the first green bit and ending at the end should be triggered once every frame

My code works verry well. Maybe it just does it in a diffrent way to yours, I am trying to replicate you'rs now, please try and replicate mine and you will see.

I am using VB 6.0 btw

One diffrence is that I am drawing using a command that draws onto the window and using two variable arrays x(planet number) and y(planey number), you are using an object that has .x and .y built into it already.

Is LBound the name of an array or a command?

What the heck is .Xvec, is it x cordatate vector? If it means vector then are vectors a data structure e.g 'an array' or a matmatical teqneek e.g 'modulo'?

Why do you have -= and +=? is this a languege spacific like how 'c' used ++ to increment a variable.

We both are finding the difrence between each planet to every other planet by taking one away from the other

You are doing it in this line

Dist_X = body_x(i) - body_x(i2)

Dist_Y = body_y(i) - body_y(i2)

my FX and FY are force on X cordanate and force on Y cordanate

for example If whitin the main loop of the program I put this

fx(n) = fx(n) * 0.9

fy(n) = fy(n) * 0.9

Then it would act as friction slowing everything down and causing all the masses to spiral into one another.

fx(n) = fx(n) + 0.01 speeds up planet n by 1 in the x (across) direction.

The program basicly finds the difrence between each planet and each other planet, simerly to your's, your's seems to use a lot less code which is good and mine works out (because of how it does it) the distance as if planets can only be an integer value of pixels appart from each other.

You'rs also seems to turn into a spinning disk quite quickly and contains planets of diffrent weights and sizes, wheras mine will become more messy and cayotic if large numbers of planets are started at random positions and velocitys and in mine all planets weigh and are the same size and drawn as single points.

Sorry I forgot to say that the code starting after the first green bit and ending at the end should be triggered once every frame

My code works verry well. Maybe it just does it in a diffrent way to yours, I am trying to replicate you'rs now, please try and replicate mine and you will see.

I am using VB 6.0 btw

One diffrence is that I am drawing using a command that draws onto the window and using two variable arrays x(planet number) and y(planey number), you are using an object that has .x and .y built into it already.

Is LBound the name of an array or a command?

What the heck is .Xvec, is it x cordatate vector? If it means vector then are vectors a data structure e.g 'an array' or a matmatical teqneek e.g 'modulo'?

Why do you have -= and +=? is this a languege spacific like how 'c' used ++ to increment a variable.

We both are finding the difrence between each planet to every other planet by taking one away from the other

You are doing it in this line

Dist_X = body_x(i) - body_x(i2)

Dist_Y = body_y(i) - body_y(i2)

my FX and FY are force on X cordanate and force on Y cordanate

for example If whitin the main loop of the program I put this

fx(n) = fx(n) * 0.9

fy(n) = fy(n) * 0.9

Then it would act as friction slowing everything down and causing all the masses to spiral into one another.

fx(n) = fx(n) + 0.01 speeds up planet n by 1 in the x (across) direction.

The program basicly finds the difrence between each planet and each other planet, simerly to your's, your's seems to use a lot less code which is good and mine works out (because of how it does it) the distance as if planets can only be an integer value of pixels appart from each other.

You'rs also seems to turn into a spinning disk quite quickly and contains planets of diffrent weights and sizes, wheras mine will become more messy and cayotic if large numbers of planets are started at random positions and velocitys and in mine all planets weigh and are the same size and drawn as single points.

Here is the Program

http://www.with-logic.co.uk/a/Grav.exe

Here are the VB Files

http://www.with-logic.co.uk/a/Project1.vbp

http://www.with-logic.co.uk/a/Project1.vbw

http://www.with-logic.co.uk/a/Form1.frm

You will see that single points moving quickly towards a mass of points are effected by the compound force of all the points over larger distances than just one on one, as you would expect in a realistic simulation.

##### Share on other sites

My code works verry well. Maybe it just does it in a diffrent way to yours, I am trying to replicate you'rs now, please try and replicate mine and you will see

Indeed it works fine. I'll make my own version, too - that'd be fun! Gravity is an "inverse sqare" force, and since I didn't seem to find any inverse sqare stuff in your code I though it didn't work - but now I see it!

I must say that my program is probably a lot more accurate than yours, since your code only have 40 different values of gravitational attraction depending on distance. Since gravity depends on both mass and distance, this limits your program to only simulate bodies that all have the same mass.

Is LBound the name of an array or a command?

If I dim an array named "A" from 0 to 100, then Lbound(A) refers to the lowest number in the array = 0. Opposite is Ubound(A) = 100.

What the heck is .Xvec, is it x cordatate vector? If it means vector then are vectors a data structure e.g 'an array' or a matmatical teqneek e.g 'modulo'?

I'm using the so called types, which is a sort of variable with any number of "sub-variables" you like - like this:

Type Planet
X As Single
Y As Single
XVec As Single
YVec As Single
Mass As Single
Col as Integer
End Type 

All values are dimmed as singles = single precision floaing point. X and Y are the cordinates of the planet. Xvec and Yvec are the horisontal and vertical velocity vector. Mass, radius and color are self-explanatory.

So, when it says

Body(i).X += Body(i).XVec

this simply means add the horisontal velocity vector to the x coordinate of planet number i - this then defines the new x value, in which the planet will be drawn in the next loop.

Why do you have -= and +=? is this a languege spacific like how 'c' used ++ to increment a variable.

X += 5 is just a short form of writing X = X + 5. A lot of programming languages allows you to use this short form.

You'rs also seems to turn into a spinning disk quite quickly and contains planets of diffrent weights and sizes, wheras mine will become more messy and cayotic if large numbers of planets are started at random positions and velocitys and in mine all planets weigh and are the same size and drawn as single points.

The planets in my program doesn't shape a spinning disc by it self - I make them do that. I place a large object (the sun) in the center, then I place all objects in a circle around it, and then I "push" them exactly enough to obtain a circular orbit around the sun. I wouldn't mind explaining it in greater detail, if you want that.

Search wikipedia for all the explanations and formulas. It's a great resource if you want to know more about this stuff.

Michael

##### Share on other sites

Got it!

Btw there is an error in your code. When finding the correct distance using pythagoras your code shouldn't look like this:

        dist_x = x(n) - x(nn)
If dist_x < 0 Then dist_x = -dist_x
dist_y = y(n) - y(nn)
If dist_y < 0 Then dist_y = -dist_y
dist = dist_x + dist_y
dist = Int(dist)

Instead it should look like this:

        dist_x = x(n) - x(nn)
dist_y = y(n) - y(nn)
dist = Sqr(dist_x^2 + dist_y^2)
dist = Int(dist)


Since the hypothenuse of the triangle - your distance - is sqare root ( dist x sqared + dist y sqared), not just dist x + dist y.

I also cut out this:

If dist_x < 0 Then dist_x = -dist_x
If dist_y < 0 Then dist_y = -dist_y

since -X sqared = X sqared and of course -Y sqared = Y sqared

Oh, and please notice how much faster the freeBasic version of the program is running. Mind you, it is an exact copy of your program - I haven't changed a thing! (except fixing the distance error mentioned before, and that actually slows it down a bit.)

This is why you should definitely drop visual basic and go FreeBasic! The syntax is almost the same, but the compiler is much, much better. The newest version also supports object oriented programming and whatnot.

Edit: the attached file contains all programs + source code discussed so far + a new version I made of your program. You can more or less just copy - paste my code into VB since the languages are much alike.

Michael

Gravity programs.zip

##### Share on other sites

Thank you for the optimisations and corections, they have been verry helpfull.

I didn't know I was using pythagoras, i'll take it that it is more accurate now than how I was doing it.

Im looking into FreeBasic now, If I wasn't going to be lurning BlueJ (Java) and "VB 2005" next turm in college having just lurned FPC (pascell) this turm I would defenatly be downloading it.

##### Share on other sites

you can optimise the program down to the point where you do (n^2 + n)/2 iterations per interval. where n=the number of significant feilds

it's basically the "handshake problem" how many handshakes if everyone in a room were to shake hands. it'll end up as triangular numbers.

you can cut a few more lines out by writing

distance (d) = root(x^2+y^2+z^2)

force(x) = x/(distance^3)

because you're essentailly writing

x component of the unit vector = x/d

force = 1/(d^2)

x component of force = force * x component of unit vector

so it follows that

1/(d^2) * x/d = x/(d^3)

what do you do to prevent the objects superimposing? i often get readings like d=0 so i get undefined values for force. should i add 1 to the distance before i determine the force?

##### Share on other sites

Your'e welcome, alan2here. It's been fun, and I*m glad you showed me your version of the program, too. I'll definitely keep working on that idea!

you can optimise the program down to the point where you do (n^2 + n)/2 iterations per interval. where n=the number of significant feilds

it's basically the "handshake problem" how many handshakes if everyone in a room were to shake hands. it'll end up as triangular numbers.

Yup, didn't get around to mentioning that. This is what I do when coding:

For n = 1 To np
For nn = n+1 To np
...
...
Next
...
Next

For n = 1 To np
For nn = 1 To np
...
...
Next
...
Next`

Oh, and its (n^2 - n)/2 ;-)

It can also be written as (n*(n-1))/2

you can cut a few more lines out by writing

distance (d) = root(x^2+y^2+z^2)

force(x) = x/(distance^3)

because you're essentailly writing

x component of the unit vector = x/d

force = 1/(d^2)

x component of force = force * x component of unit vector

so it follows that

1/(d^2) * x/d = x/(d^3)

Yes, I did this too when trying to optimize - but it slowed the program down! Go figure ?!? I suppose the ^3 part somehow means more calculations, but I honestly dont know.

what do you do to prevent the objects superimposing? i often get readings like d=0 so i get undefined values for force. should i add 1 to the distance before i determine the force?

I use a code snippet that prevents the "distance" variable to never be smaller than the sum of the two objects' radiuses. But its really a pain... You prevent d ~ 0 related problems by doing so, but the objects stop moving in realistic orbits when they overlap. Alan2heres program actually gets aroun this very nicely, since al gravitational attraction is precalculated and can never exceed the initially given max value. Nice and simple, although a tad inaccurate.

best regards,

Michael

##### Share on other sites

I use a code snippet that prevents the "distance" variable to never be smaller than the sum of the two objects' radiuses. But its really a pain... You prevent d ~ 0 related problems by doing so, but the objects stop moving in realistic orbits when they overlap. Alan2heres program actually gets aroun this very nicely, since al gravitational attraction is precalculated and can never exceed the initially given max value. Nice and simple, although a tad inaccurate.

hmm, so it's the uncertainty principle all over again... i like it.

i just programmed this in vb6, i used a fudge factor added to the distance, if the two objects are actually superimposing it reads force as if the distance were equal to the fudge factor. the only extra code is " + factor" on the end of the calculation for the distance

we havent really answered the questions on this thread, "can you do one equation to map the gravity gradient?" the answer is yes BUT it's a 4 dimesional equation and must be recalculated for each position and rewritten for each interval. it is the sum of the forces applied by each other particle felt differently in all locations.

i'm thinking about doing this 3d, any tips on simple foreshortenning or stereoscopics?

##### Share on other sites

we havent really answered the questions on this thread, "can you do one equation to map the gravity gradient?" the answer is yes BUT it's a 4 dimesional equation and must be recalculated for each position and rewritten for each interval. it is the sum of the forces applied by each other particle felt differently in all locations.

Yeah, as already mentioned in this thread we'll probably never see the n-body problem for n > 3 compressed into a single equation. And even IF we could get an exact quantitative description of mutiple-body systems, it'd still only give those exact results if it were to simulate a system of perfectly spherical and homogenous bodies that didn't change mass or density over time or even rotated around their axes - simulating a real system over long stretches of time would still not work, since a lot of different phenomenons like continuous destruction and formation of bodies, tidal effects and solar wind wouldn't be taken into consideration within such an equation. If we'd just wanted a qualitative simulation then that'd be ok - but in an exact quantitative simulation that'd change success into miserable failure.

i'm thinking about doing this 3d, any tips on simple foreshortenning or stereoscopics?

Hmm... I'm not entirely sure what you mean (englis isn't my mother tongue). If you're searching for a way to represent the result graphically, then three windows showing XZ - XY - ZY works pretty well. Some time back I saw a very nice and short openGL code snippet that allowed you to ghost around in a 3D representation of a system of moving bodies. It looked really nice, but I don't know if VB6 supports openGL? Actually I might want to look into that myself too - I'll post a link once I dig it up.

Edit: There you go. Scroll right down to the bottom. http://www.freebasic.net/forum/viewtopic.php?t=1664

##### Share on other sites

Hmm... I'm not entirely sure what you mean (englis isn't my mother tongue). If you're searching for a way to represent the result graphically, then three windows showing XZ - XY - ZY works pretty well. Some time back I saw a very nice and short openGL code snippet that allowed you to ghost around in a 3D representation of a system of moving bodies. It looked really nice, but I don't know if VB6 supports openGL? Actually I might want to look into that myself too - I'll post a link once I dig it up.

i mean perspective, like when you look down a tunnel and the back end appears smaller.

stereoscopics basically means projecting the image differently for each eye so you actually look into the screen, much the same as 3d movies, i'm working on the algorythm now. i'll try to go with just looking down the z axis then maybe allow the user to look in different directions.

##### Share on other sites

Sounds interesting, please drop a hint when you've got something. Anyway, this IS what graphics libraries like openGL and Direct3D were made for. But of course it takes time to learn - just started out myself - but on the other hand there are loads of ressources, tutorials, forums and the like out there!

##### Share on other sites

^^ I can recommend the NeHe OpenGL tutorials. You don´t have to read them in order; try getting the rotating cube to work and then stick your nose in whichever tutorial topic seems interesting and experiment around with it.

##### Share on other sites

i've got the algorythms working so far, with red/green glasses you actually perceive depth.

it's basically a right triangle using the z and x axes to place a dot where the hypotenuse intersects the screen. taken from two vantage points, depth perception becomes childs play.

it's not a mobile perspective, you can't move around in the frame. i'll take a look at openGL, sounds interesting.

## Create an account

Register a new account