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

   1 // Copyright Paul A. Bristow 2017
   2 // Copyright John Z. Maddock 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 /*! \brief Graph showing use of Lambert W function.
   9
  10 \details
  11
  12 Both Lambert W0 and W-1 branches can be shown on one graph.
  13 But useful to have another graph for larger values of argument z.
  14 Need two separate graphs for Lambert W0 and -1 prime because
  15 the sensible ranges and axes are too different.
  16
  17 One would get too small LambertW0 in top right and W-1 in bottom left.
  18
  19 */
  20 #ifndef BOOST_MATH_STANDALONE
  21
  22 #include <boost/math/special_functions/lambert_w.hpp>
  23 using boost::math::lambert_w0;
  24 using boost::math::lambert_wm1;
  25 using boost::math::lambert_w0_prime;
  26 using boost::math::lambert_wm1_prime;
  27
  28 #include <boost/math/special_functions.hpp>
  29 using boost::math::isfinite;
  30 #include <boost/svg_plot/svg_2d_plot.hpp>
  31 using namespace boost::svg;
  32 #include <boost/svg_plot/show_2d_settings.hpp>
  33 using boost::svg::show_2d_plot_settings;
  34
  35 #include <iostream>
  36 // using std::cout;
  37 // using std::endl;
  38 #include <exception>
  39 #include <stdexcept>
  40 #include <string>
  41 #include <array>
  42 #include <vector>
  43 #include <utility>
  44 using std::pair;
  45 #include <map>
  46 using std::map;
  47 #include <set>
  48 using std::multiset;
  49 #include <limits>
  50 using std::numeric_limits;
  51 #include <cmath> //
  52
  53   /*!
  54   */
  55 int main()
  56 {
  57   try
  58   {
  59     std::cout << "Lambert W graph example." << std::endl;
  60
  61 //[lambert_w_graph_1
  62 //] [/lambert_w_graph_1]
  63     {
  64       std::map<const double, double> wm1s;   // Lambert W-1 branch values.
  65       std::map<const double, double> w0s;   // Lambert W0 branch values.
  66
  67       std::cout.precision(std::numeric_limits<double>::max_digits10);
  68
  69       int count = 0;
  70       for (double z = -0.36787944117144232159552377016146086744581113103176804; z < 2.8; z += 0.001)
  71       {
  72         double w0 = lambert_w0(z);
  73         w0s[z] = w0;
  74    //     std::cout << "z " << z << ", w = " << w0 << std::endl;
  75         count++;
  76       }
  77       std::cout << "points " << count << std::endl;
  78
  79       count = 0;
  80       for (double z = -0.3678794411714423215955237701614608727; z < -0.001; z += 0.001)
  81       {
  82         double wm1 = lambert_wm1(z);
  83         wm1s[z] = wm1;
  84         count++;
  85       }
  86       std::cout << "points " << count << std::endl;
  87
  88       svg_2d_plot data_plot;
  89       data_plot.title("Lambert W function.")
  90         .x_size(400)
  91         .y_size(300)
  92         .legend_on(true)
  93         .legend_lines(true)
  94         .x_label("z")
  95         .y_label("W")
  96         .x_range(-1, 3.)
  97         .y_range(-4., +1.)
  98         .x_major_interval(1.)
  99         .y_major_interval(1.)
 100         .x_major_grid_on(true)
 101         .y_major_grid_on(true)
 102         //.x_values_on(true)
 103         //.y_values_on(true)
 104         .y_values_rotation(horizontal)
 105         //.plot_window_on(true)
 106         .x_values_precision(3)
 107         .y_values_precision(3)
 108         .coord_precision(4) // Needed to avoid stepping on curves.
 109         .copyright_holder("Paul A. Bristow")
 110         .copyright_date("2018")
 111         //.background_border_color(black);
 112         ;
 113       data_plot.plot(w0s, "W0 branch").line_color(red).shape(none).line_on(true).bezier_on(false).line_width(1);
 114       data_plot.plot(wm1s, "W-1 branch").line_color(blue).shape(none).line_on(true).bezier_on(false).line_width(1);
 115       data_plot.write("./lambert_w_graph");
 116
 117       show_2d_plot_settings(data_plot); // For plot diagnosis only.
 118
 119     } // small z Lambert W
 120
 121     {  // bigger argument z Lambert W
 122
 123       std::map<const double, double> w0s_big;   // Lambert W0 branch values for large z and W.
 124       std::map<const double, double> wm1s_big;   // Lambert W-1 branch values for small z and large -W.
 125       int count = 0;
 126       for (double z = -0.3678794411714423215955237701614608727; z < 10000.; z += 50.)
 127       {
 128         double w0 = lambert_w0(z);
 129         w0s_big[z] = w0;
 130         count++;
 131       }
 132       std::cout << "points " << count << std::endl;
 133
 134       count = 0;
 135       for (double z = -0.3678794411714423215955237701614608727; z < -0.001; z += 0.001)
 136       {
 137         double wm1 = lambert_wm1(z);
 138         wm1s_big[z] = wm1;
 139         count++;
 140       }
 141      std::cout << "Lambert W0 large z argument points = " << count << std::endl;
 142
 143      svg_2d_plot data_plot2;
 144      data_plot2.title("Lambert W0 function for larger z.")
 145       .x_size(400)
 146       .y_size(300)
 147       .legend_on(false)
 148       .x_label("z")
 149       .y_label("W")
 150       //.x_label_on(true)
 151       //.y_label_on(true)
 152       //.xy_values_on(false)
 153       .x_range(-1, 10000.)
 154       .y_range(-1., +8.)
 155       .x_major_interval(2000.)
 156       .y_major_interval(1.)
 157       .x_major_grid_on(true)
 158       .y_major_grid_on(true)
 159       //.x_values_on(true)
 160       //.y_values_on(true)
 161       .y_values_rotation(horizontal)
 162       //.plot_window_on(true)
 163       .x_values_precision(3)
 164       .y_values_precision(3)
 165       .coord_precision(4) // Needed to avoid stepping on curves.
 166       .copyright_holder("Paul A. Bristow")
 167       .copyright_date("2018")
 168       //.background_border_color(black);
 169     ;
 170
 171     data_plot2.plot(w0s_big, "W0 branch").line_color(red).shape(none).line_on(true).bezier_on(false).line_width(1);
 172     // data_plot2.plot(wm1s_big, "W-1 branch").line_color(blue).shape(none).line_on(true).bezier_on(false).line_width(1);
 173     // This wouldn't show anything useful.
 174     data_plot2.write("./lambert_w_graph_big_w");
 175    } // Big argument z Lambert W
 176
 177     { //  Lambert W0 Derivative plots
 178
 179     //  std::map<const double, double> wm1ps;   // Lambert W-1 prime branch values.
 180       std::map<const double, double> w0ps;   // Lambert W0 prime branch values.
 181
 182       std::cout.precision(std::numeric_limits<double>::max_digits10);
 183
 184       int count = 0;
 185       for (double z = -0.36; z < 3.; z += 0.001)
 186       {
 187         double w0p = lambert_w0_prime(z);
 188         w0ps[z] = w0p;
 189         // std::cout << "z " << z << ", w0 = " << w0 << std::endl;
 190         count++;
 191       }
 192       std::cout << "points " << count << std::endl;
 193
 194       //count = 0;
 195       //for (double z = -0.36; z < -0.1; z += 0.001)
 196       //{
 197       //  double wm1p = lambert_wm1_prime(z);
 198       //  std::cout << "z " << z << ", w-1 = " << wm1p << std::endl;
 199       //  wm1ps[z] = wm1p;
 200       //  count++;
 201       //}
 202       //std::cout << "points " << count << std::endl;
 203
 204       svg_2d_plot data_plotp;
 205       data_plotp.title("Lambert W0 prime function.")
 206         .x_size(400)
 207         .y_size(300)
 208         .legend_on(false)
 209         .x_label("z")
 210         .y_label("W0'")
 211         .x_range(-0.3, +1.)
 212         .y_range(0., +5.)
 213         .x_major_interval(0.2)
 214         .y_major_interval(2.)
 215         .x_major_grid_on(true)
 216         .y_major_grid_on(true)
 217         .y_values_rotation(horizontal)
 218         .x_values_precision(3)
 219         .y_values_precision(3)
 220         .coord_precision(4) // Needed to avoid stepping on curves.
 221         .copyright_holder("Paul A. Bristow")
 222         .copyright_date("2018")
 223         ;
 224
 225       // derivative of N[productlog(0, x), 55]  at x=0 to 10
 226       // Plot[D[N[ProductLog[0, x], 55], x], {x, 0, 10}]
 227       // Plot[ProductLog[x]/(x + x ProductLog[x]), {x, 0, 10}]
 228       data_plotp.plot(w0ps, "W0 prime branch").line_color(red).shape(none).line_on(true).bezier_on(false).line_width(1);
 229       data_plotp.write("./lambert_w0_prime_graph");
 230   } // Lambert W0 Derivative plots
 231
 232     { //  Lambert Wm1 Derivative plots
 233
 234     std::map<const double, double> wm1ps;   // Lambert W-1 prime branch values.
 235
 236     std::cout.precision(std::numeric_limits<double>::max_digits10);
 237
 238     int count = 0;
 239     for (double z = -0.3678; z < -0.00001; z += 0.001)
 240     {
 241       double wm1p = lambert_wm1_prime(z);
 242       // std::cout << "z " << z << ", w-1 = " << wm1p << std::endl;
 243       wm1ps[z] = wm1p;
 244       count++;
 245     }
 246     std::cout << "Lambert W-1 prime points = " << count << std::endl;
 247
 248     svg_2d_plot data_plotp;
 249     data_plotp.title("Lambert W-1 prime function.")
 250       .x_size(400)
 251       .y_size(300)
 252       .legend_on(false)
 253       .x_label("z")
 254       .y_label("W-1'")
 255       .x_range(-0.4, +0.01)
 256       .x_major_interval(0.1)
 257       .y_range(-20., -5.)
 258       .y_major_interval(5.)
 259       .x_major_grid_on(true)
 260       .y_major_grid_on(true)
 261       .y_values_rotation(horizontal)
 262       .x_values_precision(3)
 263       .y_values_precision(3)
 264       .coord_precision(4) // Needed to avoid stepping on curves.
 265       .copyright_holder("Paul A. Bristow")
 266       .copyright_date("2018")
 267       ;
 268
 269       // derivative of N[productlog(0, x), 55]  at x=0 to 10
 270       // Plot[D[N[ProductLog[0, x], 55], x], {x, 0, 10}]
 271       // Plot[ProductLog[x]/(x + x ProductLog[x]), {x, 0, 10}]
 272       data_plotp.plot(wm1ps, "W-1 prime branch").line_color(blue).shape(none).line_on(true).bezier_on(false).line_width(1);
 273       data_plotp.write("./lambert_wm1_prime_graph");
 274     } // Lambert W-1 prime graph
 275  } // try
 276   catch (std::exception& ex)
 277   {
 278     std::cout << ex.what() << std::endl;
 279   }
 280 }  // int main()
 281
 282    /*
 283
 284    //[lambert_w_graph_1_output
 285
 286    //] [/lambert_w_graph_1_output]
 287    */
 288
 289 #endif // BOOST_MATH_STANDALONE