If you've been programming in Java, you actually probably have a better intuitive understanding of pointers than you realize, though you don't know the C syntax. But that's neither here nor there, as I've been procrastinating on some work that I have to do, and I ended up translating the C++ implementation from UTexas into Java so I could keep from getting back to it for a little longer...
First, I noticed you asked about 6D phase space. All this means is that you consider your array of "positions" to also include velocities, i.e. for each particle, you consider its "location" to be (x,y,z,vx,vy,vz). That means that "velocities" end up being (vx,vy,vz,ax,ay,az). These integrators work with first order differential equations in general, and it so happens that physical equations are often better solved as first order ODEs in 6 dimensions than second order ODEs in 3 (first order would be dx/dt = v, second order would be d^2x/dt^2 = a - physics usually defines equations using the second order ones, but you can always "unroll" second order ODEs into twice as many first order ones; it so happens that second order solvers are considerably nastier than first order ones, so we tend to do this unrolling when simulating). Below I'll show you what this actually means to you in code, in case that doesn't make sense.
Also, I'm sorry that the MathWorld folks were unhelpful; mathematicians in general can be real jerks (the word I'm really thinking of is a lot more colorful, but I'll keep it PG-13 here...), especially to programmers. But I do suggest that learning the Greek alphabet probably will be worth your time, especially if you're doing physics stuff; there's only so many letters in the regular alphabet, so we tend to run out pretty quickly! And it's hard to follow equations when you can't subvocalize what you're reading, so it helps to know the difference between a zeta and a xi. Eventually when you see them in equations you'll stop panicking and realize that the equations are not really any more difficult than the normal looking ones you're used to seeing.
Anyhow, here is the code, translated from C++ into the form I would probably use if I were using RK4 in Java (I also changed x->t and y->x to conform to my expectation that t=time and x=location). This won't run as is, because there are some things you need to fill in (mainly you need to provide something to set the accelerations of the particles), and possibly because I've made a mistake somewhere (I didn't actually try to get this to run yet). But this should at least give you an idea what's going on in the algorithm.
Code:
float[] x = coordinates of all particles in PHASE SPACE, i.e. including velocities:
x[0] = particle[0].position.x, x[1] = particle[0].position.y,
x[2] = particle[0].velocity.vx, x[3] = particle[0].velocity.vy,
x[4] = particle[1].position.x, etc.
float t = current time
float h = time step
//evaluate derivatives wrt t of each variable in x array, load into dxdt (assumed pre-initialized to correct length)
void rhs_eval(float t, float[] x, float[] dxdt){
if (dxdt.length != x.length) { return; } //add some error-handling code here, because array lengths must be equal
int POSX = 0; //defined merely for clarity below
int POSY = 1;
int VELX = 2;
int VELY = 3;
for (int i=0; i<x.length; i += 4){
dxdt[i+POSX] = x[i+VELX]; //derivative of position is velocity
dxdt[i+POSY] = x[i+VELY];
dxdt[i+VELX] = [x component of the acceleration for particle i/4 based on time and configuration passed (t and x)];
dxdt[i+VELY] = [same for y component];
}
}
//depends on rhs_eval function - Java doesn't have a simple mechanism for passing functions, so we treat it as a global function
void rk4_fixed (float t, float[] x, float h){
int n = x.length;
// Declare local arrays
float[] k1, k2, k3, k4, f, dxdt;
k1 = new float[n];
k2 = new float[n];
k3 = new float[n];
k4 = new float[n];
f = new float[n];
dxdt = new float[n];
// Zeroth intermediate step
rhs_eval(t, x,dxdt);
for (int j = 0; j < n; j++){
k1[j] = h * dxdt[j];
f[j] = x[j] + k1[j] / 2.;
}
// First intermediate step
rhs_eval(t + h / 2., f,dxdt);
for (int j = 0; j < n; j++){
k2[j] = h * dxdt[j];
f[j] = x[j] + k2[j] / 2.;
}
// Second intermediate step
rhs_eval(t + h / 2., f,dxdt);
for (int j = 0; j < n; j++){
k3[j] = h * dxdt[j];
f[j] = x[j] + k3[j];
}
// Third intermediate step
rhs_eval(x + h, f, dxdt);
for (int j = 0; j < n; j++){
k4[j] = h * dxdt[j];
}
// Actual step
for (int j = 0; j < n; j++){
x[j] += k1[j] / 6. + k2[j] / 3. + k3[j] / 3. + k4[j] / 6.;
}
x += h;
return;
}
It's not very OO, but it is what it is. Any questions, feel free to ask.
Also, I've spent more time messing around with Verlet integration than I probably should have, so if you have any questions there, I might be able to help, conceptually at least, if not with actual code; I've implemented many different interactions in that framework, and if you're having trouble with something I very well may have fought through the same thing before.
Good luck!
Pages: 1
Runge Kutta (RK4) for dummies (Read 5043 times)