ceph/src/boost/libs/math/example/barycentric_interpolation_example.cpp

   1
   2 // Copyright Nick Thompson, 2017
   3
   4 // Distributed under the Boost Software License, Version 1.0.
   5 // (See accompanying file LICENSE_1_0.txt or
   6 // copy at http://www.boost.org/LICENSE_1_0.txt).
   7
   8 #include <iostream>
   9 #include <limits>
  10 #include <vector>
  11
  12 //[barycentric_rational_example
  13
  14 /*`
  15 This example shows how to use barycentric rational interpolation, using Walter Kohn's classic paper
  16 "Solution of the Schrodinger Equation in Periodic Lattices with an Application to Metallic Lithium"
  17 In this paper, Kohn needs to repeatedly solve an ODE (the radial Schrodinger equation) given a potential
  18 which is only known at non-equally samples data.
  19
  20 If he'd only had the barycentric rational interpolant of boost::math!
  21
  22 References: Kohn, W., and N. Rostoker. "Solution of the Schrodinger equation in periodic lattices with an application to metallic lithium." Physical Review 94.5 (1954): 1111.
  23 */
  24
  25 #include <boost/math/interpolators/barycentric_rational.hpp>
  26
  27 int main()
  28 {
  29     // The lithium potential is given in Kohn's paper, Table I:
  30     std::vector<double> r(45);
  31     std::vector<double> mrV(45);
  32
  33     // We'll skip the code for filling the above vectors with data for now...
  34     //<-
  35
  36     r[0] = 0.02; mrV[0] = 5.727;
  37     r[1] = 0.04, mrV[1] = 5.544;
  38     r[2] = 0.06, mrV[2] = 5.450;
  39     r[3] = 0.08, mrV[3] = 5.351;
  40     r[4] = 0.10, mrV[4] = 5.253;
  41     r[5] = 0.12, mrV[5] = 5.157;
  42     r[6] = 0.14, mrV[6] = 5.058;
  43     r[7] = 0.16, mrV[7] = 4.960;
  44     r[8] = 0.18, mrV[8] = 4.862;
  45     r[9] = 0.20, mrV[9] = 4.762;
  46     r[10] = 0.24, mrV[10] = 4.563;
  47     r[11] = 0.28, mrV[11] = 4.360;
  48     r[12] = 0.32, mrV[12] = 4.1584;
  49     r[13] = 0.36, mrV[13] = 3.9463;
  50     r[14] = 0.40, mrV[14] = 3.7360;
  51     r[15] = 0.44, mrV[15] = 3.5429;
  52     r[16] = 0.48, mrV[16] = 3.3797;
  53     r[17] = 0.52, mrV[17] = 3.2417;
  54     r[18] = 0.56, mrV[18] = 3.1209;
  55     r[19] = 0.60, mrV[19] = 3.0138;
  56     r[20] = 0.68, mrV[20] = 2.8342;
  57     r[21] = 0.76, mrV[21] = 2.6881;
  58     r[22] = 0.84, mrV[22] = 2.5662;
  59     r[23] = 0.92, mrV[23] = 2.4242;
  60     r[24] = 1.00, mrV[24] = 2.3766;
  61     r[25] = 1.08, mrV[25] = 2.3058;
  62     r[26] = 1.16, mrV[26] = 2.2458;
  63     r[27] = 1.24, mrV[27] = 2.2035;
  64     r[28] = 1.32, mrV[28] = 2.1661;
  65     r[29] = 1.40, mrV[29] = 2.1350;
  66     r[30] = 1.48, mrV[30] = 2.1090;
  67     r[31] = 1.64, mrV[31] = 2.0697;
  68     r[32] = 1.80, mrV[32] = 2.0466;
  69     r[33] = 1.96, mrV[33] = 2.0325;
  70     r[34] = 2.12, mrV[34] = 2.0288;
  71     r[35] = 2.28, mrV[35] = 2.0292;
  72     r[36] = 2.44, mrV[36] = 2.0228;
  73     r[37] = 2.60, mrV[37] = 2.0124;
  74     r[38] = 2.76, mrV[38] = 2.0065;
  75     r[39] = 2.92, mrV[39] = 2.0031;
  76     r[40] = 3.08, mrV[40] = 2.0015;
  77     r[41] = 3.24, mrV[41] = 2.0008;
  78     r[42] = 3.40, mrV[42] = 2.0004;
  79     r[43] = 3.56, mrV[43] = 2.0002;
  80     r[44] = 3.72, mrV[44] = 2.0001;
  81     //->
  82
  83     // Now we want to interpolate this potential at any r:
  84     boost::math::barycentric_rational<double> b(r.data(), mrV.data(), r.size());
  85
  86     for (size_t i = 1; i < 8; ++i)
  87     {
  88         double r = i*0.5;
  89         std::cout <<  "(r, V) = (" << r << ", " << -b(r)/r << ")\n";
  90     }
  91 }
  92 //] [/barycentric_rational_example]