Fixed coefficients in IAS15, thanks to Piotr A. Dybczynski (dybol@amu…

….edu.pl)

rebound/tests/test_interruble_pool.py

@@ -15,13 +15,15 @@ def test_nopool(self):
         pool = InterruptiblePool(2)
         params = [1.,1.1]
         res = [runsim(params[0]), runsim(params[1])]
-        self.assertAlmostEqual(res,[0.9950041652780258,1.095870355119381],delta=1e-15)
+        self.assertAlmostEqual(res[0],0.9950041652780258,delta=1e-15)
+        self.assertAlmostEqual(res[1],1.095870355119381,delta=1e-15)
 
     def test_pool(self):
         pool = InterruptiblePool(2)
         params = [1.,1.1]
         res = pool.map(runsim,params)
-        self.assertAlmostEqual(res,[0.9950041652780258,1.095870355119381],delta=1e-15)
+        self.assertAlmostEqual(res[0],0.9950041652780258,delta=1e-15)
+        self.assertAlmostEqual(res[1],1.095870355119381,delta=1e-15)
 
 if __name__ == "__main__":
     unittest.main()

rebound/tests/test_simulation.py

@@ -88,7 +88,7 @@ def test_calculate_orbits(self):
     def test_com(self):
         self.sim.move_to_com()
         com = self.sim.calculate_com()
-        self.assertEqual(com.x, 0.)
+        self.assertAlmostEqual(com.x, 0., delta=1e-15)
 
     def test_init_megno(self):
         self.sim.init_megno()

src/integrator_ias15.c

@@ -77,11 +77,12 @@ static void predict_next_step(double ratio, int N3,  const struct reb_dpconst7 _
 static const double safety_factor           = 0.25; /**< Maximum increase/deacrease of consecutve timesteps. */
 
 // Gauss Radau spacings
-static const double h[8]    = { 0.0, 0.0562625605369221464656521910, 0.1802406917368923649875799428, 0.3526247171131696373739077702, 0.5471536263305553830014485577, 0.7342101772154105410531523211, 0.8853209468390957680903597629, 0.9775206135612875018911745004}; 
+static const double h[8]    = { 0.0, 0.0562625605369221464656521910318, 0.180240691736892364987579942780, 0.352624717113169637373907769648, 0.547153626330555383001448554766, 0.734210177215410531523210605558, 0.885320946839095768090359771030, 0.977520613561287501891174488626};
 // Other constants
-static const double rr[28] = {0.0562625605369221464656522, 0.1802406917368923649875799, 0.1239781311999702185219278, 0.3526247171131696373739078, 0.2963621565762474909082556, 0.1723840253762772723863278, 0.5471536263305553830014486, 0.4908910657936332365357964, 0.3669129345936630180138686, 0.1945289092173857456275408, 0.7342101772154105410531523, 0.6779476166784883945875001, 0.5539694854785181760655724, 0.3815854601022409036792446, 0.1870565508848551580517038, 0.8853209468390957680903598, 0.8290583863021736216247076, 0.7050802551022034031027798, 0.5326962297259261307164520, 0.3381673205085403850889112, 0.1511107696236852270372074, 0.9775206135612875018911745, 0.9212580530243653554255223, 0.7972799218243951369035946, 0.6248958964481178645172667, 0.4303669872307321188897259, 0.2433104363458769608380222, 0.0921996667221917338008147};
-static const double c[21] = {-0.0562625605369221464656522, 0.0101408028300636299864818, -0.2365032522738145114532321, -0.0035758977292516175949345, 0.0935376952594620658957485, -0.5891279693869841488271399, 0.0019565654099472210769006, -0.0547553868890686864408084, 0.4158812000823068616886219, -1.1362815957175395318285885, -0.0014365302363708915610919, 0.0421585277212687082291130, -0.3600995965020568162530901, 1.2501507118406910366792415, -1.8704917729329500728817408, 0.0012717903090268677658020, -0.0387603579159067708505249, 0.3609622434528459872559689, -1.4668842084004269779203515, 2.9061362593084293206895457, -2.7558127197720458409721005};
-static const double d[21] = {0.0562625605369221464656522, 0.0031654757181708292499905, 0.2365032522738145114532321, 0.0001780977692217433881125, 0.0457929855060279188954539, 0.5891279693869841488271399, 0.0000100202365223291272096, 0.0084318571535257015445000, 0.2535340690545692665214616, 1.1362815957175395318285885, 0.0000005637641639318207610, 0.0015297840025004658189490, 0.0978342365324440053653648, 0.8752546646840910912297246, 1.8704917729329500728817408, 0.0000000317188154017613665, 0.0002762930909826476593130, 0.0360285539837364596003871, 0.5767330002770787313544596, 2.2485887607691598182153473, 2.7558127197720458409721005};
+static const double rr[28] = {0.0562625605369221464656522, 0.1802406917368923649875799, 0.1239781311999702185219278, 0.3526247171131696373739078, 0.2963621565762474909082556, 0.1723840253762772723863278, 0.5471536263305553830014486, 0.4908910657936332365357964, 0.3669129345936630180138686, 0.1945289092173857456275408, 0.7342101772154105315232106, 0.6779476166784883850575584, 0.5539694854785181665356307, 0.3815854601022408941493028, 0.1870565508848551485217621, 0.8853209468390957680903598, 0.8290583863021736216247076, 0.7050802551022034031027798, 0.5326962297259261307164520, 0.3381673205085403850889112, 0.1511107696236852365671492, 0.9775206135612875018911745, 0.9212580530243653554255223, 0.7972799218243951369035945, 0.6248958964481178645172667, 0.4303669872307321188897259, 0.2433104363458769703679639, 0.0921996667221917338008147};
+static const double c[21] = {-0.0562625605369221464656522, 0.0101408028300636299864818, -0.2365032522738145114532321, -0.0035758977292516175949345, 0.0935376952594620658957485, -0.5891279693869841488271399, 0.0019565654099472210769006, -0.0547553868890686864408084, 0.4158812000823068616886219, -1.1362815957175395318285885, -0.0014365302363708915424460, 0.0421585277212687077072973, -0.3600995965020568122897665, 1.2501507118406910258505441, -1.8704917729329500633517991, 0.0012717903090268677492943, -0.0387603579159067703699046, 0.3609622434528459832253398, -1.4668842084004269643701553, 2.9061362593084293014237913, -2.7558127197720458314421588};
+static const double d[21] = {0.0562625605369221464656522, 0.0031654757181708292499905, 0.2365032522738145114532321, 0.0001780977692217433881125, 0.0457929855060279188954539, 0.5891279693869841488271399, 0.0000100202365223291272096, 0.0084318571535257015445000, 0.2535340690545692665214616, 1.1362815957175395318285885, 0.0000005637641639318207610, 0.0015297840025004658189490, 0.0978342365324440053653648, 0.8752546646840910912297246, 1.8704917729329500633517991, 0.0000000317188154017613665, 0.0002762930909826476593130, 0.0360285539837364596003871, 0.5767330002770787313544596, 2.2485887607691597933926895, 2.7558127197720458314421588};
+
 
 // Weights for integration of a first order differential equation (Note: interval length = 2) 
 static const double w[8] = {0.03125, 0.185358154802979278540728972807180754479812609, 0.304130620646785128975743291458180383736715043, 0.376517545389118556572129261157225608762708603, 0.391572167452493593082499533303669362149363727, 0.347014795634501068709955597003528601733139176, 0.249647901329864963257869294715235590174262844, 0.114508814744257199342353731044292225247093225};
@@ -736,13 +737,13 @@ void integrator_generate_constants(void){
         mpf_init(_d[i]);
     }
     mpf_set_str(_h[0],"0.0",10);
-    mpf_set_str(_h[1],"0.0562625605369221464656521910",10);
-    mpf_set_str(_h[2],"0.1802406917368923649875799428",10);
-    mpf_set_str(_h[3],"0.3526247171131696373739077702",10);
-    mpf_set_str(_h[4],"0.5471536263305553830014485577",10);
-    mpf_set_str(_h[5],"0.7342101772154105410531523211",10);
-    mpf_set_str(_h[6],"0.8853209468390957680903597629",10);
-    mpf_set_str(_h[7],"0.9775206135612875018911745004",10);
+    mpf_set_str(_h[1],"0.0562625605369221464656521910318",10);
+    mpf_set_str(_h[2],"0.180240691736892364987579942780",10);
+    mpf_set_str(_h[3],"0.352624717113169637373907769648",10);
+    mpf_set_str(_h[4],"0.547153626330555383001448554766",10);
+    mpf_set_str(_h[5],"0.734210177215410531523210605558", 10);
+    mpf_set_str(_h[6],"0.885320946839095768090359771030",10);
+    mpf_set_str(_h[7],"0.977520613561287501891174488626",10);
 
     int l=0;
     for (int j=1;j<8;++j) {