/[ascend]/trunk/models/johnpye/dopri5/dopri5.c
ViewVC logotype

Diff of /trunk/models/johnpye/dopri5/dopri5.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1445 by jpye, Sat May 26 14:43:48 2007 UTC revision 1451 by jpye, Mon May 28 10:20:09 2007 UTC
# Line 102  static double max_d (double a, double b) Line 102  static double max_d (double a, double b)
102    
103  static double hinit (unsigned n, FcnEqDiff *fcn, double x, double* y,  static double hinit (unsigned n, FcnEqDiff *fcn, double x, double* y,
104                       double posneg, double* f0, double* f1, double* yy1, int iord,                       double posneg, double* f0, double* f1, double* yy1, int iord,
105                       double hmax, double* atoler, double* rtoler, int itoler){                       double hmax, double* atoler, double* rtoler, int itoler
106            ,void *user_data
107    ){
108      double   dnf, dny, atoli, rtoli, sk, h, h1, der2, der12, sqr;      double   dnf, dny, atoli, rtoli, sk, h, h1, der2, der12, sqr;
109      unsigned i;      unsigned i;
110    
# Line 141  static double hinit (unsigned n, FcnEqDi Line 143  static double hinit (unsigned n, FcnEqDi
143      /* perform an explicit Euler step */      /* perform an explicit Euler step */
144      for (i = 0; i < n; i++)      for (i = 0; i < n; i++)
145          yy1[i] = y[i] + h * f0[i];          yy1[i] = y[i] + h * f0[i];
146      fcn (n, x+h, yy1, f1);      fcn (n, x+h, yy1, f1, user_data);
147    
148      /* estimate the second derivative of the solution */      /* estimate the second derivative of the solution */
149      der2 = 0.0;      der2 = 0.0;
# Line 176  static double hinit (unsigned n, FcnEqDi Line 178  static double hinit (unsigned n, FcnEqDi
178    
179  /* core integrator */  /* core integrator */
180  static int dopcor (unsigned n, FcnEqDiff *fcn, double x, double* y, double xend,  static int dopcor (unsigned n, FcnEqDiff *fcn, double x, double* y, double xend,
181                     double hmax, double h, double* rtoler, double* atoler,          double hmax, double h, double* rtoler, double* atoler,
182                     int itoler, FILE* fileout, SolTrait *solout, int iout,          int itoler, FILE* fileout, SolTrait *solout, int iout,
183                     long nmax, double uround, int meth, long nstiff, double safe,          long nmax, double uround, int meth, long nstiff, double safe,
184                     double beta, double fac1, double fac2, unsigned* icont){          double beta, double fac1, double fac2, unsigned* icont
185            , void *user_data  
186    ){
187      double   facold, expo1, fac, facc1, facc2, fac11, posneg, xph;      double   facold, expo1, fac, facc1, facc2, fac11, posneg, xph;
188      double   atoli, rtoli, hlamb, err, sk, hnew, yd0, ydiff, bspl;      double   atoli, rtoli, hlamb, err, sk, hnew, yd0, ydiff, bspl;
189      double   stnum, stden, sqr;      double   stnum, stden, sqr;
# Line 224  static int dopcor (unsigned n, FcnEqDiff Line 228  static int dopcor (unsigned n, FcnEqDiff
228      last  = 0;      last  = 0;
229      hlamb = 0.0;      hlamb = 0.0;
230      iasti = 0;      iasti = 0;
231      fcn(n, x, y, k1);      fcn(n, x, y, k1, user_data);
232      hmax = fabs (hmax);      hmax = fabs (hmax);
233      iord = 5;      iord = 5;
234      if (h == 0.0)      if (h == 0.0)
235          h = hinit (n, fcn, x, y, posneg, k1, k2, k3, iord, hmax, atoler, rtoler, itoler);          h = hinit (n, fcn, x, y, posneg, k1, k2, k3, iord, hmax, atoler, rtoler, itoler, user_data);
236      nfcn += 2;      nfcn += 2;
237      reject = 0;      reject = 0;
238      xold = x;      xold = x;
# Line 237  static int dopcor (unsigned n, FcnEqDiff Line 241  static int dopcor (unsigned n, FcnEqDiff
241          irtrn = 1;          irtrn = 1;
242          hout = h;          hout = h;
243          xout = x;          xout = x;
244          solout (naccpt+1, xold, x, y, n, &irtrn);          solout (naccpt+1, xold, x, y, n, &irtrn, user_data);
245          if (irtrn < 0)          if (irtrn < 0)
246          {          {
247              if (fileout)              if (fileout)
# Line 278  static int dopcor (unsigned n, FcnEqDiff Line 282  static int dopcor (unsigned n, FcnEqDiff
282          /* the first 6 stages */          /* the first 6 stages */
283          for (i = 0; i < n; i++)          for (i = 0; i < n; i++)
284              yy1[i] = y[i] + h * a21 * k1[i];              yy1[i] = y[i] + h * a21 * k1[i];
285          fcn (n, x+c2*h, yy1, k2);          fcn (n, x+c2*h, yy1, k2, user_data);
286          for (i = 0; i < n; i++)          for (i = 0; i < n; i++)
287              yy1[i] = y[i] + h * (a31*k1[i] + a32*k2[i]);              yy1[i] = y[i] + h * (a31*k1[i] + a32*k2[i]);
288          fcn (n, x+c3*h, yy1, k3);          fcn (n, x+c3*h, yy1, k3, user_data);
289          for (i = 0; i < n; i++)          for (i = 0; i < n; i++)
290              yy1[i] = y[i] + h * (a41*k1[i] + a42*k2[i] + a43*k3[i]);              yy1[i] = y[i] + h * (a41*k1[i] + a42*k2[i] + a43*k3[i]);
291          fcn (n, x+c4*h, yy1, k4);          fcn (n, x+c4*h, yy1, k4, user_data);
292          for (i = 0; i <n; i++)          for (i = 0; i <n; i++)
293              yy1[i] = y[i] + h * (a51*k1[i] + a52*k2[i] + a53*k3[i] + a54*k4[i]);              yy1[i] = y[i] + h * (a51*k1[i] + a52*k2[i] + a53*k3[i] + a54*k4[i]);
294          fcn (n, x+c5*h, yy1, k5);          fcn (n, x+c5*h, yy1, k5, user_data);
295          for (i = 0; i < n; i++)          for (i = 0; i < n; i++)
296              ysti[i] = y[i] + h * (a61*k1[i] + a62*k2[i] + a63*k3[i] + a64*k4[i] + a65*k5[i]);              ysti[i] = y[i] + h * (a61*k1[i] + a62*k2[i] + a63*k3[i] + a64*k4[i] + a65*k5[i]);
297          xph = x + h;          xph = x + h;
298          fcn (n, xph, ysti, k6);          fcn (n, xph, ysti, k6, user_data);
299          for (i = 0; i < n; i++)          for (i = 0; i < n; i++)
300              yy1[i] = y[i] + h * (a71*k1[i] + a73*k3[i] + a74*k4[i] + a75*k5[i] + a76*k6[i]);              yy1[i] = y[i] + h * (a71*k1[i] + a73*k3[i] + a74*k4[i] + a75*k5[i] + a76*k6[i]);
301          fcn (n, xph, yy1, k2);          fcn (n, xph, yy1, k2, user_data);
302          if (iout == 2){          if (iout == 2){
303              if (nrds == n){              if (nrds == n){
304                  for (i = 0; i < n; i++)                  for (i = 0; i < n; i++)
# Line 418  static int dopcor (unsigned n, FcnEqDiff Line 422  static int dopcor (unsigned n, FcnEqDiff
422              {              {
423                  hout = h;                  hout = h;
424                  xout = x;                  xout = x;
425                  solout (naccpt+1, xold, x, y, n, &irtrn);                  solout (naccpt+1, xold, x, y, n, &irtrn, user_data);
426                  if (irtrn < 0)                  if (irtrn < 0)
427                  {                  {
428                      if (fileout)                      if (fileout)
# Line 464  int dopri5( Line 468  int dopri5(
468      double* atoler, int itoler, SolTrait *solout, int iout, FILE* fileout, double uround,      double* atoler, int itoler, SolTrait *solout, int iout, FILE* fileout, double uround,
469      double safe, double fac1, double fac2, double beta, double hmax, double h,      double safe, double fac1, double fac2, double beta, double hmax, double h,
470      long nmax, int meth, long nstiff, unsigned nrdens, unsigned* icont, unsigned licont      long nmax, int meth, long nstiff, unsigned nrdens, unsigned* icont, unsigned licont
471        , void *user_data
472  ){  ){
473      int       arret, idid;      int       arret, idid;
474      unsigned  i;      unsigned  i;
# Line 660  int dopri5( Line 665  int dopri5(
665      }      }
666      else      else
667      {      {
668          idid = dopcor (n, fcn, x, y, xend, hmax, h, rtoler, atoler, itoler, fileout,          idid = dopcor (n, fcn, x, y, xend, hmax, h, rtoler, atoler, itoler, fileout
669                         solout, iout, nmax, uround, meth, nstiff, safe, beta, fac1, fac2, icont);              , solout, iout, nmax, uround, meth, nstiff, safe, beta, fac1, fac2, icont
670                , user_data
671            );
672          free (ysti);          free (ysti);
673          free (k6);          free (k6);
674          free (k5);    /* reverse order freeing too increase chances */          free (k5);    /* reverse order freeing too increase chances */
Legend:
Removed from v.1445  
changed lines
  Added in v.1451
john.pye@anu.edu.au
 ViewVC Help
Powered by ViewVC 1.1.22