boost - odeint

boost - odeint - C++ - Programmation

Marsh Posté le 06-07-2014 à 14:24:01    

Bonjour,
J'ai un petit souci que je n'arrive pas à résoudre avec la bibliothèque odeint (qui permet de résoudre un système d'équations différentielles) de boost (version 1.55).
En fait, j'arrive à résoudre le système, mais je n'obtiens pas le même résultat sous linux et sous windows (les deux étant en 64 bits). Sous linux, le programme sort même parfois la valeur NaN alors que ça fonctionne très bien sous windows.
Voilà un exemple de code qui ne fonctionne pas :

Code :
  1. #include <iostream>
  2. #include <stdlib.h>
  3. #include <boost/numeric/odeint.hpp>
  5. #define eps0 8.85418782E-12
  6. #define mu0 1.25663706E-6
  8. using namespace std;
  9. using namespace boost::numeric::odeint;
  11. typedef boost::array< double , 4 > state_type;
  12. typedef controlled_runge_kutta< runge_kutta_cash_karp54< state_type > > stepper_type;
  14. class equation{
  15. public:
  16.     double k0;
  17.     double neff;
  18.     double nu;
  19.     equation( double ak0) {k0=ak0;neff=1.4602;nu=1.0;}
  20.     void operator() ( const state_type &y , state_type &dydt , const double  r  )
  21.     {
  22.       double n=0.03*(1-pow(r/2E-6,2))+1.45;
  23. double np=-0.03*2*r/pow(2E-6,2);
  24. double p=pow(k0*n,2)-pow(k0*neff,2);
  25.       dydt[0] = y[2];
  26.       dydt[1] = y[3];
  27.       dydt[2] = y[0]*(pow(nu/r,2)-p)+y[2]*(2*pow(k0*neff,2)/p*np/n-1/r)-y[1]*sqrt(mu0/eps0)*2.0*nu*k0*neff*k0/r/p*np/n;
  28.       dydt[3] = y[1]*(pow(nu/r,2)-p)+y[3]*(2*k0*k0*n*np/p-1/r)-y[0]*sqrt(eps0/mu0)*2*nu*k0*neff*k0*n*np/r/p;
  29.     }
  30. };
  32. void write_cout(const state_type &x, const double t)
  33. {
  34. cout << t << endl;
  35. cout << x[0] << " ; " << x[1] << " ; " << x[2] << " ; " << x[3] << endl;
  36. }
  38. int main()
  39. {
  40. double k0=2*M_PI/1E-6;
  41.     state_type y;
  42.     y[0]=7.58003E-7;
  43.     y[1]=y[0];
  44.     y[2]=758003.0;
  45.     y[3]=y[2];
  46.     equation eq1(k0);
  47.     cout << y[0]<< " ; "  << y[1] << " ; " << y[2] <<  "; " << y[3] << " ; " << eq1.neff << " ; " << eq1.nu << " ; " << eq1.k0 << endl;
  48.     integrate_adaptive(stepper_type(default_error_checker< double >( 1E-15 , 1E-15 ) ),eq1,y,1E-12,2E-6,2E-10,write_cout);
  49.     cout << y[0]<< " ; "  << y[1] << " ; " << y[2] <<  "; " << y[3] << " ; " << eq1.neff << " ; " << eq1.nu << " ; " << eq1.k0 << endl;
  50.     return 0;
  51. }


 
Pour la compilation, je ne mets aucune option particulière mis à part le chemin vers boost.
Sous linux, j'ai essayé avec gcc 4.1,4.4 et 4.8 et aussi icc, le résultat est le même : y[0] vaut NaN à la fin.
Sous windows j'ai compilé avec mingw 4.7 et 4.8 et y[0] vaut environ 105.
Si je mets 1E-14 en précision (au lieu de 1E-15), ça fonctionne aussi sous linux, mais je n'obtiens pas tout-à fait les mêmes valeurs : par exemple y[4] vaut 542116.37 sous linux et 542116.52 sous windows... Idem avec 1E-6.
 
 
Quelqu'un aurait-il une idée sur la question?
Merci!

Reply

Marsh Posté le 06-07-2014 à 14:24:01   

Reply

Sujets relatifs:

Leave a Replay

Make sure you enter the(*)required information where indicate.HTML code is not allowed