[Boost-commit] svn:boost r82749 - trunk/libs/math/test

Subject: [Boost-commit] svn:boost r82749 - trunk/libs/math/test
From: pbristow_at_[hidden]
Date: 2013-02-06 20:26:16

Author: pbristow
Date: 2013-02-05 07:59:09 EST (Tue, 05 Feb 2013)
New Revision: 82749
URL: http://svn.boost.org/trac/boost/changeset/82749

Log:
More bessel and Lots more neumann tests.  I failure at double epsilon airy TODO.
Text files modified: 
   trunk/libs/math/test/test_bessel_airy_zeros.cpp |   218 ++++++++++++++++++++++++++++++++++++--- 
   1 files changed, 197 insertions(+), 21 deletions(-)

Modified: trunk/libs/math/test/test_bessel_airy_zeros.cpp
==============================================================================
--- trunk/libs/math/test/test_bessel_airy_zeros.cpp	(original)
+++ trunk/libs/math/test/test_bessel_airy_zeros.cpp	2013-02-05 07:59:09 EST (Tue, 05 Feb 2013)
@@ -36,10 +36,40 @@
 // using the online special function calculator at functions.wolfram.com,
 // and values generated with Boost.Multiprecision at about 1000-bit or 100 decimal digits precision.
 
+// We are most grateful for the invaluable
 // Weisstein, Eric W. "Bessel Function Zeros." From MathWorld--A Wolfram Web Resource.
 // http://mathworld.wolfram.com/BesselFunctionZeros.html 
+// and the newer http://www.wolframalpha.com/
 
-// See also http://dlmf.nist.gov/10.21
+// See also NIST Handbook of Mathetmatical Function http://dlmf.nist.gov/10.21
+/*
+
+The algorithms for estimating the roots of both cyl. Bessel
+as well as cyl. Neumann have the same cross-over points,
+and also use expansions that have the same order of approximation.
+
+Therefore, tests will be equally effective for both functions
+in the regions of order.
+
+I have recently changed a critical cross-over in the algorithms
+from a value of order of 1.2 to a value of order of 2.2.
+In addition, there is a critical cross-over in the rank of the
+zero from rank 1 to rank 2 and above. The first zero is
+treated differently than the remaining ones.
+
+So I would be most interested in various regions of order,
+each one tested with about 20 zeros should suffice:
+•	Order 219/100: This checks a region just below a critical cutoff
+•	Order 221/100: This checks a region just above a critical cutoff
+•	Order 0: Something always tends to go wrong at zero.
+•	Order 1/1000: A small order
+•	Order 71/19: Merely an intermediate order.
+•	Order 7001/19: A medium-large order, small enough to retain moderate efficiency of calculation.
+ 
+If we would like, we could add a few selected high zeros
+such as the 1000th zero for a few modest orders such as 71/19, etc.
+
+*/
 
 template <class RealType>
 void test_bessel_zeros(RealType)
@@ -98,21 +128,43 @@
   BOOST_CHECK_CLOSE_FRACTION(cyl_bessel_j_zero(static_cast<RealType>(0), 4U), static_cast<RealType>(11.791534439014281613743044911925458922022924699695L), tolerance);
   BOOST_CHECK_CLOSE_FRACTION(cyl_bessel_j_zero(static_cast<RealType>(0), 5U), static_cast<RealType>(14.930917708487785947762593997388682207915850115633L), tolerance);
 
-  {
+  { // Same test using the multiple zeros version.
     std::vector<RealType> zeros;
-  
     cyl_bessel_j_zero(static_cast<RealType>(0.0), 3, 1U, std::back_inserter(zeros) );
-
     BOOST_CHECK_CLOSE_FRACTION(zeros[0], static_cast<RealType>(2.4048255576957727686216318793264546431242449091460L), tolerance);
     BOOST_CHECK_CLOSE_FRACTION(zeros[1], static_cast<RealType>(5.5200781102863106495966041128130274252218654787829L), tolerance);
     BOOST_CHECK_CLOSE_FRACTION(zeros[2], static_cast<RealType>(8.6537279129110122169541987126609466855657952312754L), tolerance);
+  }
+  // 1/1000 a small order.
+  /* Table[N[BesselJZero[1/1000, n], 50], {n, 1, 4, 1}]
+  n | 
+1 | 2.4063682720422009275161970278295108254321633626292
+2 | 5.5216426858401848664019464270992222126391378706092
+3 | 8.6552960859298799453893840513333150237193779482071
+4 | 11.793103797689738596231262077785930962647860975357
 
+Table[N[BesselJZero[1/1000, n], 50], {n, 10, 20, 1}]
+n | 
+10 | 30.636177039613574749066837922778438992469950755736
+11 | 33.777390823252864715296422192027816488172667994611
+12 | 36.918668992567585467000743488690258054442556198147
+13 | 40.059996426251227493370316149043896483196561190610
+14 | 43.201362392820317233698309483240359167380135262681
+15 | 46.342759065846108737848449985452774243376260538634
+16 | 49.484180603489984324820981438067325210499739716337
+17 | 52.625622557085775090390071484188995092211215108718
+18 | 55.767081479279692992978326069855684800673801918763
+19 | 58.908554657366270044071505013449016741804538135905
+20 | 62.050039927521244984641179233170843941940575857282
 
-     //std::copy(zeros.begin(), zeros.end(),
-     //           std::ostream_iterator<double>(std::cout, "\n"));
-  }
+*/
 
-  /*
+  BOOST_CHECK_CLOSE_FRACTION(cyl_bessel_j_zero(static_cast<RealType>(1)/1000, 1U), static_cast<RealType>(2.4063682720422009275161970278295108254321633626292L), tolerance);
+  BOOST_CHECK_CLOSE_FRACTION(cyl_bessel_j_zero(static_cast<RealType>(1)/1000, 4U), static_cast<RealType>(11.793103797689738596231262077785930962647860975357L), tolerance);
+  BOOST_CHECK_CLOSE_FRACTION(cyl_bessel_j_zero(static_cast<RealType>(1)/1000, 10U), static_cast<RealType>(30.636177039613574749066837922778438992469950755736L), tolerance);
+  BOOST_CHECK_CLOSE_FRACTION(cyl_bessel_j_zero(static_cast<RealType>(1)/1000, 20U), static_cast<RealType>(62.050039927521244984641179233170843941940575857282L), tolerance);
+
+    /*
   Table[N[BesselJZero[1, n], 50], {n, 1, 4, 1}]
   n | 
   1 | 3.8317059702075123156144358863081607665645452742878
@@ -142,15 +194,62 @@
   BOOST_CHECK_CLOSE_FRACTION(cyl_bessel_j_zero(static_cast<RealType>(5), 4U), static_cast<RealType>(18.980133875179921120770736748466932306588828411497L), tolerance);
   BOOST_CHECK_CLOSE_FRACTION(cyl_bessel_j_zero(static_cast<RealType>(5), 5U), static_cast<RealType>(22.217799896561267868824764947529187163096116704354L), tolerance);
 
-  // Some none integral tests.
-  BOOST_CHECK_CLOSE_FRACTION(cyl_bessel_j_zero(static_cast<RealType>(3.736842105263157894736842105263157894736842105263157894736842105263157894736842105263157894736842105L), 1U), static_cast<RealType>(7.273175193831648950318569426229076558896319670162279791988152000556091140599946365217211157877052381L), tolerance);
-  BOOST_CHECK_CLOSE_FRACTION(cyl_bessel_j_zero(static_cast<RealType>(3.736842105263157894736842105263157894736842105263157894736842105263157894736842105263157894736842105L), 20U), static_cast<RealType>(67.81514561969629092555679137555595116511146058545787883557679231060644931096494584364894743334132014L), tolerance);
+  // An intermediate order
+  /*
+  Table[N[BesselJZero[71/19, n], 50], {n, 1, 20, 1}]
+
+  7.27317519383164895031856942622907655889631967016227,
+  10.7248583088831417325361727458514166471107495990849, 
+  14.0185045994523881061204595580426602824274719315813, 
+  17.2524984591704171821624871665497773491959038386104, 
+  20.4566788740445175951802340838942858854605020778141, 
+  23.6436308971423452249455142271473195998540517250404, 
+  26.8196711402550877454213114709650192615223905192969, 
+  29.9883431174236747426791417966614320438788681941419, 
+  33.1517968976905208712508624699734452654447919661140, 
+  36.3114160002162074157243540350393860813165201842005, 
+  39.4681324675052365879451978080833378877659670320292, 
+  42.6225978013912364748550348312979540188444334802274, 
+  45.7752814645368477533902062078067265814959500124386, 
+  48.9265304891735661983677668174785539924717398947994, 
+  52.0766070453430027942797460418789248768734780634716, 
+  55.2257129449125713935942243278172656890590028901917, 
+  58.3740061015388864367751881504390252017351514189321, 
+  61.5216118730009652737267426593531362663909441035715, 
+  64.6686310537909303683464822148736607945659662871596, 
+  67.8151456196962909255567913755559511651114605854579
+  */
+  BOOST_CHECK_CLOSE_FRACTION(cyl_bessel_j_zero(static_cast<RealType>(71)/19, 1U), static_cast<RealType>(7.27317519383164895031856942622907655889631967016227L), tolerance);
+  BOOST_CHECK_CLOSE_FRACTION(cyl_bessel_j_zero(static_cast<RealType>(71)/19, 4U), static_cast<RealType>(17.2524984591704171821624871665497773491959038386104L), tolerance);
+  BOOST_CHECK_CLOSE_FRACTION(cyl_bessel_j_zero(static_cast<RealType>(71)/19, 10U), static_cast<RealType>(36.3114160002162074157243540350393860813165201842005L), tolerance);
+  BOOST_CHECK_CLOSE_FRACTION(cyl_bessel_j_zero(static_cast<RealType>(71)/19, 20U), static_cast<RealType>(67.8151456196962909255567913755559511651114605854579L), tolerance);
+  /*
+
+  Table[N[BesselJZero[7001/19, n], 50], {n, 1, 2, 1}]
 
-  // Some none integral tests in 'tough' regions.
+1 | 381.92201523024489386917204470434842699154031135348
+2 | 392.17508657648737502651299853099852567001239217724
+
+Table[N[BesselJZero[7001/19, n], 50], {n, 19, 20, 1}]
+ 
+19 | 491.67809669154347398205298745712766193052308172472
+20 | 496.39435037938252557535375498577989720272298310802
+
+  */
+  BOOST_CHECK_CLOSE_FRACTION(cyl_bessel_j_zero(static_cast<RealType>(7001)/19, 1U), static_cast<RealType>(381.92201523024489386917204470434842699154031135348L), tolerance);
+  BOOST_CHECK_CLOSE_FRACTION(cyl_bessel_j_zero(static_cast<RealType>(7001)/19, 2U), static_cast<RealType>(392.17508657648737502651299853099852567001239217724L), tolerance);
+  BOOST_CHECK_CLOSE_FRACTION(cyl_bessel_j_zero(static_cast<RealType>(7001)/19, 20U), static_cast<RealType>(496.39435037938252557535375498577989720272298310802L), tolerance);
+
+  // Some non-integral tests.
+  BOOST_CHECK_CLOSE_FRACTION(cyl_bessel_j_zero(static_cast<RealType>(3.73684210526315789473684210526315789473684210526315789L), 1U), static_cast<RealType>(7.273175193831648950318569426229076558896319670162279791988152000556091140599946365217211157877052381L), tolerance);
+  BOOST_CHECK_CLOSE_FRACTION(cyl_bessel_j_zero(static_cast<RealType>(3.73684210526315789473684210526315789473684210526315789L), 20U), static_cast<RealType>(67.81514561969629092555679137555595116511146058545787883557679231060644931096494584364894743334132014L), tolerance);
+
+  // Some non-integral tests in 'tough' regions. 
+  // Order 219/100: This checks a region just below a critical cutoff.
   BOOST_CHECK_CLOSE_FRACTION(cyl_bessel_j_zero(static_cast<RealType>(219)/100, 1U), static_cast<RealType>(5.37568854370623186731066365697341253761466705063679L), tolerance);
   BOOST_CHECK_CLOSE_FRACTION(cyl_bessel_j_zero(static_cast<RealType>(219)/100, 2U), static_cast<RealType>(8.67632060963888122764226633146460596009874991130394L), tolerance);
   BOOST_CHECK_CLOSE_FRACTION(cyl_bessel_j_zero(static_cast<RealType>(219)/100, 20U), static_cast<RealType>(65.4517712237598926858973399895944886397152223643028L), tolerance);
-
+  // Order 221/100: This checks a region just above a critical cutoff.
   BOOST_CHECK_CLOSE_FRACTION(cyl_bessel_j_zero(static_cast<RealType>(221)/100, 1U), static_cast<RealType>(5.40084731984998184087380740054933778965260387203942L), tolerance);
   BOOST_CHECK_CLOSE_FRACTION(cyl_bessel_j_zero(static_cast<RealType>(221)/100, 2U), static_cast<RealType>(8.70347906513509618445695740167369153761310106851599L), tolerance);
   BOOST_CHECK_CLOSE_FRACTION(cyl_bessel_j_zero(static_cast<RealType>(221)/100, 20U), static_cast<RealType>(65.4825314862621271716158606625527548818843845600782L), tolerance);
@@ -176,7 +275,8 @@
      BOOST_CHECK_THROW(cyl_bessel_j_zero(static_cast<RealType>(std::numeric_limits<RealType>::infinity()), 1U), std::domain_error);
   }
   
-  // Tests of cyc_neumann function
+  // Tests of cyc_neumann zero function (BesselKZero in Wolfram).
+
   /*
   Table[N[BesselKZero[0, n], 50], {n, 1, 5, 1}]
 n | 
@@ -203,27 +303,103 @@
 5 | 16.378966558947456561726714466123708444627678549687
 
 */
+  // Some simple integer values.
 
   using boost::math::cyl_neumann_zero;
 
   BOOST_CHECK_CLOSE_FRACTION(cyl_neumann_zero(static_cast<RealType>(0), 1U), static_cast<RealType>(0.89357696627916752158488710205833824122514686193001L), tolerance);
   BOOST_CHECK_CLOSE_FRACTION(cyl_neumann_zero(static_cast<RealType>(1), 2U), static_cast<RealType>(5.4296810407941351327720051908525841965837574760291L), tolerance);
   BOOST_CHECK_CLOSE_FRACTION(cyl_neumann_zero(static_cast<RealType>(2), 3U), static_cast<RealType>(10.023477979360037978505391792081418280789658279097L), tolerance);
-
   
-  {
+  { // Repeat rest using multiple zeros version.
     std::vector<RealType> zeros;
-  
     cyl_neumann_zero(static_cast<RealType>(0.0), 3, 1U, std::back_inserter(zeros) );
-
     BOOST_CHECK_CLOSE_FRACTION(zeros[0], static_cast<RealType>(0.89357696627916752158488710205833824122514686193001L), tolerance);
     BOOST_CHECK_CLOSE_FRACTION(zeros[1], static_cast<RealType>(3.9576784193148578683756771869174012814186037655636L), tolerance);
     BOOST_CHECK_CLOSE_FRACTION(zeros[2], static_cast<RealType>(7.0860510603017726976236245968203524689715103811778L), tolerance);
+  }
+  // Order 0: Something always tends to go wrong at zero.
 
+  // TODO ???
 
-     //std::copy(zeros.begin(), zeros.end(),
-     //           std::ostream_iterator<double>(std::cout, "\n"));
-  }
+  /* Order 219/100: This checks a region just below a critical cutoff.
+
+  Table[N[BesselKZero[219/100, n], 50], {n, 1, 20, 4}]
+
+1 | 3.6039149425338727979151181355741147312162055042157
+5 | 16.655399111666833825247894251535326778980614938275
+9 | 29.280564448169163756478439692311605757712873534942
+13 | 41.870269811145814760551599481942750124112093564643
+17 | 54.449180021209532654553613813754733514317929678038
+  */
+  BOOST_CHECK_CLOSE_FRACTION(cyl_neumann_zero(static_cast<RealType>(219)/100, 1U), static_cast<RealType>(3.6039149425338727979151181355741147312162055042157L), tolerance);
+  BOOST_CHECK_CLOSE_FRACTION(cyl_neumann_zero(static_cast<RealType>(219)/100, 5U), static_cast<RealType>(16.655399111666833825247894251535326778980614938275L), tolerance);
+  BOOST_CHECK_CLOSE_FRACTION(cyl_neumann_zero(static_cast<RealType>(219)/100, 17U), static_cast<RealType>(54.449180021209532654553613813754733514317929678038L), tolerance);
+
+  /* Order 221/100: This checks a region just above a critical cutoff.
+
+  Table[N[BesselKZero[220/100, n], 50], {n, 1, 20, 5}]
+
+  1 | 3.6154383428745996706772556069431792744372398748425
+  6 | 19.833435100254138641131431268153987585842088078470
+  11 | 35.592602956438811360473753622212346081080817891225
+  16 | 51.320322762482062633162699745957897178885350674038
+  */
+
+  //BOOST_CHECK_CLOSE_FRACTION(cyl_neumann_zero(static_cast<RealType>(220)/100, 1U), static_cast<RealType>(3.6154383428745996706772556069431792744372398748425L), tolerance);
+  // cyl_neumann_zero(static_cast<RealType>(220)/100, 1U){3.6154383428746009} 
+  //and                             static_cast<RealType>(3.6154383428745996706772556069431792744372398748425L)
+  //                                                     {3.6154383428745995}
+  BOOST_CHECK_CLOSE_FRACTION(cyl_neumann_zero(static_cast<RealType>(220)/100, 6U), static_cast<RealType>(19.833435100254138641131431268153987585842088078470L), tolerance);
+  BOOST_CHECK_CLOSE_FRACTION(cyl_neumann_zero(static_cast<RealType>(220)/100, 11U), static_cast<RealType>(35.592602956438811360473753622212346081080817891225L), tolerance);
+  BOOST_CHECK_CLOSE_FRACTION(cyl_neumann_zero(static_cast<RealType>(220)/100, 16U), static_cast<RealType>(51.320322762482062633162699745957897178885350674038L), tolerance);
+
+  /*  Order 1/1000: A small order.
+
+  Table[N[BesselKZero[1/1000, n], 50], {n, 1, 20, 5}]
+
+  1 | 0.89502371604431360670577815537297733265776195646969
+  6 | 16.502492490954716850993456703662137628148182892787
+  11 | 32.206774708309182755790609144739319753463907110990
+  16 | 47.913467031941494147962476920863688176374357572509
+  */
+
+  BOOST_CHECK_CLOSE_FRACTION(cyl_neumann_zero(static_cast<RealType>(1)/1000, 1U), static_cast<RealType>(0.89502371604431360670577815537297733265776195646969L), tolerance);
+  BOOST_CHECK_CLOSE_FRACTION(cyl_neumann_zero(static_cast<RealType>(1)/1000, 6U), static_cast<RealType>(16.5024924909547168509934567036621376281481828927870L), tolerance);
+  BOOST_CHECK_CLOSE_FRACTION(cyl_neumann_zero(static_cast<RealType>(1)/1000, 11U), static_cast<RealType>(32.206774708309182755790609144739319753463907110990L), tolerance);
+  BOOST_CHECK_CLOSE_FRACTION(cyl_neumann_zero(static_cast<RealType>(1)/1000, 16U), static_cast<RealType>(47.913467031941494147962476920863688176374357572509L), tolerance);
+
+  /* Order 71/19: Merely an intermediate order.
+
+  Table[N[BesselKZero[71/19, n], 50], {n, 1, 20, 5}]
+
+  1 | 5.3527167881149432911848659069476821793319749146616
+  6 | 22.051823727778538215953091664153117627848857279151
+  11 | 37.890091170552491176745048499809370107665221628364
+  16 | 53.651270581421816017744203789836444968181687858095
+  */
+    BOOST_CHECK_CLOSE_FRACTION(cyl_neumann_zero(static_cast<RealType>(71)/19, 1U), static_cast<RealType>(5.3527167881149432911848659069476821793319749146616L), tolerance);
+    BOOST_CHECK_CLOSE_FRACTION(cyl_neumann_zero(static_cast<RealType>(71)/19, 6U), static_cast<RealType>(22.051823727778538215953091664153117627848857279151L), tolerance);
+    BOOST_CHECK_CLOSE_FRACTION(cyl_neumann_zero(static_cast<RealType>(71)/19, 11U), static_cast<RealType>(37.890091170552491176745048499809370107665221628364L), tolerance);
+    BOOST_CHECK_CLOSE_FRACTION(cyl_neumann_zero(static_cast<RealType>(71)/19, 16U), static_cast<RealType>(53.651270581421816017744203789836444968181687858095L), tolerance);
+
+
+  /* Order 7001/19: A medium-large order, small enough to retain moderate efficiency of calculation.
+
+  Table[N[BesselKZero[7001/19, n], 50], {n, 1}]
+
+   1 | 375.18866334770357669101711932706658671250621098115
+
+  Table[N[BesselKZero[7001/19, n], 50], {n, 2}]
+  Standard computation time exceeded :-(
+  */
+  BOOST_CHECK_CLOSE_FRACTION(cyl_neumann_zero(static_cast<RealType>(7001)/19, 1U), static_cast<RealType>(375.18866334770357669101711932706658671250621098115L), tolerance);
+
+ /* A high zero such as the 1000th zero for a modest order such as 71/19.
+
+   Table[N[BesselKZero[71/19, n], 50], {n, 1000}]
+   Standard computation time exceeded :-(
+ */
 
 // Tests of Airy zeros.