Could someone pseudo-code it out for me.

First, I'll share my non-energy conservation system that I used in my
flocking algorithm ( http://www.geocities.com/pterandon/boids.html ), and in
planetary orbiting simulations.

dt=1  (we're stepping along povray frame by povray frame)
p= vector for current position
v=velocity
a=acceleration, or sum of forces at current position (based on gravity of
sun, "repulsion" of neighbors in flock),


----- for each frame----
p2=p1+v1 *dt
a2= constant * (sum of new forces based on new location).
v2=v1+a2
--- repeat ad nauseum ---


As I've shown before, (
http://news.povray.org/povray.binaries.images/thread/%3C3A65E429.18F2DC62%40my-dejanews.com%3E/?ttop=184334&toff=4400
)

it's impossible to have a "stable orbit" with this because this model
doesn't conserve energy.  I could however probably do a 1-D case for energy
conservation quite easily.

So, in povray, how would one pseudo-code out a model that conserves energy
and accounts for multiple fields (gravity of sun *plus* other planets,
repulsion of *all* neighbors in a flock).

----- for each frame----

p2=p1+v1 *dt
a2= constant * (sum of new forces based on new location).
v2=v1+a2
KE2=constant(vlength(v2)^2)

... I'm stumped.

Post a reply to this message

> So, in povray, how would one pseudo-code out a model that conserves energy
> and accounts for multiple fields (gravity of sun *plus* other planets,
> repulsion of *all* neighbors in a flock).
>
> ----- for each frame----
>
> p2=p1+v1 *dt
> a2= constant * (sum of new forces based on new location).
> v2=v1+a2
> KE2=constant(vlength(v2)^2)
>
> ... I'm stumped.

I've don't have a real clue on how you might solve this. Some interesting
links while googling where these:

About Energy as a whole, very straightforward and simple:
<http://rabi.phys.virginia.edu/105/2003/ps2s.html>

About some ballmenu for a Website, but right at the end, there's something
about energy conservation:
http://hep.itp.tuwien.ac.at/~ipp/aboutballmenuphysics.html

I'm not all too deep into these energy-theories etc and their scripting. But
if I'm not mistaken, you should just look to it that when two particles
exchange forces, their sum is zero. So, for example, you might just take all
movement vectors, and add them up. If the particles started from a resting
state, then they should add up to zero (unless there are things like
stationary particles or so). The resulting vector is the extra amount of
energy added to the system, and needs to be diminished. In turn, you might
just add a fraction of that vector to every particle, so that, in the end,
you get zero total energy.

But I'm really quiet fuzzy about this all and would just experiment with the
above technique to see if it works. If it doesn't, then google is back! :-)

-- 
"Tim Nikias v2.0"
Homepage: <http://www.nolights.de>
Email: tim.nikias (@) nolights.de

Post a reply to this message

Among other things, Greg M. Johnson wrote:

> Could someone pseudo-code it out for me.
> 
> First, I'll share my non-energy conservation system that I used in my
> flocking algorithm ( http://www.geocities.com/pterandon/boids.html ), and
> in planetary orbiting simulations.
> 
> dt=1  (we're stepping along povray frame by povray frame)
> p= vector for current position
> v=velocity
> a=acceleration, or sum of forces at current position (based on gravity of
> sun, "repulsion" of neighbors in flock),
> 
> ----- for each frame----
> p2=p1+v1 *dt
> a2= constant * (sum of new forces based on new location).
> v2=v1+a2
> --- repeat ad nauseum ---
> 

Hmm... I thought I had replied to this... Was my post lost somehow? I can't 
find it anywhere :/

Well, I'll write it again.

First, shouldn't it be:

  v2=v1+a2*dt ?

Second, as I mentioned, there are several algorithms used for integrating 
the equations of motion (that's what you're trying). An algorithm that 
minimizes numerical errors is the so-called "velocity Verlet" [Swope et al. 
J. Chem. Phys. 76 (1982) 637]:

-------------
  Calculate all forces on all particles, for each particle:
    a1 = (total force at p1)/(real or fictitious mass)
1 Calculate new positions for all particles:
    p2 = p1 + v1*dt + 1/2*a1*(dt^2)
  Calculate intermediate velocities (at t1+1/2*dt) for all particles:
    v' = v1 + 1/2*a1*dt
  Calculate all forces on all particles at their new position:
    a2 = (total force at p2)/(real or fictitious mass)
  Caluclate the new velocities for all particles:
    v2 = v' + 1/2*a2*dt
  Goto 1 (p2 -> p1, v2 -> v1, a2 -> a1)
-------------

This should conserve energy. What is not conserved is kinetic energy 
(temperature). If a simulation at constant temperature is desired, there 
are some "tricks" available, the simplest one is just re-scaling all 
velocities so that kinetic energy remains constant. After calculating the 
new velocities, do:

  f = sqrt(2*Ek/sum(mass*v2^2))
  v2(scaled) = f*v2

where Ek is the desired constant kinetic energy.

-- 
light_source{9+9*x,1}camera{orthographic look_at(1-y)/4angle 30location
9/4-z*4}light_source{-9*z,1}union{box{.9-z.1+x clipped_by{plane{2+y-4*x
0}}}box{z-y-.1.1+z}box{-.1.1+x}box{.1z-.1}pigment{rgb<.8.2,1>}}//Jellby

Post a reply to this message

<"Greg M. Johnson" <gregj;-)565### [at] aolcom>> wrote:

> Could someone pseudo-code it out for me.

AFAIK there is no simple "energy conservation algorithm".
You and some others around here seem to have the misconception
that the energy of each particle is conserved.
However this is only the case if you only have one central 
force (the sun for example), but no additional forces between 
your particles (but in this case you would need no numerical
integration at all).
If you have 3 or more particles interacting with forces, only the
sum of the energy of all particles is conserved. So there is no easy
way to use it directly.

Using a better integration algorithm, like runge-kutta, should solve
your problem.

Lutz-Peter

Post a reply to this message

First of all thanks Jelby for that link in the other thread.  I went to it,
but it was a program in some other language.  I'm one of those idiots who's
not only unable to program in anything in but povray but also
philosophically chooses to do number crunching in the program. (Yeah, I
know, povray is a raytracer, not a physics software language [take it easy,
McCoy  ]).

I think some of you were answering the question about the conservation of
energy after collisions  based on an entire system of particles, or after a
three-way collision.  Thanks, but that was not actually the problem of my
interest.

When I tried two years ago, I wasn't able to set up a correct model of the
solar system, (perhaps it might help to say to do an "AI" model?),  without
the planets either sinking into the sun or leaving the system altogether.  I
think I tried pretty hard, and i could not find a gravitational constant (or
a sufficiently small time slice) where particles didn't do one or the other.
I was wondering if it were the case that one cannot do it with a computer.
By that I mean that as a travel any finite step obliquely across a gravity
well,  there's a difference between real life and my model.  My model just
re-computes new force at position 2,  whereas a real-life particle would
have the gravity acting on it both "instantaneously" and "continuously".
Someone noted way back then that I wasn't conserving energy. I think others
with nice intentions pointed me to theories that were over my head, and I
suspect, not applicable.

Then I saw the concern in Tim's thread about conserving energy and wondered
if there were some trick.


Jelby, in your second post,  you offer an interesting idea.
>  p2 = p1 + v1*dt + 1/2*a1*(dt^2)
But is fundamentally any different from going to a smaller time slice?  And
I was hoping for a solution where dt= one povray frame.

So I'm falling back to three ideas:
1) it is in fact impossible

2) there's some way to do this where you do an actual integration of the
forces as one moves from p1 to p1+v1dt.   I hope I'm not too dense and
haven't missed Jelby explaining exactly how to do this.

3) Perhaps there was an algorithmic error on my part. In my earlier attempt,
http://news.povray.org/povray.binaries.scene-files/thread/%3C3A65E45C.907F6BB7%40my-dejanews.com%3E/?ttop=184097&toff=400

                #declare actorp[n]=actorp[n]+.05*actorv[n];
                #declare actorv[n]=actorv[n]+Grav*(actoravoid[n]);


I was only moving the particles 5 percent of their velocity, because I'd
seen that it made for a smoother curve, and in my head I was thinking it was
like miles/s and the frames were just seconds.  But perhaps I lost
responsiveness (and gained energy loss) from this technique.    I'm doing
some experiments now with this set at 1.0

Thanks for all your consideration and interest.

Post a reply to this message

In article <4041693c$1@news.povray.org>,
 "Greg M. Johnson" <gregj;-)565### [at] aolcom> wrote:

> I think some of you were answering the question about the conservation of
> energy after collisions  based on an entire system of particles, or after a
> three-way collision.  Thanks, but that was not actually the problem of my
> interest.

Any 3-way interaction, actually. Even just gravitational interaction.


> By that I mean that as a travel any finite step obliquely across a gravity
> well,  there's a difference between real life and my model.  My model just
> re-computes new force at position 2,  whereas a real-life particle would
> have the gravity acting on it both "instantaneously" and "continuously".
> Someone noted way back then that I wasn't conserving energy. I think others
> with nice intentions pointed me to theories that were over my head, and I
> suspect, not applicable.

There is no way to compute an exact solution for 3 or more 
gravitationally interacting bodies. However, you can compute a very 
close approximation. I did a simulation (in C++) that was able to 
simulate the orbits of the planets from Mercury to Jupiter for several 
Jovian years (several dozen Earth years), with very little drift. (I 
simply haven't entered the data for the other planets)
There are some rather complex algorithms needed to do this with both 
accuracy and efficiency, however. And there will always be cumulative 
errors, the simulation will diverge from reality over time, there is no 
way around that.


> Jelby, in your second post,  you offer an interesting idea.
> >  p2 = p1 + v1*dt + 1/2*a1*(dt^2)
> But is fundamentally any different from going to a smaller time slice?  And

This basically assumes that acceleration will be constant over the time 
interval, and takes that into account for the new point. I take it you 
were originally using p2 = p1 + v1*dt?
This does give a more accurate result, but is still not exact. You can 
get even better results by doing this one step, then backing up and 
using the average of the acceleration at the original position and the 
first approximation of the new position to compute the real new 
position, thus taking into account the fact that the acceleration 
changes over the time step. This isn't as good as you can get either, 
but it's easy to compute.


> I was hoping for a solution where dt= one povray frame.

That is probably a bad idea. The accuracy of the simulation will then 
depend on the number of frames in the animation. You may get very 
different results when you do a longer render.


> So I'm falling back to three ideas:
> 1) it is in fact impossible

Exact solutions are impossible. "Close enough" solutions are quite 
possible.

-- 
Christopher James Huff <cja### [at] earthlinknet>
http://home.earthlink.net/~cjameshuff/
POV-Ray TAG: <chr### [at] tagpovrayorg>
http://tag.povray.org/

Post a reply to this message

Okay cool.

Here's a question related to my ability to detect a "failed model".

Take a "sun".
Take any mass and throw it in any direction not exactly at its center at any
velocity.

If the mass doesn't sink into the sun,  will it:

i)  undertake an "orbit"  "right there" in the sense that it will undertake
an elliptical orbit which intersects that point of release,
OR
ii) will a particle "destined" for a stable orbit sometimes "drift" out
towards the radius of its stable orbit.

I'm ultimately trying to interpret my :
http://news.povray.org/povray.binaries.images/thread/<3A65E429.18F2DC62@my-dejanews.com>/?ttop=184334&toff=4400

If I had  been much, much more patient, might some of those orbits which
looked like they were decaying be just wandering off to something "close
enough" to their own stable orbit?





> > So I'm falling back to three ideas:
> > 1) it is in fact impossible
>
> Exact solutions are impossible. "Close enough" solutions are quite
> possible.
>
> -- 
> Christopher James Huff <cja### [at] earthlinknet>
> http://home.earthlink.net/~cjameshuff/
> POV-Ray TAG: <chr### [at] tagpovrayorg>
> http://tag.povray.org/

Post a reply to this message

Among other things, Christopher James Huff wrote:

> There is no way to compute an exact solution for 3 or more
> gravitationally interacting bodies. However, you can compute a very
> close approximation. I did a simulation (in C++) that was able to
> simulate the orbits of the planets from Mercury to Jupiter for several
> Jovian years (several dozen Earth years), with very little drift. (I
> simply haven't entered the data for the other planets)
> There are some rather complex algorithms needed to do this with both
> accuracy and efficiency, however. And there will always be cumulative
> errors, the simulation will diverge from reality over time, there is no
> way around that.

Exactly. In fact, such systems have a chaotic behaviour. By "chaotic" I 
don't mean they're undeterministic, but they show a strong dependence on 
initial conditions and divergence. If you don't mimic the *exact* initial 
conditions and don't solve *exactly* the movement equations, after some 
(maybe short) time, the solution you get is completely different from the 
real solution. It's the same thing that happens with weather forecast, 
although the equations look simpler.

>> Jelby, in your second post,  you offer an interesting idea.
>> >  p2 = p1 + v1*dt + 1/2*a1*(dt^2)
>> But is fundamentally any different from going to a smaller time slice? 
>> And
> 
> This basically assumes that acceleration will be constant over the time
> interval, and takes that into account for the new point. I take it you
> were originally using p2 = p1 + v1*dt?
> This does give a more accurate result, but is still not exact. You can
> get even better results by doing this one step, then backing up and
> using the average of the acceleration at the original position and the
> first approximation of the new position to compute the real new
> position, thus taking into account the fact that the acceleration
> changes over the time step. This isn't as good as you can get either,
> but it's easy to compute.

Yep, like all others this is only an approximation. It is designed to be 
cheap in computation power and storage. Position is computed using a simple 
second-order approximation. Velocity is computed using a first-order 
approximation, but taking the average of the acceleration. There are other 
algorithms, called "predictor-corrector", which basically do what you 
suggest.

Incidentally, the program I mentioned (Moldy) uses one of these algorithms. 
Even if it's written in C(++?), the algorithm is outlined in its manual, 
it's worth a look, I believe.

>> I was hoping for a solution where dt= one povray frame.
> 
> That is probably a bad idea. The accuracy of the simulation will then
> depend on the number of frames in the animation. You may get very
> different results when you do a longer render.

In the field of molecular dynamics, you usually don't store the data every 
time-step. To get a precise enough solution you need small dt, but then the 
configurations are too similar to be interesting. I'd aim for a frame every 
100*dt, more or less, and see what happens.

-- 
light_source{9+9*x,1}camera{orthographic look_at(1-y)/4angle 30location
9/4-z*4}light_source{-9*z,1}union{box{.9-z.1+x clipped_by{plane{2+y-4*x
0}}}box{z-y-.1.1+z}box{-.1.1+x}box{.1z-.1}pigment{rgb<.8.2,1>}}//Jellby

Post a reply to this message

In article <404178f6$1@news.povray.org>,
 "Greg M. Johnson" <gregj;-)565### [at] aolcom> wrote:

> Take a "sun".
> Take any mass and throw it in any direction not exactly at its center at any
> velocity.
> 
> If the mass doesn't sink into the sun,  will it:
> 
> i)  undertake an "orbit"  "right there" in the sense that it will undertake
> an elliptical orbit which intersects that point of release,

In a 2-mass system, just the sun and the orbiting mass, it will either 
do exactly that or escape the system entirely.


> OR
> ii) will a particle "destined" for a stable orbit sometimes "drift" out
> towards the radius of its stable orbit.

2 body systems are stable. Maybe you're mistakenly thinking of 
elliptical orbits as unstable?


> If I had  been much, much more patient, might some of those orbits which
> looked like they were decaying be just wandering off to something "close
> enough" to their own stable orbit?

If a 2-body system drifts significantly, it is your simulation which is 
unstable. In the real world, none of the planets in our solar system 
follow perfect elliptical orbits, and they do drift around and transfer 
energy to each other, but the level of instability is too small to 
matter. When the sun goes into its red giant stage and swallows up the 
inner planets, they will still be in very similar orbits to what they 
have today.

The p2 = p1 + v1*dt + 1/2*a1*(dt^2) equation will give a more accurate 
simulation than p2 = p1 + v1*dt, but for my simulation of the inner 5 
planets, I found it insufficient. The two-step method I described was 
sufficient, the simulation was stable for at least several decades 
without trying very hard.

-- 
Christopher James Huff <cja### [at] earthlinknet>
http://home.earthlink.net/~cjameshuff/
POV-Ray TAG: <chr### [at] tagpovrayorg>
http://tag.povray.org/

Post a reply to this message

In article <cjameshuff-186E57.00003729022004@news.povray.org>,
 Christopher James Huff <cja### [at] earthlinknet> wrote:

> > So I'm falling back to three ideas:
> > 1) it is in fact impossible
> 
> Exact solutions are impossible. "Close enough" solutions are quite 
> possible.

BTW, here's a demo of several different integrators and orbit types:
http://www.princeton.edu/~rvdb/JAVA/astro/galaxy/Galaxy.html

-- 
Christopher James Huff <cja### [at] earthlinknet>
http://home.earthlink.net/~cjameshuff/
POV-Ray TAG: <chr### [at] tagpovrayorg>
http://tag.povray.org/

Post a reply to this message