This is why mathematics makes my head hurt... - by catbarf
*Zaccheus* on 17/10/2007 at 17:16
Ok.
:)
dvrabel, I'm not convinced I am using a the Euler method (which seems to be about integration, i.e. relatively simple regular curves) and I cannot see how I could use the 4th-order Runge-Kutta method.
Here's the algorithm in pseudo code (without the existing optimisations):
Quote:
/* The following calculations will be done for every combination of all the simulated objects: */ FOREACH thisObj IN objects
BEGIN LOOP
FOREACH thatObj IN objects
BEGIN LOOP
IF thisObj == thatObj CONTINUE LOOP /* objects don't accelerate themselves. */
/* Make a directional vector pointing from 'thisObj' to 'thatObj' */Vector delta = thatObj.location - thisObj.location
/* Find the length of the vector, i.e. the distance between the objects */Real d = length_of(delta)
/* Make a normalised vector pointing towards the source of the acceleration. */Vector nv = delta / d
/* Calculate the amount of acceleration (G is newton's constant of gravity). */Real absAcc = (G * thatObj.mass) / (d * d)
/* Extend the directional vector to be 'as long as' the amount of calculated acceleration. */
Vector acc = nv * absAcc
/* Add the directional acceleration to the object's directional acceleration. */thisObj.acceleration += acc
END LOOP
END LOOP
FOREACH thisObj IN objects
BEGIN LOOP
thisObj.velocity += thisObj.acceleration
thisObj.location += thisObj.velocity
thisObj.acceleration = (0,0,0)
END LOOP
The above is repeated for each simulated second.
dvrabel on 17/10/2007 at 22:25
You're correct in that you're not using the Euler method. But that's because you have a mistake. You have:
s<sub>n+1</sub> = s<sub>n</sub> + hv<sub>n+1</sub> = s<sub>n</sub> + hv<sub>n</sub> + h<sup>2</sup>a<sub>n</sub>
where h is time interval, s is displacement/position, v is velocity, a is acceleration.
For the Euler method you would have:
s<sub>n+1</sub> = s<sub>n</sub> + hv<sub>n</sub>
However, your method is almost a correct improvement of Euler. Try (from the Taylor series):
s<sub>n+1</sub> = s<sub>n</sub> +hv<sub>n</sub> + h<sup>2</sup>v'<sub>n</sub>/2 = s<sub>n</sub> + h<sub>n</sub> + h<sup>2</sup>a<sub>n</sub>/2
(Which may look familiar as one of the standard formula for constant acceleration motion, s = ut + 1/2.at<sup>2</sup>.)
Marecki on 17/10/2007 at 23:54
Try the leap-frog algorithm, Zacch. I've used this a lot in my old work of this sort, namely for molecular dynamics simulations with a large number of participants, and for this kind of stuff it was highly effective. Given the nature of your simulation, it might work for you too. Here is some reading material: (
http://phycomp.technion.ac.il/~david/thesis/node34.html)
Pyrian on 17/10/2007 at 23:57
Quote Posted by *Zaccheus*
So my question is, why would gravity be able to stop that in the future when it cannot stop it now.
I think someone's conflating the old theory with the new findings. In the old model, gravity was slowing down the expansion of the universe; the question was whether or not the ratio of gravity to expansion velocity was going to cause the universe to keep flying apart or eventually crunch back together. Turns out, the expansion of space is accelerating - nobody's really sure why. There's really no good reason to think that gravity alone is going to somehow counteract
that based on current evidence.
Quote Posted by *Zaccheus*
However, as I understand it, matter is not thought to have an outward velocity which gravity can slow down; rather space being created everywhere causes matter to end up further apart.
Well, space and mass are directly connected in general relativity. It's worth noting that the current acceleration of expansion was
not an expected finding, and nobody's really sure what's causing it ("Dark Energy").
Quote Posted by *Zaccheus*
That would imply that space itself has a velocity. This makes little sense to me.
Heh, it may not make sense to you but it's the analogy Einstein used to explain general relativity. Gravity does not directly exert a force - it is instead a result of the acceleration of space which objects are in.
*Zaccheus* on 18/10/2007 at 11:23
Marecki, thanks for that, I'll have a read. :)
Quote Posted by dvrabel
You're correct in that you're not using the Euler method. But that's because you have a mistake. You have:
s(n+1) = s(n) + hv(n+1) = s(n) + hv(n) + h^2.a(n)
For the Euler method you would have:
s(n+1) = s(n) + hv(n)
However, your method is almost a correct improvement of Euler. Try (from the Taylor series):
s(n+1) = s(n) +hv(n) + h^2.v'(n)/2 = s(n) + hv(n) + h^2.a(n)/2
(Which may look familiar as one of the standard formula for constant acceleration motion, s = ut + 1/2.at^2.)
That looks great ... what does it mean (what are "s", "h", "v", ".", and so forth) ?
By the way, the one error I do know about is this:
As you can see I calculate the acceleration over a certain time period, eg. one second. I model that as a single event, which seems wrong to me: During the one second of acceleration, the velocity will be increasing slowly and will reach its new value at the end of the second. But in my code I am adding the whole acceleration to the velocity as if the new velocity had been reached at the beginning of the second. That is wrong. Also, as the object starts moving it will get closer to the object which is accelerating it, therefore the acceleration in the second half of the second will be greater because the distance has already been reduced. But in my code I am calculating the acceleration as if the object had been at the starting distance over the whole second. That is wrong.
PeeperStorm on 21/10/2007 at 16:17
Quote Posted by *Zaccheus*
While we are on the topic of flight, I have a question:
If you take a table spoon and hold its outward curved side against a stream of running water coming out of a tap, the spoon will be pulled into the stream of water. Is that the same effect which keep a plane in the air?
Yes, it's that doofus Bernoulli and his annoying (
http://en.wikipedia.org/wiki/Bernoulli_effect) principle again.
Marecki on 21/10/2007 at 23:14
Quote Posted by *Zaccheus*
Marecki, thanks for that, I'll have a read. :)
You're welcome, always glad to help a fellow geek :)
Quote:
That looks great ... what does it mean (what are "s", "h", "v", ".", and so forth) ?
Quantities: s - position, h - time difference between two steps (i.e. the "length" of each step), v - velocity, v' (= a) - acceleration.
Operators: . - multiply, ^2 - squared.
*Zaccheus* on 22/10/2007 at 11:41
Thanks.
:)
Quote:
s(n+1) = s(n) + hv(n) + h^2.
a(n)/2That should correct the known error I was talking about when I said (
http://www.ttlg.com/forums/showthread.php?p=1658311#post1658311) "adding the whole acceleration to the velocity" is wrong because the object would not have been travelling at the new velocity during the whole time interval.
I had actually tried something very similar myself; it caused problems with objects in circular orbits. However, I think I was calculating the initial orbital velocity the wrong way back then (using trigonometry!), so I'll have to go back and try it again. One thing to be aware of is that while the new location only has to include half the acceleration, the new velocity has to include the whole acceleration.
dvrabel on 22/10/2007 at 12:39
Note that all these numerical methods will have errors since you're approximating a smooth curve with small curves (or straight lines) that don't quite match.
For example with:
s<sub>n+1</sub> = s<sub>n</sub> + hv<sub>n</sub> + h<sup>2</sup>a<sub>n</sub>/2
You're still assuming that the acceleration is constant during each time period which is generally not the case (except for objects in perfectly circular orbits). But given sufficiently small h the errors are not significant.
To look at this sort of problem in more a rigorous, mathematical way requires knowledge of maths at the A-level Further Mathematics level (simple differential equations, Taylor series and basic numerical methods).
*Zaccheus* on 22/10/2007 at 14:18
Thanks again for confirming
a<sub>n</sub>/2, I'll see if I can get that working in my code.
The changing acceleration is almost impossible to know because the objects are all (
http://www.rclsoftware.org.uk/images/gravel/chaotic.png) moving during the time interval, which in my code is always
one. Even in circular orbits the direction of the acceleration is constantly changing.
I have a degree in computing, but I've been thinking about doing A-level maths in evening class.
---
Nope, still have the same massive deterioration I saw when I first tried it ages ago.
The orbit is a lot less circular if I use
a<sub>n</sub>/2:
Quote:
FOREACH thisObj IN objects
BEGIN LOOP
thisObj.location += thisObj.velocity + (thisObj.acceleration / 2)
thisObj.velocity += thisObj.acceleration
thisObj.acceleration = (0,0,0)
END LOOP
I agree that
a<sub>n</sub>/2 should give me a better approximation, but for some reason it doesn't.
This is how I calculate the initial speed, assuming the central mass is much larger than the satellite's mass.
Quote:
Sqrt( (ConstantOfGravity * centralMass) / orbitRadius )
