 |
|
|
|
|
|
| |
| |
|
|
|
|
| |
| |
|
|
...but entirely deterministic! ;-)
OK, so check this craziness out:
http://www.youtube.com/watch?v=duy8s8C7-Uc
http://www.youtube.com/watch?v=58_s6r7PaKo
Each video is approximately 5 hours of computer time. And trust me, when
you watch this full-screen in high resolution, it looks insane! :-D
Pitty YouTube doesn't permit such things. (The file I uploaded was about
35MB, encoded with "H.264". I might put it online somewhere else for you
folks to enjoy...)
So, let me know what you think. :-D
PS. I still need to fix the rendering program. I wrote it while at work,
and it worked OK, but was rather inflexible and brittle. So when I got
home, I wrote it again. But for some reason, the version I wrote at home
is about 10x slower - and I can't figure out why! o_O
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
> http://www.youtube.com/watch?v=duy8s8C7-Uc
> http://www.youtube.com/watch?v=58_s6r7PaKo
Very nice.
> Each video is approximately 5 hours of computer time.
You're just begging me to try and write a GPU version that runs at 100fps
aren't you :-)
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
Invisible <voi### [at] dev null> wrote:
> Each video is approximately 5 hours of computer time.
I wonder how much faster it would be if it was properly implemented in
a *real* programming language... ;)
--
- Warp
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
scott wrote:
>> http://www.youtube.com/watch?v=duy8s8C7-Uc
>> http://www.youtube.com/watch?v=58_s6r7PaKo
>
> Very nice.
Thanks.
YouTube hardly does it justice, but you can see what I mean about the
way it moves, eh?
>> Each video is approximately 5 hours of computer time.
>
> You're just begging me to try and write a GPU version that runs at
> 100fps aren't you :-)
If you think you can do adaptive RK4 numerical integration of a few
thousand particles at realtime speeds, be my freakin' guest! ;-)
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
Warp wrote:
> I wonder how much faster it would be if it was properly implemented in
> a *real* programming language... ;)
Is that a challenge? :-P
Do you know how to numerically integrate a physical simulation? ;-)
You're welcome to try if you like...
You're not totally wrong though: When I initially implemented this, I
discovered, to my astonishment, that converting the state of the system
into an image took vastly longer than all that complicated Runge-Kutta
numerical integration stuff. Which, obviously, is absurd.
On closer investigation, the program was spending 80% of its time in one
function: floor. Looking at the source code, I discovered that floor is
implemented by converting a floating-point double into an
infinite-precision rational and then converting that into an integer.
Like, WTF?
It seems that while *most* numeric operations are very fast, one or two
of the more unusual ones are jaw-droppingly slow.
Anyway, I replaced floor with a direct call to the compiler's low-level
internals library and now converting to an image is near instantaneous.
(Like you'd expect!) It takes longer to suck the data off disk than to
process it.
Now if I could just find a way to make the integration go faster... ;-)
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
http://blog.orphi.me.uk/archives/365
Some more background for any interested souls.
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
>>> Each video is approximately 5 hours of computer time.
>>
>> You're just begging me to try and write a GPU version that runs at 100fps
>> aren't you :-)
>
> If you think you can do adaptive RK4 numerical integration of a few
> thousand particles at realtime speeds, be my freakin' guest! ;-)
I wonder how fast it would be in C++ :-)
Are the three attractors just applying a 1/r^2 force on each particle?
What resolution are you running this at? "a few thousand" sounds like 64x64
or something.
When I get time I will give it a shot, but it might not be a while I'm
preparing today and tomorrow for two weeks of business travel.
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
>> If you think you can do adaptive RK4 numerical integration of a few
>> thousand particles at realtime speeds, be my freakin' guest! ;-)
>
> I wonder how fast it would be in C++ :-)
Go ahead, knock yourself out. ;-)
> Are the three attractors just applying a 1/r^2 force on each particle?
Nope. Currently just 1/r for the "magnets", plus an r force for "gravity".
> What resolution are you running this at? "a few thousand" sounds like
> 64x64 or something.
It's 1,728 x 972 pixels. (And one particle per pixel, obviously.)
I *said* YouTube doesn't do it justice. ;-)
> When I get time I will give it a shot, but it might not be a while I'm
> preparing today and tomorrow for two weeks of business travel.
Yeah, sure, whatever. It's pretty fun stuff to play with though...
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
>> What resolution are you running this at? "a few thousand" sounds like
>> 64x64 or something.
>
> It's 1,728 x 972 pixels. (And one particle per pixel, obviously.)
Ah ok, so "a few" means "over 1600" :-)
> I *said* YouTube doesn't do it justice. ;-)
I can well imagine.
> Yeah, sure, whatever. It's pretty fun stuff to play with though...
Also I can imagine, this is the sort of stuff I like playing about with too.
I wonder if the plane seat will have a power socket for my laptop ... hmmmm.
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
That is brilliant! Well done you. :)
Add this to the end of the URL: &fmt=22
Not only hi-res, but bigger too... ;)
~Steve~
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
>> I wonder how fast it would be in C++ :-)
>
> Go ahead, knock yourself out. ;-)
>
>> Are the three attractors just applying a 1/r^2 force on each particle?
>
> Nope. Currently just 1/r for the "magnets", plus an r force for "gravity".
>
>> What resolution are you running this at? "a few thousand" sounds like
>> 64x64 or something.
>
> It's 1,728 x 972 pixels. (And one particle per pixel, obviously.)
OK so I quickly modified my existing RK4 code to do integration steps on
1728*972 particles in 2D with 4x 1/r attractors. So it takes about 2.5
seconds per integration step, does this compare with your results?
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
scott wrote:
>>> What resolution are you running this at? "a few thousand" sounds
>>> like 64x64 or something.
>>
>> It's 1,728 x 972 pixels. (And one particle per pixel, obviously.)
>
> Ah ok, so "a few" means "over 1600" :-)
1,679,616 pixels, baby! :-D
(Or what you DSLR people would call "1.7 megapixels". :-P )
>> I *said* YouTube doesn't do it justice. ;-)
>
> I can well imagine.
When I get home I'll put up some still frames in full resolution.
> Also I can imagine, this is the sort of stuff I like playing about with
> too. I wonder if the plane seat will have a power socket for my laptop
> ... hmmmm.
Muhuhuhu! >:-D
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
scott wrote:
> OK so I quickly modified my existing RK4 code to do integration steps on
> 1728*972 particles in 2D with 4x 1/r attractors. So it takes about 2.5
> seconds per integration step, does this compare with your results?
It's roughly comparable, certainly.
(What hardware are you running it on?)
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
>> OK so I quickly modified my existing RK4 code to do integration steps on
>> 1728*972 particles in 2D with 4x 1/r attractors. So it takes about 2.5
>> seconds per integration step, does this compare with your results?
>
> It's roughly comparable, certainly.
>
> (What hardware are you running it on?)
That was using 1 core of an Intel E8500 @ 3.16 GHz.
Of course someone has already done particle system simulation on the GPU
before:
http://2ld.de/gdc2004/
Running option 2 of that demo on my GPU gives 15fps for a million particles,
so for 1728x972 I make that about 9fps. If you restrict it to 2D rather
than 3D it should go a little faster.
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
scott wrote:
>>> OK so I quickly modified my existing RK4 code to do integration steps
>>> on 1728*972 particles in 2D with 4x 1/r attractors. So it takes
>>> about 2.5 seconds per integration step, does this compare with your
>>> results?
>>
>> It's roughly comparable, certainly.
>>
>> (What hardware are you running it on?)
>
> That was using 1 core of an Intel E8500 @ 3.16 GHz.
OK. My hardware is considerably less powerful.
Is that RK4 with a single step, or with adaptive step subdivision?
> Of course someone has already done particle system simulation on the GPU
> before:
>
> http://2ld.de/gdc2004/
>
> Running option 2 of that demo on my GPU gives 15fps for a million
> particles, so for 1728x972 I make that about 9fps. If you restrict it
> to 2D rather than 3D it should go a little faster.
Heh, nice... Coming to a gaming station near you, eh?
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
>> That was using 1 core of an Intel E8500 @ 3.16 GHz.
>
> OK. My hardware is considerably less powerful.
>
> Is that RK4 with a single step, or with adaptive step subdivision?
Just one complete integration step of plain RK4. I didn't do any checking
to see if sub-steps are needed or not.
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
scott wrote:
>>> That was using 1 core of an Intel E8500 @ 3.16 GHz.
>>
>> OK. My hardware is considerably less powerful.
>>
>> Is that RK4 with a single step, or with adaptive step subdivision?
>
> Just one complete integration step of plain RK4. I didn't do any
> checking to see if sub-steps are needed or not.
That's probably why my version is slower then.
It computes one step forward. It then recomputes that step as two
half-steps. Finally, it computes the difference between the two answers.
If the difference is too large, the first result is thrown away, and
further subdivision is done to check the second result is OK. (There is
a maxmimum number of subdivisions though - or rather, a minimum step size.)
I can try turning that off to see what speed I get - but I'm still
running on different hardware to you, so it's probably not very meaningful.
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
> It computes one step forward. It then recomputes that step as two
> half-steps.
That sounds like a waste of time. Within the RK4 algorithm you already have
4 estimates of slope, just check the variance of these to determine whether
you need to use a smaller time step or not.
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
scott wrote:
>> It computes one step forward. It then recomputes that step as two
>> half-steps.
>
> That sounds like a waste of time. Within the RK4 algorithm you already
> have 4 estimates of slope, just check the variance of these to determine
> whether you need to use a smaller time step or not.
Well, the stability depends on how sharply the force acting on the
particle changes as a function of position...
I have to say, I don't really comprehend how RK4 works. I copied the
formulas, and it seems to produce the right answers, but I have no idea
*why* it works or what it's doing.
Still, I could try what you suggest and see what the accuracy /
performance is like. (I already had to decrease the error tollerance
once to get crisp images...)
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
> I have to say, I don't really comprehend how RK4 works. I copied the
> formulas, and it seems to produce the right answers, but I have no idea
> *why* it works or what it's doing.
It's easier to understand if you consider a single function to integrate
first, rather than doing acceleration/velocity and velocity/position in one
go. There is a clever trick that lets you do both of these together without
much extra work, but it's not too obvious and usually confuses things.
Let's suppose you have a function that calculates the instantaneous velocity
of an object based on its position (x) and a time (t):
v = F(x,t)
Now using Euler integration you can easily track the position of the object
like this:
> t = 0
> x = 0
> repeat
> average_velocity = F(x,t)
> x = x + average_velocity * dt
> t = t + dt
Where dt is some time step.
As you are aware, the problem is that with large time-steps big errors
appear because you are assuming the average velocity during the time step is
equal to the instantaneous velocity at the beginning of the time step - this
is only approximately correct!
All RK4 does is to find a much better estimate of the average_velocity
during the time step, then you can use that to advance x with much higher
accuracy:
> t = 0
> x = 0
> repeat
> average_velocity = RK4_Algorithm( x , t , dt )
> x = x + average_velocity * dt
> t = t + dt
The RK4 algorithm looks like this:
> Function RK4_Algorithm( x , t , dt )
> AVE1 = F(x,t)
First step there is to compute an Average Velocity Estimate by simply using
F(x,t) like we did in Euler.
> AVE2 = F( x + AVE1 * dt/2 , t + dt/2)
This is the first clever bit, we are going to calculate another estimate for
the velocity *half* a time step forward. We use the first estimate of
velocity (AVE1) to work out what the position will be half a time step
later, then we then use this new position in our velocity finding function
to find the velocity there. This gives us our estimate of velocity at the
middle of the time step, call it AVE2.
> AVE3 = F( x + AVE2 * dt/2 , t + dt/2)
This does the same as before, but instead of using AVE1 as the assumption
for average velocity, we use AVE2.
> AVE4 = F( x + AVE3 , t + dt )
This is the 4th estimate, it uses the 3rd one to do a *full* time step
forward, and calculate the velocity there (at the *end* of the time step).
Now we have four estimates for the velocity, one at the beginning of the
time step (standard Euler), two in the middle, and one at the end of the
time step.
The RK4 algorithm then tells us to take a weighted average of these
velocities:
> AVE = (AVE1 + 2*AVE2 + 2*AVE3 + AVE4) / 6
That is it! That is the estimate of average velocity over the entire time
step. You will have noticed we called F four times, but doing this is MORE
accurate than simply doing four dt/4 Euler time steps.
> return AVE
You should be able to simply compare the values of AVE1,AVE2,AVE3,AVE4 to
determine if you need to go back and do a smaller time step.
Now to make things more general, you can replace the position "x" by the
generic "state" of the system, and the velocities by the "derivative of
state" of the system. This way, the functions F and RK4_Algorithm will
return objects of type "derivative of state".
If you are dealing with forces acting on particles, then the "state" of each
particle is its position and momentum (4 variables in 2D), and the
"derivative of state" is its velocity and force (another 4 variables).
Now, your function F must return the velocity and force. The force it
calculates based on your magnet formula or whatever, and the velocity it can
easily calculate by dividing the supplied momentum value by the mass. That
last part is the key that links the two differential equations together, and
why it is sometimes complicated how RK4 can do both at the same time for
acceleration and velocity.
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
Invisible wrote:
> http://www.youtube.com/watch?v=duy8s8C7-Uc
> http://www.youtube.com/watch?v=58_s6r7PaKo
A giant fly staring at me through my monitor freaks me out!
--
"Smoking helps you lose weight -- one lung at a time!"
/\ /\ /\ /
/ \/ \ u e e n / \/ a w a z
>>>>>>mue### [at] nawazorg<<<<<<
anl
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
Mueen Nawaz wrote:
> A giant fly staring at me through my monitor freaks me out!
More or less so than a squashed bug surrounded by smaller squashed bugs? ;-)
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
Invisible wrote:
> 1,679,616 pixels, baby! :-D
>
> (Or what you DSLR people would call "1.7 megapixels". :-P )
Actually that would be five megapixels. Almost all the digital cameras count
each cell of the CCD as a pixel, so a RGB triad takes three pixels.
Some of the cameras you see where they have a full three colors per pixel
triple the megapixel count in their marketing.
--
Darren New, San Diego CA, USA (PST)
There's no CD like OCD, there's no CD I knoooow!
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
Darren New wrote:
> Invisible wrote:
>> 1,679,616 pixels, baby! :-D
>>
>> (Or what you DSLR people would call "1.7 megapixels". :-P )
>
> Actually that would be five megapixels. Almost all the digital cameras
> count each cell of the CCD as a pixel, so a RGB triad takes three pixels.
>
> Some of the cameras you see where they have a full three colors per
> pixel triple the megapixel count in their marketing.
Epic fail. O_O
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
scott wrote:
> It's easier to understand if you consider a single function to integrate
> first, rather than doing acceleration/velocity and velocity/position in
> one go.
Er, yes.
OTOH, almost all physical systems are of this type.
(Unfortunately, all the texts explaining this stuff use abstract
mathematical functions, which makes it extremely difficult to figure out
what's going on.)
> There is a clever trick that lets you do both of these together
> without much extra work, but it's not too obvious and usually confuses
> things.
I'll bet...
> As you are aware, the problem is that with large time-steps big errors
> appear because you are assuming the average velocity during the time
> step is equal to the instantaneous velocity at the beginning of the time
> step - this is only approximately correct!
>
> All RK4 does is to find a much better estimate of the average_velocity
> during the time step, then you can use that to advance x with much
> higher accuracy.
So it gets a more accurate average velocity by sampling it more?
> Now we have four estimates for the velocity, one at the beginning of the
> time step (standard Euler), two in the middle, and one at the end of the
> time step.
Right. I don't really understand why it has to be these exact points...
> The RK4 algorithm then tells us to take a weighted average of these
> velocities:
...or why this specific weighting is the correct one. But hey, I guess
somebody had a reason for that.
> That is it! That is the estimate of average velocity over the entire
> time step. You will have noticed we called F four times, but doing this
> is MORE accurate than simply doing four dt/4 Euler time steps.
Clearly I'm going to have to meditate on exactly why that's the case. I
would have thought that taking 4 seperate steps would be more accurate
than computing 4 steps and then blurring them all into one step, but hey.
> You should be able to simply compare the values of AVE1,AVE2,AVE3,AVE4
> to determine if you need to go back and do a smaller time step.
Mmm, OK. I'll take a look at some numerical work and see what happens...
> Now to make things more general, you can replace the position "x" by the
> generic "state" of the system, and the velocities by the "derivative of
> state" of the system. This way, the functions F and RK4_Algorithm will
> return objects of type "derivative of state".
>
> If you are dealing with forces acting on particles, then the "state" of
> each particle is its position and momentum (4 variables in 2D), and the
> "derivative of state" is its velocity and force (another 4 variables).
>
> Now, your function F must return the velocity and force. The force it
> calculates based on your magnet formula or whatever, and the velocity it
> can easily calculate by dividing the supplied momentum value by the
> mass. That last part is the key that links the two differential
> equations together, and why it is sometimes complicated how RK4 can do
> both at the same time for acceleration and velocity.
Hmm... How is momentum and velocity different?
Oh, wait - the mass could be non-unital. I forgot about that...
Thinking about Euler integration, you end up doing something like
x' = x' + dt x''
x = x + dt x'
In other words, you compute the new velocity from the old one, and then
assume that velocity is constant while you compute the new position from
the old one. But if you assume that acceleration is constant, you can
actually *exactly solve* for position. I can't figure out the equation
off the top of my head, but it would be some kind of quadratic curve. I
wonder - would this give you any more of an accurate result?
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
Invisible wrote:
> Darren New wrote:
>> Invisible wrote:
>>> 1,679,616 pixels, baby! :-D
>>>
>>> (Or what you DSLR people would call "1.7 megapixels". :-P )
>>
>> Actually that would be five megapixels. Almost all the digital cameras
>> count each cell of the CCD as a pixel, so a RGB triad takes three pixels.
>>
>> Some of the cameras you see where they have a full three colors per
>> pixel triple the megapixel count in their marketing.
>
> Epic fail. O_O
Not really. Pixel means "picture element". Who's to say whether that's one
CCD cell or one tri-color cell? :-)
--
Darren New, San Diego CA, USA (PST)
There's no CD like OCD, there's no CD I knoooow!
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
Darren New <dne### [at] sanrrcom> wrote:
> Not really. Pixel means "picture element". Who's to say whether that's one
> CCD cell or one tri-color cell? :-)
I'm wondering: Given that the blue component is significantly less
important than the other components (as demonstrated in another thread),
wouldn't it be cheaper to produce CCD cells where there are only half of
the blue cells? In other words, each second RGB triplet would actually lack
the B detector (the blue component can then be interpolated from adjacent
pixels in software). This would most probably not produce pictures of
significantly less quality.
--
- Warp
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
On Thu, 07 May 2009 20:13:08 +0200, Warp <war### [at] tagpovrayorg> wrote:
> I'm wondering: Given that the blue component is significantly less
> important than the other components (as demonstrated in another thread),
> wouldn't it be cheaper to produce CCD cells where there are only half of
> the blue cells? In other words, each second RGB triplet would actually
> lack the B detector
The sensor itself (be it CCD or CMOS) is monochrome in the vast majority
of digital cameras, the notable exceptions being the handful of cameras
based on Foveon sensors.
Most digital cameras use a Bayer filter to get colour; these have half the
cells dedicated to green and a quarter each to red and blue.
> the blue component can then be interpolated from
> adjacent pixels in software
This already happens with all three colour components.
--
FE
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
> OTOH, almost all physical systems are of this type.
Indeed, and the exact solution of that is a damped sine wave, which is why
so much stuff vibrates and resonates in the real world.
>> Now we have four estimates for the velocity, one at the beginning of the
>> time step (standard Euler), two in the middle, and one at the end of the
>> time step.
>
> Right. I don't really understand why it has to be these exact points...
>
>> The RK4 algorithm then tells us to take a weighted average of these
>> velocities:
>
> ...or why this specific weighting is the correct one. But hey, I guess
> somebody had a reason for that.
It matches the Taylor Series approximation of the function:
http://grb.physics.unlv.edu/~zbb/files/RungeKuttaProof.pdf
> Clearly I'm going to have to meditate on exactly why that's the case. I
> would have thought that taking 4 seperate steps would be more accurate
> than computing 4 steps and then blurring them all into one step, but hey.
RK4 is not just taking 4 blind steps though, it's using the result of the
previous ones to get the next ones, which is how it gets to be more accurate
than a simple linear estimate.
>> Now, your function F must return the velocity and force. The force it
>> calculates based on your magnet formula or whatever, and the velocity it
>> can easily calculate by dividing the supplied momentum value by the mass.
>> That last part is the key that links the two differential equations
>> together, and why it is sometimes complicated how RK4 can do both at the
>> same time for acceleration and velocity.
>
> Hmm... How is momentum and velocity different?
>
> Oh, wait - the mass could be non-unital. I forgot about that...
Yes that's rather key :-) I also find that keeping momentum as a state
variable rather than velocity (which would also work) avoids some of the
confusion with the differential of position in the algorithm.
> Thinking about Euler integration, you end up doing something like
>
> x' = x' + dt x''
> x = x + dt x'
>
> In other words, you compute the new velocity from the old one, and then
> assume that velocity is constant while you compute the new position from
> the old one. But if you assume that acceleration is constant, you can
> actually *exactly solve* for position. I can't figure out the equation off
> the top of my head,
x = x + dt x' + 1/2 x'' dt^2
That looks a bit like a Taylor Series expansion, doesn't it :-) You can
extend that to:
x = x + dt x' + 1/2 x'' dt^2 + 1/6 x''' dt^3
and also write it for x' instead of x:
x' = x' + dt x'' + 1/2 x''' dt^2 + 1/6 x'''' dt^3
And that is exactly what your RK4 algorithm is estimating :-)
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
scott wrote:
>> OTOH, almost all physical systems are of this type.
>
> Indeed, and the exact solution of that is a damped sine wave, which is
> why so much stuff vibrates and resonates in the real world.
Well, I don't know about that, but physics is mostly about forces, and
forces generate acceleration, and force is usually a function of
position (or perhaps velocity). Hence differential equations seem to be
kind of ubiquitous in physics.
(My memory may be defective, but I was under the impression that indeed
differential and integral calculus were *invented* by the eniment
physacist Newton...)
Interestingly, a trivial differential equation has exponential growth as
its solution, but you never see this in mechanical systems.
>> Right. I don't really understand why it has to be these exact points...
>> ...or why this specific weighting is the correct one. But hey, I guess
>> somebody had a reason for that.
>
> It matches the Taylor Series approximation of the function:
Taylor...series...?
I thought a Taylor series is just a way of constructing a polynomial
approximation to a function that you don't have any better way to
compute. (Assuming the series is even convergent in the first place.)
I'm not sure how you use that to construct an approximation to a
completely *unknown* function...
> http://grb.physics.unlv.edu/~zbb/files/RungeKuttaProof.pdf
"Hence the idea in Theorem 9.7 is that if the step size in the RK4
method is reduced by a factor of 1/2 we can expect that the overall
F.G.E. will be reduced by a factor of 1/16."
...OK, well *that* would certainly explain why it's drastically more
stable then. (Although it still doesn't explain how this phenominon
comes about in the first place.)
>> Clearly I'm going to have to meditate on exactly why that's the case.
>> I would have thought that taking 4 seperate steps would be more
>> accurate than computing 4 steps and then blurring them all into one
>> step, but hey.
>
> RK4 is not just taking 4 blind steps though, it's using the result of
> the previous ones to get the next ones, which is how it gets to be more
> accurate than a simple linear estimate.
And taking 4 linear steps uses the result of the previous steps for the
next one. I don't get why taking those steps, throwing the result away
and just keeping the intermediate results, and then mixing that all
together gives a better result. (Altough it clearly *does*.)
>> Hmm... How is momentum and velocity different?
>>
>> Oh, wait - the mass could be non-unital. I forgot about that...
>
> Yes that's rather key :-)
Only if the particle has non-unital mass. ;-)
> I also find that keeping momentum as a state
> variable rather than velocity (which would also work) avoids some of the
> confusion with the differential of position in the algorithm.
Perhaps.
>> But if you assume that acceleration is
>> constant, you can actually *exactly solve* for position. I can't
>> figure out the equation off the top of my head,
>
> x = x + dt x' + 1/2 x'' dt^2
Yeah, something like that. :-}
> That looks a bit like a Taylor Series expansion, doesn't it :-)
I wouldn't know.
> You can extend that to:
>
> x = x + dt x' + 1/2 x'' dt^2 + 1/6 x''' dt^3
>
> and also write it for x' instead of x:
>
> x' = x' + dt x'' + 1/2 x''' dt^2 + 1/6 x'''' dt^3
>
> And that is exactly what your RK4 algorithm is estimating :-)
I don't follow. (I.e., I don't see how the final formula is remotely
related to RK4.)
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
> Well, I don't know about that, but physics is mostly about forces, and
> forces generate acceleration, and force is usually a function of position
> (or perhaps velocity). Hence differential equations seem to be kind of
> ubiquitous in physics.
>
> Interestingly, a trivial differential equation has exponential growth as
> its solution, but you never see this in mechanical systems.
Diesel engines would destroy themselves without a control system to keep the
speed down, that's the best I can come up with now of a mechanical system
with exponential growth. I think the key point is that most things have a
force that acts in the opposite direction of the movement, which usually
keeps things stable and also produces the oscillations.
> I thought a Taylor series is just a way of constructing a polynomial
> approximation to a function that you don't have any better way to compute.
Indeed, so if we can estimate the derivatives by some numerical method, we
should be able to estimate how the actual function goes forwards, at least
for a little while.
> And taking 4 linear steps uses the result of the previous steps for the
> next one.
But it doesn't use the results to go back and work on the *same* part of the
function to improve the accuracy, it just blindly moves onto the next time
step ignoring whatever errors might have been present.
> I don't get why taking those steps, throwing the result away and just
> keeping the intermediate results, and then mixing that all together gives
> a better result. (Altough it clearly *does*.)
Another way to look at it I guess is that RK4 is effectively estimating
"deeper" derivatives, like x'' and x''' to better estimate how the function
goes.
>> I also find that keeping momentum as a state variable rather than
>> velocity (which would also work) avoids some of the confusion with the
>> differential of position in the algorithm.
>
> Perhaps.
It's also crucial when you come do doing rotations as well as linear
movements. For linear movements you can simply add on the accerlation times
dt to the velocity, then the same for position. But you cannot do this for
angular movements, you need to add on the torque to angluar momentum, use
that to calculate a thing called spin, then add that onto the angular
orientation.
> I don't follow. (I.e., I don't see how the final formula is remotely
> related to RK4.)
Try writing a function to estimate x' x'' and x''' of a function at a
certain time and position. Then use it to predict what the value of x will
be at the next time step based on the Taylor Series expansion. Then,
assuming you are estimating the derivatives in the "correct" way, you can
simplify your solution and you'll end up with the RK4 algorithm :-) I have
no idea what mathematical jiggery pokery you'll need to do that, but that is
essentially what RK4 is doing.
Euler integration estimates the next value by assuming:
x = x + dt x'
ie it will perfectly match a function like x = 4 + 5t
RK4 estimates the next value by assuming:
x = x + dt x' + 1/2 x'' dt^2 + 1/6 x''' dt^3
ie it will perfectly match a function like x = 4 + 5t - 3t^2 + 45t^4
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
> I think the key point is that most
> things have a force that acts in the opposite direction of the movement,
> which usually keeps things stable and also produces the oscillations.
Yeah.
Of course, the fun part happens when objects hit each other. Then you
have a sudden and drastic change in the forces applied. Good luck
integrating that...
>> I thought a Taylor series is just a way of constructing a polynomial
>> approximation to a function that you don't have any better way to
>> compute.
>
> Indeed, so if we can estimate the derivatives by some numerical method,
> we should be able to estimate how the actual function goes forwards, at
> least for a little while.
So... well, we don't need to "estimate" the acceleration at all. We can
compute it directly. (Assuming we have an exact position...) But you're
saying you compute the acceleration several different times in order to
estimate its integral or something?
> Another way to look at it I guess is that RK4 is effectively estimating
> "deeper" derivatives, like x'' and x''' to better estimate how the
> function goes.
...so you compute acceleration several times to try to estimate the
derivative of acceleration? (I think the term is "jerk".) That makes
slightly more sense.
> It's also crucial when you come do doing rotations as well as linear
> movements. For linear movements you can simply add on the accerlation
> times dt to the velocity, then the same for position. But you cannot do
> this for angular movements, you need to add on the torque to angluar
> momentum, use that to calculate a thing called spin, then add that onto
> the angular orientation.
o_O
Let us pray I never need to handle anything other than point masses...
>> I don't follow. (I.e., I don't see how the final formula is remotely
>> related to RK4.)
>
> Try writing a function to estimate x' x'' and x''' of a function at a
> certain time and position.
Well, if you have a way to compute f(x), you could try using
f'(x) ~= f(x + dx) - f(x)
By analogy you have
f''(x) ~= f'(x + dx) - f'(x)
and substituting the previous you eventually arrive at
f''(x) ~= f(x + 2 dx) - 2 f(x + dx) + f(x)
By a similar process you have
f'''(x) ~= f(x + 3 dx) - 3 f(x + 2 dx) + 3 f(x + dx) - f(x)
Of course, none of this does you any good if you can't compute f in the
first place.
> Then use it to predict what the value of x
> will be at the next time step based on the Taylor Series expansion.
Ouch. That sounds like fun...
> Then, assuming you are estimating the derivatives in the "correct" way,
> you can simplify your solution and you'll end up with the RK4 algorithm
> :-) I have no idea what mathematical jiggery pokery you'll need to do
> that, but that is essentially what RK4 is doing.
Mmm, right.
Of course, now I'm trying to comprehend how this *explicit* method
differs from the more accurate *implicit* method I keep hearing about...
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
> Of course, the fun part happens when objects hit each other. Then you have
> a sudden and drastic change in the forces applied. Good luck integrating
> that...
Usually in real-time simulations you can just apply an impulse to the body
and bypass the integrator altogether. Your collision response algorithm
will calculate exactly what impulse to apply to each body so that they
"bounce" realistically, you then directly modify the momentum values, and in
the next step RK4 doesn't know what happened :-)
A more accurate way is to actually simulate the contact deformation with
forces, but usually due to the big forces involved you need very small time
steps which is beyond real time simulation.
> So... well, we don't need to "estimate" the acceleration at all. We can
> compute it directly. (Assuming we have an exact position...) But you're
> saying you compute the acceleration several different times in order to
> estimate its integral or something?
By computing x' several times (which you can exactly) you are in effect
estimating x'' and x''', which gives you a better estimate for x than just
using x' alone. Remember that here "x" can be position or velocity as far
as RK4 is concerned.
> ...so you compute acceleration several times to try to estimate the
> derivative of acceleration? (I think the term is "jerk".) That makes
> slightly more sense.
Yes.
> Let us pray I never need to handle anything other than point masses...
If you have set up your RK4 code nice and object orientatedly, it is really
very simple to add in.
> Well, if you have a way to compute f(x), you could try using
You don't though, you are trying to find f(x+dt) given that you can easily
calculate f'(x).
f''(x) ~= (f'(x+dt) - f'(x)) / dt
f'''(x) ~= (f''(x+dt) - f''(x)) / dt
f''''(x) ~= (f'''(x+dt) - f'''(x)) / dt
If you substitute those expressions into the Taylor series expansion and
write it in terms of just f'(x) then you will probably end up with something
similar to RK4. Maybe :-)
f(x) = f(x) + f'(x)*dt + 1/2*f''(x)*dt^2 + ...
> Of course, now I'm trying to comprehend how this *explicit* method differs
> from the more accurate *implicit* method I keep hearing about...
You can usually rewrite your equations of motion for your whole simulation
as one big matrix equation, something like Ax + Bx' + Cx'' = D where A-D are
matrices and x is a vector of all your state variables. C would typically
contain masses and inertias, A stiffnesses, B damping and D external forces.
An implicit integrator works on this entire formula in one stab, solving
usually by matrix inversion, iterative methods or whatever. For big
interacting systems it is much better as the constraints (eg two masses
connected by a spring) are forced throughout the calculation, and not just
adjusted at the end of each time step. This is how mechanical finite
element simulation software works, essentially solving a massive matrix
equation for each time step.
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
Invisible wrote:
> PS. I still need to fix the rendering program. I wrote it while at work,
> and it worked OK, but was rather inflexible and brittle. So when I got
> home, I wrote it again. But for some reason, the version I wrote at home
> is about 10x slower - and I can't figure out why! o_O
I just spent about 2 hours peering over the source code for the two
programs, back to back, side by side. Considering that they were written
on different ways on different computers, it's remarkable how many of
the files are identical character-by-character. (Most of them differ by
variable names, function names, order of arguments, etc.)
I even tried replacing source files from the defective version with
files from the optimal version, but to no avail. The program is big, but
not *that* damned big! There's only so many places to look for bugs...
I tried everything I could think of to fix the problem. *Finally*, I
discovered that replacing one block of code with an apparently
equivilent block causes a 10x speedup. Whiskey Tango Foxtrot?
At least, I know the bug is inside one tiny function. But where?
Well, I just nailed it. It was in the adaptive subsampling code. When
computing the error in the current step, I wrote
error = mag (position1 - position3) + mag (velocity1 - position3)
Looks fine, right?
...until you put it on consecutive lines:
errorP = mag (position1 - position3)
errorV = mag (velocity1 - position3)
PWN3D!
Change the second line to
errorV = mag (velocity1 - velocity3)
and suddenly the program speeds up massively.
D'OH!! >_<
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
And now for some speed tuning.
A simple way to measure the program's speed is to count how many frames
per minute it can render. So, set it to 320x480 so it's moderately fast
[on this ancient single-core 32-bit 1.5 GHz CPU] and count.
After fixing the stupid bug: ~24 frames/minute.
Changed array indexing from tuples to integers: ~24 frames/minute.
Merged the four arrays into one interleaved array: ~26 frames/minute.
[I'm thinking this might improve cache behaviour slightly... but
apparently not much. Or just the cache isn't the bottleneck.]
Removed array bounds checks: ~27 frames/minute.
Removed adaptive time step: ~80 frames/minute.
[Well there's a surprise. Pitty the frames look awful...]
Added adaptive time step within RK4 as per Scott's suggestion: ~1
frame/minute.
[Ah, but now the "error" is being calculated in a totally different way,
so the error tolerance required to get a good image may have changed.
What happens if I use a less strict tolerance?]
Turned down error bound: ~68 frames/minute.
[I made the maximum permissible error 10,000x larger, and the images
still come out looking flawless. But they render vastly faster...!]
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
scott wrote:
> I think the key point is that most
> things have a force that acts in the opposite direction of the movement,
I've always found this to be quite fascinating. Why should that be?
Why is back EMF always in the direction that reduces the EMF? Why does
friction work in a direction that reduces the effects of friction?
Almost everything one can think of works this way. I only remember
encountering two or three effects in the last few decades (other than
economics and some odd quantum stuff, IIRC) that don't damp themselves.
Unless you get something like thermal runaway, where you have two separate
effects reinforcing each other that would individually damp each other out.
Funkily enough, if you look at wikipedia under "Negative_mass", there's some
stuff in there about what happens if mass is negative, with fun observations
like "negative masses would produce an attractive force on one another, but
would be repelled because of their negative inertial masses". I.e., the two
minus signs on the mass cancel out, so two negative masses are attracted to
each other, *but* they'll move apart, because negative mass moves towards
the direction you push it from. :-)
</ramble>
--
Darren New, San Diego CA, USA (PST)
There's no CD like OCD, there's no CD I knoooow!
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
Invisible wrote:
> Change the second line to
> errorV = mag (velocity1 - velocity3)
> and suddenly the program speeds up massively.
Shouldn't velocity and position be incompatible types?
--
Darren New, San Diego CA, USA (PST)
There's no CD like OCD, there's no CD I knoooow!
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
Darren New wrote:
> Invisible wrote:
>> Change the second line to
>> errorV = mag (velocity1 - velocity3)
>> and suddenly the program speeds up massively.
>
> Shouldn't velocity and position be incompatible types?
Well, if you wanted to go that far you could. I just had them both as 2D
vectors...
--
http://blog.orphi.me.uk/
http://www.zazzle.com/MathematicalOrchid*
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
Orchid XP v8 wrote:
> Well, if you wanted to go that far you could.
Well, now you know why people do, when it's important to be correct. :-)
Isn't it pretty trivial to do that in Haskell, tho?
--
Darren New, San Diego CA, USA (PST)
There's no CD like OCD, there's no CD I knoooow!
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
Invisible wrote:
> And now for some speed tuning.
>
> A simple way to measure the program's speed is to count how many frames
> per minute it can render. So, set it to 320x480 so it's moderately fast
> [on this ancient single-core 32-bit 1.5 GHz CPU] and count.
>
> After fixing the stupid bug: ~24 frames/minute.
>
> Changed array indexing from tuples to integers: ~24 frames/minute.
>
> Merged the four arrays into one interleaved array: ~26 frames/minute.
>
> [I'm thinking this might improve cache behaviour slightly... but
> apparently not much. Or just the cache isn't the bottleneck.]
>
> Removed array bounds checks: ~27 frames/minute.
>
> Removed adaptive time step: ~80 frames/minute.
>
> [Well there's a surprise. Pitty the frames look awful...]
>
> Added adaptive time step within RK4 as per Scott's suggestion: ~1
> frame/minute.
>
> [Ah, but now the "error" is being calculated in a totally different way,
> so the error tolerance required to get a good image may have changed.
> What happens if I use a less strict tolerance?]
>
> Turned down error bound: ~68 frames/minute.
>
> [I made the maximum permissible error 10,000x larger, and the images
> still come out looking flawless. But they render vastly faster...!]
A clearer (for a human) way to see the speed is to show the inverse of that:
seconds per frame.
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
Darren New wrote:
> Orchid XP v8 wrote:
>> Well, if you wanted to go that far you could.
>
> Well, now you know why people do, when it's important to be correct. :-)
>
> Isn't it pretty trivial to do that in Haskell, tho?
Yeah, something along the lines of
newtype Position = Position Vector2 deriving (Eq, Show, Vector, Num)
newtype Velocity = Velocity Vector2 deriving (Eq, Show, Vector, Num)
Now just replace all occurrances of Vector2 with Position or Velocity as
appropriate.
The only problem is that *now* you can't say
p1 = p0 * dt v0
without sprinkling it with typecasts.
Actually, what you'd *probably* do is define
step :: TimeStep -> Velocity -> Position
step dt (Velocity v) = Position (dt * v)
and then have
p1 = p0 + step dt v0
Only problem now is that you'd have to define another such function for
acceleration - assuming you're going to have a seperate type for that.
There's probably some better way of encoding all this. (Possibly
assosiated types, or maybe just phantom types...)
--
http://blog.orphi.me.uk/
http://www.zazzle.com/MathematicalOrchid*
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
Nicolas Alvarez wrote:
> A clearer (for a human) way to see the speed is to show the inverse of that:
> seconds per frame.
Well, I can run the program and count the number of frames rendered per
minute (approximately). The number of seconds per frame isn't something
I can directly measure. (Of course, you can compute it yourself from my
numbers quite easily...)
--
http://blog.orphi.me.uk/
http://www.zazzle.com/MathematicalOrchid*
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
Darren New wrote:
> scott wrote:
>> I think the key point is that most things have a force that acts in
>> the opposite direction of the movement,
>
> I've always found this to be quite fascinating. Why should that be?
It seems like conservation of energy and the second law of
thermodynamics would essentially require this (barring situations where
there's some external source of energy).
If you had an effect that reinforced itself without using some external
source of energy then any measurable consequence of this effect would
have to increase as a system approached a high-entropy state. Since
there's generally not much interesting to measure in a maximal entropy
state it shouldn't be a surprise that there's dampening forces everywhere.
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
Orchid XP v8 wrote:
> Nicolas Alvarez wrote:
>
>> A clearer (for a human) way to see the speed is to show the inverse of
>> that: seconds per frame.
>
> Well, I can run the program and count the number of frames rendered per
> minute (approximately). The number of seconds per frame isn't something
> I can directly measure. (Of course, you can compute it yourself from my
> numbers quite easily...)
Imagine a 3D renderer. You put one sphere: 500fps. You add a second sphere:
250fps. What? It dropped a *lot* by just adding one more sphere! How is
this going to scale to a hundred objects?
Now, you put one sphere, it takes 2 milliseconds. You add another, 4
milliseconds. That sounds better.
Even though it's the *same numbers*.
Time per frame is clearer than frames per second because the brain can't get
around the inverse curve so easily.
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
> The only problem is that *now* you can't say
>
> p1 = p0 * dt v0
>
> without sprinkling it with typecasts.
>
> Actually, what you'd *probably* do is define
>
> step :: TimeStep -> Velocity -> Position
> step dt (Velocity v) = Position (dt * v)
>
> and then have
>
> p1 = p0 + step dt v0
Couldn't you just define a * function that takes a velocity and a time and
returns a position? You can do the same for one that takes an acceleration
and a time and returns velocity. Then you wouldn't need to change the code,
and * looks better than "step".
In fact, maybe it's possible to do all this automatically if you keep track
of the units somehow? Could you do this in Haskell? Somehow tell it that
distance is units L, velocity is units L*T, etc, and then it would work it
all out automatically when you did multiplication and division. Would be
pretty cool for phyical simulations.
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
> Of course, now I'm trying to comprehend how this *explicit* method
> differs from the more accurate *implicit* method I keep hearing about...
Useful link:
http://www.pixar.com/companyinfo/research/pbm2001/
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
scott wrote:
> In fact, maybe it's possible to do all this automatically if you keep
> track of the units somehow?
Yes.
> Could you do this in Haskell?
That was kind of what I'd meant to ask. How easy is it? Is Haskell high
level enough to make a library for units that works like physical units?
--
Darren New, San Diego CA, USA (PST)
There's no CD like OCD, there's no CD I knoooow!
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
scott wrote:
> Couldn't you just define a * function that takes a velocity and a time
> and returns a position?
No.
I rephrase: Not if you want those to be different types.
The functions +, -, * and / are only defined to operate on values of the
exact same type (and return a result of the same type also).
You can, however, define new functions with the appropriate types. (Say,
have functions named |+|, |-|, etc. or something like that.)
Alternatively, you can provide new definitions for +, -, etc., but then
you have to manually remove access to the standard versions every time
you want to use them.
--
http://blog.orphi.me.uk/
http://www.zazzle.com/MathematicalOrchid*
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
Darren New wrote:
> That was kind of what I'd meant to ask. How easy is it? Is Haskell high
> level enough to make a library for units that works like physical units?
I know that several such libraries exist. I've never tried to use one,
so I couldn't say how they work.
--
http://blog.orphi.me.uk/
http://www.zazzle.com/MathematicalOrchid*
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
scott wrote:
>> Of course, now I'm trying to comprehend how this *explicit* method
>> differs from the more accurate *implicit* method I keep hearing about...
>
> Useful link:
>
> http://www.pixar.com/companyinfo/research/pbm2001/
...those guys! o_O
--
http://blog.orphi.me.uk/
http://www.zazzle.com/MathematicalOrchid*
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
