gravity
Graham Dennis
Example system of gravitationally-attracted particles
t
q
p_1 p_2 p_3 p_4
x y vx vy
initialMomentaNoises
x = 2.0*p_1; y = 2.0*p_2; vx = 0.1*p_3; vy = 0.1*p_4;
m
initialMassNoises
// Set the mass equal to a poissonian noise. // but don't let the mass be zero m = q; if (m < 0.01) m = 1.0;
1000 100
motion
mass
dx_dt = vx; dy_dt = vy; for (long k = 0; k < j; k++) { real inverseSeparationCubed = pow((x(j => k) - x(j => j))*(x(j => k) - x(j => j)) + (y(j => k) - y(j => j))*(y(j => k) - y(j => j)), -3.0/2.0); // printf("j: %li k: %li x[j]: %e x[k]: %e y[j]: %e y[k]: %e\n", j, k, x[j], x[k], y[j], y[k]); // printf("separationSquared: %e inverseSeparationCubed: %e\n", (x[k] - x[j])*(x[k] - x[j]) + (y[k] - y[j])*(y[k] - y[j]), inverseSeparationCubed); dvx_dt(j => j) += m(j => k)*(x(j => k) - x(j => j))*inverseSeparationCubed; dvx_dt(j => k) += m(j => j)*(x(j => j) - x(j => k))*inverseSeparationCubed; dvy_dt(j => j) += m(j => k)*(y(j => k) - y(j => j))*inverseSeparationCubed; dvy_dt(j => k) += m(j => j)*(y(j => j) - y(j => k))*inverseSeparationCubed; }
xR yR vxR vyR
motion
xR = x; yR = y; vxR = vx; vyR = vy;
energy px py
mass motion
// Check conserved quantities energy = 0.5*m*(vx*vx + vy*vy); for (long k = 0; k < j; k++) { energy += -m(j => j)*m(j => k)*pow((x(j => j)-x(j => k))*(x(j => j)-x(j => k)) + (y(j => j)-y(j => k))*(y(j => j)-y(j => k)), -0.5); } px = m*vx; py = m*vy;