RACIPELIB.c
// solve the ODE by RK-45
solve the ODE by RK-45
solve the ODE by RK - 45
void solve_ODE_rk45 (int j, struct opts *simu_opts, struct topo *topoinfo, struct rlt *tmprlt)
{
  double abserr;
  double relerr;
  int    flag;
  int    n_step     = 100;
  int    i_step     = 1;
  int    i          = 0;
  double testdelta  = 0.0;
  double t          = 0.0;
  double t_out      = 0.0;
  double t_start    = 0.0;
  double t_stop     = 0.0;
  double *y;
  double *yp;
  double *ytmp;
  abserr = r4_epsilon ( ); 
  relerr = r4_epsilon ( ); 
  y     = (double *)calloc(topoinfo->numG,      sizeof(double));
  yp    = (double *)calloc(topoinfo->numG,      sizeof(double));
  ytmp  = (double *)calloc(topoinfo->numG,      sizeof(double));
  for (i = 0; i < topoinfo->numG; i++){
    ytmp[i] = 2000.0;
  }
  int cnt_loop = 0;
  RIVs(y, ytmp, tmprlt->paras, topoinfo);
  testdelta = sumdelta(y, ytmp, topoinfo->numG);
  while (testdelta != 0 && cnt_loop < simu_opts->maxiters) {
    t_start = t_stop;
    t_stop  = t_stop + 100;
    cnt_loop = cnt_loop + 1;
    model_ODE ( t, ytmp, yp, tmprlt->paras, topoinfo);
    flag = 1;
    for ( i_step = 1; i_step <= n_step; i_step++ )
    { 
        for (i = 0; i < topoinfo->numG; i++){
            ytmp[i] = y[i];
        }
        t = ( ( double ) ( n_step - i_step + 1 ) * t_start  
            + ( double ) (          i_step - 1 ) * t_stop )
            / ( double ) ( n_step              );
        t_out = ( ( double ) ( n_step - i_step ) * t_start  
                + ( double ) (          i_step ) * t_stop )
                / ( double ) ( n_step          );
        flag = r4_rkf45 (model_ODE, topoinfo->numG, y, yp, tmprlt->paras, topoinfo, &t, t_out, &relerr, abserr, flag );
    }
    testdelta = sumdelta(y, ytmp, topoinfo->numG);
  }
  for (i = 0; i < topoinfo->numG; i++){
    tmprlt->y_store[topoinfo->numG*j + i] = y[i];
  }
  free(y);
  free(yp);
  free(ytmp);
}
void solve_ODE_rk45 ( int j , struct opts * simu_opts , struct topo * topoinfo , struct rlt * tmprlt ) { double abserr ; double relerr ; int flag ; int n_step = 100 ; int i_step = 1 ; int i = 0 ; double testdelta = 0.0 ; double t = 0.0 ; double t_out = 0.0 ; double t_start = 0.0 ; double t_stop = 0.0 ; double * y ; double * yp ; double * ytmp ; abserr = r4_epsilon ( ) ; relerr = r4_epsilon ( ) ; y = ( double * ) calloc ( topoinfo -> numG , sizeof ( double ) ) ; yp = ( double * ) calloc ( topoinfo -> numG , sizeof ( double ) ) ; ytmp = ( double * ) calloc ( topoinfo -> numG , sizeof ( double ) ) ; for ( i = 0 ; i < topoinfo -> numG ; i ++ ) { ytmp [ i ] = 2000.0 ; } int cnt_loop = 0 ; RIVs ( y , ytmp , tmprlt -> paras , topoinfo ) ; testdelta = sumdelta ( y , ytmp , topoinfo -> numG ) ; while ( testdelta != 0 && cnt_loop < simu_opts -> maxiters ) { t_start = t_stop ; t_stop = t_stop + 100 ; cnt_loop = cnt_loop + 1 ; model_ODE ( t , ytmp , yp , tmprlt -> paras , topoinfo ) ; flag = 1 ; for ( i_step = 1 ; i_step <= n_step ; i_step ++ ) { for ( i = 0 ; i < topoinfo -> numG ; i ++ ) { ytmp [ i ] = y [ i ] ; } t = ( ( double ) ( n_step - i_step + 1 ) * t_start + ( double ) ( i_step - 1 ) * t_stop ) / ( double ) ( n_step ) ; t_out = ( ( double ) ( n_step - i_step ) * t_start + ( double ) ( i_step ) * t_stop ) / ( double ) ( n_step ) ; flag = r4_rkf45 ( model_ODE , topoinfo -> numG , y , yp , tmprlt -> paras , topoinfo , & t , t_out , & relerr , abserr , flag ) ; } testdelta = sumdelta ( y , ytmp , topoinfo -> numG ) ; } for ( i = 0 ; i < topoinfo -> numG ; i ++ ) { tmprlt -> y_store [ topoinfo -> numG * j + i ] = y [ i ] ; } free ( y ) ; free ( yp ) ; free ( ytmp ) ; }
solve the ODE by RK-45
solve the ODE by RK - 45
//sqrt ( r4_epsilon ( ) ); //sqrt ( r4_epsilon ( ) ); // printf("RACIPELIB 2008"); // printf ( "%4d  %12f\n", flag, t);  // if (i_step == n_step - 1) { //     printf ( "%4d  %12f  %12f  %12f\n", flag, t, y[0], y[1]);  // } // printf ( "%4d  %12f  %12f  %12f\n", flag, t, y[0], y[1]);