Re: [Help-gsl] differential equation system - problem

[Pau Cervera Badia]

The function

gsl_odeiv_evolve_apply (e, c, s, &sys, &t, t_end, &h, y) evolves the system 
from the actual time t to t_end. You have supplied t_end = 1 = t, so the system does not 
evolve at all. If you change this the code runs.




Christian Ost wrote:

> 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 |
>
>
> _______________________________________________
> Help-gsl mailing list
> address@hidden
> http://lists.gnu.org/mailman/listinfo/help-gsl
>
>
>  

--
Pau Cervera i Badia       e-mail: address@hidden

***************************************************
Departament de Física Fonamental
Universitat de Barcelona
Martí i Franqués, 1
Planta 3, despatx 346 bis
08028 Barcelona
Spain

tel: +34 934 921 155
***************************************************
"To err is human, but to really foul things up requires a computer."
http://www.outpost9.com/reference/jargon/jargon_26.html#TAG996