[Pomp-commits] r44 - in pkg: . src

[noreply at r-forge.r-project.org]

Author: kingaa
Date: 2008-08-15 12:08:00 +0200 (Fri, 15 Aug 2008)
New Revision: 44

Modified:
   pkg/DESCRIPTION
   pkg/src/euler.c
Log:
fix bug in Euler step

Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION	2008-08-14 16:48:12 UTC (rev 43)
+++ pkg/DESCRIPTION	2008-08-15 10:08:00 UTC (rev 44)
@@ -1,8 +1,8 @@
 Package: pomp
 Type: Package
 Title: Statistical inference for partially observed Markov processes
-Version: 0.21-2
-Date: 2008-08-14
+Version: 0.21-3
+Date: 2008-08-15
 Author: Aaron A. King, Edward L. Ionides, Carles Martinez Breto
 Maintainer: Aaron A. King <kingaa at umich.edu>
 Description: Inference methods for partially-observed Markov processes

Modified: pkg/src/euler.c
===================================================================
--- pkg/src/euler.c	2008-08-14 16:48:12 UTC (rev 43)
+++ pkg/src/euler.c	2008-08-15 10:08:00 UTC (rev 44)
@@ -20,10 +20,12 @@
   int nzero = ndim[6];
   double covar_fn[covdim];
   int j, k, p, step, neuler;
-  double dt;
+  double dt, tol;
 
   struct lookup_table covariate_table = {covlen, covdim, 0, time_table, covar_table};
   
+  tol = sqrt(DOUBLE_EPS); // relative tolerance in choosing Euler stepsize
+
   // copy the start values into the result array
   for (p = 0; p < nrep; p++)
     for (k = 0; k < nvar; k++) 
@@ -36,18 +38,22 @@
 
     t = times[step-1];
     dt = *deltat;
+
     // neuler is the number of Euler steps to take in going from
     // times[step-1] to times[step].
     // this choice is meant to be conservative
-    // (i.e., so that the actual dt does not exceed the specified dt)
-    // (at least not by much) and to counteract roundoff error.
-    // it is not guaranteed: suggestions would be appreciated.
-    neuler = (int) ceil(((times[step]-t)/dt)-2*(DOUBLE_EPS/dt));
+    // (i.e., so that the actual dt does not exceed the specified dt 
+    // by more than the relative tolerance 'tol')
+    // and to counteract roundoff error.
+    // It seems to work well, but is not guaranteed: 
+    // suggestions would be appreciated.
 
-    // if the next step would carry us beyond times[step], reduce dt
     if (t+dt >= times[step]) {
       dt = times[step] - t; 
       neuler = 1;
+    } else {
+      neuler = (int) ceil((times[step]-t)/dt/(1+tol));
+      dt = (times[step]-t)/((double) neuler);
     }
 
     for (p = 0; p < nrep; p++) {
@@ -75,11 +81,11 @@
 
       }
 
+      t += dt;
+
       if (j == neuler-2) {	// penultimate step
 	dt = times[step]-t;
 	t = times[step]-dt;
-      } else {
-	t += dt;
       }
 
     }