Basic question about GSL ODE func RK4
Hi,
I am looking for some help on rk4 and coulomb repulsion.
It is a technical question but one has to go through the story to understand what I am tyring to achieve. Thank you already.
I am trying to simulate the Coulomb repulsion between charged particles. I am using RK4 to do this which requires a function. The physics formula of the function is:
dx/dt = k*qi*(xi-xj)/sqrt((xi-xj)^2+(yi-yj)^2+(zi-zj)^2)
where qi is the charge of the i-th particle, xi, yi, zi are the
coordinates of the i-th particle and xj, yj,zj are the coordinates of the j-th
particle.
As an input for the initial values of the positions and charge for each particle, I am using arrays (currently only for 4 particles):
q[4] = { 1.0, 1.0, -1.0, 1.0 } charge
x0[4] = { 1.0, 2.3, 3.2, 4.0 } init posn x
y0[4] = { 1.0, 2.4, 3.0, 4.0 } init posn y
z0[4] = { 1.0, 2.5, 3.3, 4.2 } init posn z
and by using two "for" loops I can loop through each particle and compute the value of the total effect of all the particles on eachother. My question is: since there are two values for each coordinate (i and j), how can I write the function in an acceptable form for RK4?
I did try the following:
Code:
Code:
int func (double t, const double y[], double f[], void *params){
double mu = *(double *)params;
double dx = (y[0]-y[1])/r;
double dy = (y[2]-y[3])/r;
double dz = (y[4]-y[5])/r;
f[0] = mu*y[6]*dx/(r*r);
f[1] = mu*y[6]*dy/(r*r);
f[2] = mu*y[6]*dz/(r*r);
return GSL_SUCCESS;
}
however, I am not convinced this is correct (at all).
I would appreciate any feedback.
Cheers.
