[Help-gsl] differential equation system - problem

[Christian Ost]

Hi,

currently i am trying to get familiar with GSL and as a beginning I
tried to get a numerical
solution, of the gravitational attraction and movement between two masses.

i.e.:   [latex]\ddot{r}_a=-G*m_b*\frac{1}{(r_a-r_b)^2} [/latex]
and   [latex]\ddot{r}_b=-G*m_a*\frac{1}{(r_b-r_a)^2} [/latex].

I tried to solve it as follows:

#include <gsl/gsl_sys.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>
#include <stdio.h>

int Func (double t, const double y[], double f[], void *pars) {

    double G=1;
    double m_a=3;
    double m_b=1;

    f[0]=y[1];
    f[1]=-G*m_b/((y[0]-y[2])*(y[0]-y[2]));

    f[2]=y[3];
    f[3]=-G*m_a/((y[2]-y[0])*(y[2]-y[0]));

    return GSL_SUCCESS;
}

int main() {

    double t=1;
    const gsl_odeiv_step_type * T =gsl_odeiv_step_rkf45;


    gsl_odeiv_step * s    = gsl_odeiv_step_alloc (T, 4);
    gsl_odeiv_control * c = gsl_odeiv_control_y_new (0,1.e-15);
    gsl_odeiv_evolve * e  = gsl_odeiv_evolve_alloc (4);


    gsl_odeiv_system sys = {Func, NULL, 4, NULL};

    double h = 1e-6;
    double y [4];
    int evolve=1000;

    y[0]=1;
    y[1]=0;

    y[2]=1.1;
    y[3]=0;

    while (evolve) {
        int status = gsl_odeiv_evolve_apply (e, c, s,
                                             &sys,
                                             &t, 1,
                                             &h, y);

        if (status != GSL_SUCCESS) {
            return 0;
        }

        printf ("%f - %f\n",y[2], y[0]);
        evolve--;
    }

    gsl_odeiv_evolve_free(e);
    gsl_odeiv_control_free(c);
    gsl_odeiv_step_free(s);
    
    return EXIT_SUCCESS;
}

but none of the two two masses changes its position. 

What did I do wrong?

Thanks for answering this (probably very stupid) question

Christian Ost

--

Christian Ost | address@hidden | https://lg3d-livecd.dev.java.net |