Re: [Help-gsl] problem

[caoj]

Hi Pau Cervera Badia:

Thank you for your response. I check my program again. It stuck in  t=5.790250, 
t1=6.467000.  After gsl_odeiv_evolve_apply, t always is 5.790250 and does not 
change. So the program stuck in the "while" loop. I try to change h to 1e-2 and 
 1e-10. I get the same result. Do you have anoter idea? Thanks again.

Jun Cao


Quoting Pau Cervera Badia <address@hidden>:

> gsl_odeiv_evolve_apply(e,c,s,&sys,&t,t1,&h,y) will advance from time t 
> to the next integration time with an optimum step-size, provided that 
> the next time is less than t1. Otherwise it will integrate the system to 
> t1 exactly. So after some iterations the while condition will be false.
> 
> You can check that gsl_odeiv_evolve_apply is running properly with:
> 
> while (t < t1) {
> 
>     int status = gsl_odeiv_evolve_apply(e,c,s,&sys,&t,t1,&h,y);
> 
>     printf("%.5e %.5e\n", t, t1); // this should print t less than t1
> 
>     if (status != GSL_SUCCESS)
>        break;
> 
> }
> printf("%.5e %.5e\n",t,t1); // this should print t equal to t1
> 
> If I'm not missing something, your program should work. Maybe the 
> problem is that t is always less than t1 for some other reason (maybe 
> the routine couldn't reach the prescribed precision and 
> gsl_odeiv_evolve_apply is trying to reduce the step-size more and more). 
> Maybe you can check the h values.
> 
> Maybe something will be of any help.
> 
> address@hidden wrote:
> 
> >Hi:
> >
> >I use ordinary differential equations from GSL. The part of my program is
> like 
> >example shown in the document.
> >
> >  const gsl_odeiv_step_type * T 
> >    = gsl_odeiv_step_rk8pd;
> >
> >  gsl_odeiv_step * s 
> >    = gsl_odeiv_step_alloc (T, 9);
> >  gsl_odeiv_control * c 
> >    = gsl_odeiv_control_y_new (1e-6, 0.0);
> >  gsl_odeiv_evolve * e 
> >    = gsl_odeiv_evolve_alloc (9);
> >
> >  double mu = 10;
> >  gsl_odeiv_system sys = {func, NULL, 9, &mu};
> >
> >  double t = 0.0;  
> >  double h = 1e-6;
> >
> >  for (i = 1; i <no_of_data_points; i++) {
> >     t1 = mytable[i];
> >     while (t < t1)
> >     {
> >        int status = gsl_odeiv_evolve_apply (e, c, s,
> >                                           &sys, 
> >                                           &t, t1,
> >                                           &h, y);
> >
> >      if (status != GSL_SUCCESS)
> >          break;
> >      }
> >      ...
> >   }
> >  gsl_odeiv_evolve_free (e);
> >  gsl_odeiv_control_free (c);
> >  gsl_odeiv_step_free (s);
> >
> >The weird thing is that at some i, the program stuck inside of "while" loop
> and 
> >never come out. gsl_odeiv_evolve_apply will give t=t, so t is always < t1.
> As 
> >result the program never stop. Dose anyone know how to solve this problem? 
> >Thanks.
> >
> >Jun Cao
> >
> >
> >
> >_______________________________________________
> >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                 Martí i Franqués, 1
>    Universitat de Barcelona                   Planta 3, despatx 346 bis
>                                                         08028 Barcelona
>    tel: +34 934 921 155                                           Spain
> 
>    "To err is human, but to really foul things up requires a computer."
> }
> 
> 
>