1 // (C) Copyright John Maddock 2008.
2 // Use, modification and distribution are subject to the
3 // Boost Software License, Version 1.0. (See accompanying file
4 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
5
6 #ifdef _MSC_VER
7 # pragma warning (disable : 4127) // conditional expression is constant
8 # pragma warning (disable : 4180) // qualifier applied to function type has no meaning; ignored
9 # pragma warning (disable : 4503) // decorated name length exceeded, name was truncated
10 # pragma warning (disable : 4512) // assignment operator could not be generated
11 # pragma warning (disable : 4224) // nonstandard extension used : formal parameter 'function_ptr' was previously defined as a type
12 #endif
13
14 #include <boost/math/special_functions.hpp>
15 #include <boost/math/tools/roots.hpp>
16 #include <boost/function.hpp>
17 #include <boost/bind.hpp>
18
19 #include <list>
20 #include <map>
21 #include <string>
22 #include <boost/svg_plot/svg_2d_plot.hpp>
23 #include <boost/svg_plot/show_2d_settings.hpp>
24
25 class function_arity1_plotter
26 {
27 public:
function_arity1_plotter()28 function_arity1_plotter() : m_min_x(0), m_max_x(0), m_min_y(0), m_max_y(0), m_has_legend(false) {}
29
add(boost::function<double (double)> f,double a,double b,const std::string & name)30 void add(boost::function<double(double)> f, double a, double b, const std::string& name)
31 {
32 if(name.size())
33 m_has_legend = true;
34 //
35 // Now set our x-axis limits:
36 //
37 if(m_max_x == m_min_x)
38 {
39 m_max_x = b;
40 m_min_x = a;
41 }
42 else
43 {
44 if(a < m_min_x)
45 m_min_x = a;
46 if(b > m_max_x)
47 m_max_x = b;
48 }
49 m_points.push_back(std::pair<std::string, std::map<double,double> >(name, std::map<double,double>()));
50 std::map<double,double>& points = m_points.rbegin()->second;
51 double interval = (b - a) / 200;
52 for(double x = a; x <= b; x += interval)
53 {
54 //
55 // Evaluate the function, set the Y axis limits
56 // if needed and then store the pair of points:
57 //
58 double y = f(x);
59 if((m_min_y == m_max_y) && (m_min_y == 0))
60 m_min_y = m_max_y = y;
61 if(m_min_y > y)
62 m_min_y = y;
63 if(m_max_y < y)
64 m_max_y = y;
65 points[x] = y;
66 }
67 }
68
add(const std::map<double,double> & m,const std::string & name)69 void add(const std::map<double, double>& m, const std::string& name)
70 {
71 if(name.size())
72 m_has_legend = true;
73 m_points.push_back(std::pair<std::string, std::map<double,double> >(name, m));
74 std::map<double, double>::const_iterator i = m.begin();
75
76 while(i != m.end())
77 {
78 if((m_min_x == m_min_y) && (m_min_y == 0))
79 {
80 m_min_x = m_max_x = i->first;
81 }
82 if(i->first < m_min_x)
83 {
84 m_min_x = i->first;
85 }
86 if(i->first > m_max_x)
87 {
88 m_max_x = i->first;
89 }
90
91 if((m_min_y == m_max_y) && (m_min_y == 0))
92 {
93 m_min_y = m_max_y = i->second;
94 }
95 if(i->second < m_min_y)
96 {
97 m_min_y = i->second;
98 }
99 if(i->second > m_max_y)
100 {
101 m_max_y = i->second;
102 }
103
104 ++i;
105 }
106 }
107
plot(const std::string & title,const std::string & file,const std::string & x_lable=std::string (),const std::string & y_lable=std::string ())108 void plot(const std::string& title, const std::string& file,
109 const std::string& x_lable = std::string(),
110 const std::string& y_lable = std::string())
111 {
112 using namespace boost::svg;
113
114 static const svg_color colors[5] =
115 {
116 darkblue,
117 darkred,
118 darkgreen,
119 darkorange,
120 chartreuse
121 };
122
123 svg_2d_plot plot;
124 plot.image_x_size(600);
125 plot.image_y_size(400);
126 plot.copyright_holder("John Maddock").copyright_date("2008").boost_license_on(true);
127 plot.coord_precision(4); // Could be 3 for smaller plots?
128 plot.title(title).title_font_size(20).title_on(true);
129 plot.legend_on(m_has_legend);
130
131 double x_delta = (m_max_x - m_min_x) / 50;
132 double y_delta = (m_max_y - m_min_y) / 50;
133 plot.x_range(m_min_x, m_max_x + x_delta)
134 .y_range(m_min_y, m_max_y + y_delta);
135 plot.x_label_on(true).x_label(x_lable);
136 plot.y_label_on(true).y_label(y_lable);
137 plot.y_major_grid_on(false).x_major_grid_on(false);
138 plot.x_num_minor_ticks(3);
139 plot.y_num_minor_ticks(3);
140 //
141 // Work out axis tick intervals:
142 //
143 double l = std::floor(std::log10((m_max_x - m_min_x) / 10) + 0.5);
144 double interval = std::pow(10.0, (int)l);
145 if(((m_max_x - m_min_x) / interval) > 10)
146 interval *= 5;
147 plot.x_major_interval(interval);
148 l = std::floor(std::log10((m_max_y - m_min_y) / 10) + 0.5);
149 interval = std::pow(10.0, (int)l);
150 if(((m_max_y - m_min_y) / interval) > 10)
151 interval *= 5;
152 plot.y_major_interval(interval);
153 plot.plot_window_on(true);
154 plot.plot_border_color(lightslategray)
155 .background_border_color(lightslategray)
156 .legend_border_color(lightslategray)
157 .legend_background_color(white);
158
159 int color_index = 0;
160
161 for(std::list<std::pair<std::string, std::map<double,double> > >::const_iterator i = m_points.begin();
162 i != m_points.end(); ++i)
163 {
164 plot.plot(i->second, i->first)
165 .line_on(true)
166 .line_color(colors[color_index])
167 .line_width(1.)
168 .shape(none);
169 if(i->first.size())
170 ++color_index;
171 color_index = color_index % (sizeof(colors)/sizeof(colors[0]));
172 }
173 plot.write(file);
174 }
175
clear()176 void clear()
177 {
178 m_points.clear();
179 m_min_x = m_min_y = m_max_x = m_max_y = 0;
180 m_has_legend = false;
181 }
182
183 private:
184 std::list<std::pair<std::string, std::map<double, double> > > m_points;
185 double m_min_x, m_max_x, m_min_y, m_max_y;
186 bool m_has_legend;
187 };
188
189 template <class F>
190 struct location_finder
191 {
location_finderlocation_finder192 location_finder(F _f, double t, double x0) : f(_f), target(t), x_off(x0){}
193
operator ()location_finder194 double operator()(double x)
195 {
196 try
197 {
198 return f(x + x_off) - target;
199 }
200 catch(const std::overflow_error&)
201 {
202 return boost::math::tools::max_value<double>();
203 }
204 catch(const std::domain_error&)
205 {
206 if(x + x_off == x_off)
207 return f(x_off + boost::math::tools::epsilon<double>() * x_off);
208 throw;
209 }
210 }
211
212 private:
213 F f;
214 double target;
215 double x_off;
216 };
217
218 template <class F>
find_end_point(F f,double x0,double target,bool rising,double x_off=0)219 double find_end_point(F f, double x0, double target, bool rising, double x_off = 0)
220 {
221 boost::math::tools::eps_tolerance<double> tol(50);
222 boost::uintmax_t max_iter = 1000;
223 return x_off + boost::math::tools::bracket_and_solve_root(
224 location_finder<F>(f, target, x_off),
225 x0,
226 1.5,
227 rising,
228 tol,
229 max_iter).first;
230 }
231
sqrt1pm1(double x)232 double sqrt1pm1(double x)
233 {
234 return boost::math::sqrt1pm1(x);
235 }
236
lbeta(double a,double b)237 double lbeta(double a, double b)
238 {
239 return std::log(boost::math::beta(a, b));
240 }
241
main()242 int main()
243 {
244 function_arity1_plotter plot;
245 double (*f)(double);
246 double (*f2)(double, double);
247 double (*f2u)(unsigned, double);
248 double (*f2i)(int, double);
249 double (*f3)(double, double, double);
250 double (*f4)(double, double, double, double);
251 double max_val;
252
253 f = boost::math::zeta;
254 plot.add(f, find_end_point(f, 0.1, 40.0, false, 1.0), 10, "");
255 plot.add(f, -20, find_end_point(f, -0.1, -40.0, false, 1.0), "");
256 plot.plot("Zeta Function Over [-20,10]", "zeta1.svg", "z", "zeta(z)");
257
258 plot.clear();
259 plot.add(f, -14, 0, "");
260 plot.plot("Zeta Function Over [-14,0]", "zeta2.svg", "z", "zeta(z)");
261
262 f = boost::math::tgamma;
263 max_val = f(6);
264 plot.clear();
265 plot.add(f, find_end_point(f, 0.1, max_val, false), 6, "");
266 plot.add(f, find_end_point(f, 0.1, -max_val, true, -1), find_end_point(f, -0.1, -max_val, false), "");
267 plot.add(f, find_end_point(f, 0.1, max_val, false, -2), find_end_point(f, -0.1, max_val, true, -1), "");
268 plot.add(f, find_end_point(f, 0.1, -max_val, true, -3), find_end_point(f, -0.1, -max_val, false, -2), "");
269 plot.add(f, find_end_point(f, 0.1, max_val, false, -4), find_end_point(f, -0.1, max_val, true, -3), "");
270 plot.plot("tgamma", "tgamma.svg", "z", "tgamma(z)");
271
272 f = boost::math::lgamma;
273 max_val = f(10);
274 plot.clear();
275 plot.add(f, find_end_point(f, 0.1, max_val, false), 10, "");
276 plot.add(f, find_end_point(f, 0.1, max_val, false, -1), find_end_point(f, -0.1, max_val, true), "");
277 plot.add(f, find_end_point(f, 0.1, max_val, false, -2), find_end_point(f, -0.1, max_val, true, -1), "");
278 plot.add(f, find_end_point(f, 0.1, max_val, false, -3), find_end_point(f, -0.1, max_val, true, -2), "");
279 plot.add(f, find_end_point(f, 0.1, max_val, false, -4), find_end_point(f, -0.1, max_val, true, -3), "");
280 plot.add(f, find_end_point(f, 0.1, max_val, false, -5), find_end_point(f, -0.1, max_val, true, -4), "");
281 plot.plot("lgamma", "lgamma.svg", "z", "lgamma(z)");
282
283 f = boost::math::digamma;
284 max_val = 10;
285 plot.clear();
286 plot.add(f, find_end_point(f, 0.1, -max_val, true), 10, "");
287 plot.add(f, find_end_point(f, 0.1, -max_val, true, -1), find_end_point(f, -0.1, max_val, true), "");
288 plot.add(f, find_end_point(f, 0.1, -max_val, true, -2), find_end_point(f, -0.1, max_val, true, -1), "");
289 plot.add(f, find_end_point(f, 0.1, -max_val, true, -3), find_end_point(f, -0.1, max_val, true, -2), "");
290 plot.add(f, find_end_point(f, 0.1, -max_val, true, -4), find_end_point(f, -0.1, max_val, true, -3), "");
291 plot.plot("digamma", "digamma.svg", "z", "digamma(z)");
292
293 f = boost::math::erf;
294 plot.clear();
295 plot.add(f, -3, 3, "");
296 plot.plot("erf", "erf.svg", "z", "erf(z)");
297 f = boost::math::erfc;
298 plot.clear();
299 plot.add(f, -3, 3, "");
300 plot.plot("erfc", "erfc.svg", "z", "erfc(z)");
301
302 f = boost::math::erf_inv;
303 plot.clear();
304 plot.add(f, find_end_point(f, 0.1, -3, true, -1), find_end_point(f, -0.1, 3, true, 1), "");
305 plot.plot("erf_inv", "erf_inv.svg", "z", "erf_inv(z)");
306 f = boost::math::erfc_inv;
307 plot.clear();
308 plot.add(f, find_end_point(f, 0.1, 3, false), find_end_point(f, -0.1, -3, false, 2), "");
309 plot.plot("erfc_inv", "erfc_inv.svg", "z", "erfc_inv(z)");
310
311 f = boost::math::log1p;
312 plot.clear();
313 plot.add(f, find_end_point(f, 0.1, -10, true, -1), 10, "");
314 plot.plot("log1p", "log1p.svg", "z", "log1p(z)");
315
316 f = boost::math::expm1;
317 plot.clear();
318 plot.add(f, -4, 2, "");
319 plot.plot("expm1", "expm1.svg", "z", "expm1(z)");
320
321 f = boost::math::cbrt;
322 plot.clear();
323 plot.add(f, -10, 10, "");
324 plot.plot("cbrt", "cbrt.svg", "z", "cbrt(z)");
325
326 f = sqrt1pm1;
327 plot.clear();
328 plot.add(f, find_end_point(f, 0.1, -10, true, -1), 5, "");
329 plot.plot("sqrt1pm1", "sqrt1pm1.svg", "z", "sqrt1pm1(z)");
330
331 f2 = boost::math::powm1;
332 plot.clear();
333 plot.add(boost::bind(f2, 0.0001, _1), find_end_point(boost::bind(f2, 0.0001, _1), -1, 10, false), 5, "a=0.0001");
334 plot.add(boost::bind(f2, 0.001, _1), find_end_point(boost::bind(f2, 0.001, _1), -1, 10, false), 5, "a=0.001");
335 plot.add(boost::bind(f2, 0.01, _1), find_end_point(boost::bind(f2, 0.01, _1), -1, 10, false), 5, "a=0.01");
336 plot.add(boost::bind(f2, 0.1, _1), find_end_point(boost::bind(f2, 0.1, _1), -1, 10, false), 5, "a=0.1");
337 plot.add(boost::bind(f2, 0.75, _1), -5, 5, "a=0.75");
338 plot.add(boost::bind(f2, 1.25, _1), -5, 5, "a=1.25");
339 plot.plot("powm1", "powm1.svg", "z", "powm1(a, z)");
340
341 f = boost::math::sinc_pi;
342 plot.clear();
343 plot.add(f, -10, 10, "");
344 plot.plot("sinc_pi", "sinc_pi.svg", "z", "sinc_pi(z)");
345
346 f = boost::math::sinhc_pi;
347 plot.clear();
348 plot.add(f, -5, 5, "");
349 plot.plot("sinhc_pi", "sinhc_pi.svg", "z", "sinhc_pi(z)");
350
351 f = boost::math::acosh;
352 plot.clear();
353 plot.add(f, 1, 10, "");
354 plot.plot("acosh", "acosh.svg", "z", "acosh(z)");
355
356 f = boost::math::asinh;
357 plot.clear();
358 plot.add(f, -10, 10, "");
359 plot.plot("asinh", "asinh.svg", "z", "asinh(z)");
360
361 f = boost::math::atanh;
362 plot.clear();
363 plot.add(f, find_end_point(f, 0.1, -5, true, -1), find_end_point(f, -0.1, 5, true, 1), "");
364 plot.plot("atanh", "atanh.svg", "z", "atanh(z)");
365
366 f2 = boost::math::tgamma_delta_ratio;
367 plot.clear();
368 plot.add(boost::bind(f2, _1, -0.5), 1, 40, "delta = -0.5");
369 plot.add(boost::bind(f2, _1, -0.2), 1, 40, "delta = -0.2");
370 plot.add(boost::bind(f2, _1, -0.1), 1, 40, "delta = -0.1");
371 plot.add(boost::bind(f2, _1, 0.1), 1, 40, "delta = 0.1");
372 plot.add(boost::bind(f2, _1, 0.2), 1, 40, "delta = 0.2");
373 plot.add(boost::bind(f2, _1, 0.5), 1, 40, "delta = 0.5");
374 plot.add(boost::bind(f2, _1, 1.0), 1, 40, "delta = 1.0");
375 plot.plot("tgamma_delta_ratio", "tgamma_delta_ratio.svg", "z", "tgamma_delta_ratio(delta, z)");
376
377 f2 = boost::math::gamma_p;
378 plot.clear();
379 plot.add(boost::bind(f2, 0.5, _1), 0, 20, "a = 0.5");
380 plot.add(boost::bind(f2, 1.0, _1), 0, 20, "a = 1.0");
381 plot.add(boost::bind(f2, 5.0, _1), 0, 20, "a = 5.0");
382 plot.add(boost::bind(f2, 10.0, _1), 0, 20, "a = 10.0");
383 plot.plot("gamma_p", "gamma_p.svg", "z", "gamma_p(a, z)");
384
385 f2 = boost::math::gamma_q;
386 plot.clear();
387 plot.add(boost::bind(f2, 0.5, _1), 0, 20, "a = 0.5");
388 plot.add(boost::bind(f2, 1.0, _1), 0, 20, "a = 1.0");
389 plot.add(boost::bind(f2, 5.0, _1), 0, 20, "a = 5.0");
390 plot.add(boost::bind(f2, 10.0, _1), 0, 20, "a = 10.0");
391 plot.plot("gamma_q", "gamma_q.svg", "z", "gamma_q(a, z)");
392
393 f2 = lbeta;
394 plot.clear();
395 plot.add(boost::bind(f2, 0.5, _1), 0.00001, 5, "a = 0.5");
396 plot.add(boost::bind(f2, 1.0, _1), 0.00001, 5, "a = 1.0");
397 plot.add(boost::bind(f2, 5.0, _1), 0.00001, 5, "a = 5.0");
398 plot.add(boost::bind(f2, 10.0, _1), 0.00001, 5, "a = 10.0");
399 plot.plot("beta", "beta.svg", "z", "log(beta(a, z))");
400
401 f = boost::math::expint;
402 max_val = f(4);
403 plot.clear();
404 plot.add(f, find_end_point(f, 0.1, -max_val, true), 4, "");
405 plot.add(f, -3, find_end_point(f, -0.1, -max_val, false), "");
406 plot.plot("Exponential Integral Ei", "expint_i.svg", "z", "expint(z)");
407
408 f2u = boost::math::expint;
409 max_val = 1;
410 plot.clear();
411 plot.add(boost::bind(f2u, 1, _1), find_end_point(boost::bind(f2u, 1, _1), 0.1, max_val, false), 2, "n = 1 ");
412 plot.add(boost::bind(f2u, 2, _1), find_end_point(boost::bind(f2u, 2, _1), 0.1, max_val, false), 2, "n = 2 ");
413 plot.add(boost::bind(f2u, 3, _1), 0, 2, "n = 3 ");
414 plot.add(boost::bind(f2u, 4, _1), 0, 2, "n = 4 ");
415 plot.plot("Exponential Integral En", "expint2.svg", "z", "expint(n, z)");
416
417 f3 = boost::math::ibeta;
418 plot.clear();
419 plot.add(boost::bind(f3, 9, 1, _1), 0, 1, "a = 9, b = 1");
420 plot.add(boost::bind(f3, 7, 2, _1), 0, 1, "a = 7, b = 2");
421 plot.add(boost::bind(f3, 5, 5, _1), 0, 1, "a = 5, b = 5");
422 plot.add(boost::bind(f3, 2, 7, _1), 0, 1, "a = 2, b = 7");
423 plot.add(boost::bind(f3, 1, 9, _1), 0, 1, "a = 1, b = 9");
424 plot.plot("ibeta", "ibeta.svg", "z", "ibeta(a, b, z)");
425
426 f2i = boost::math::legendre_p;
427 plot.clear();
428 plot.add(boost::bind(f2i, 1, _1), -1, 1, "l = 1");
429 plot.add(boost::bind(f2i, 2, _1), -1, 1, "l = 2");
430 plot.add(boost::bind(f2i, 3, _1), -1, 1, "l = 3");
431 plot.add(boost::bind(f2i, 4, _1), -1, 1, "l = 4");
432 plot.add(boost::bind(f2i, 5, _1), -1, 1, "l = 5");
433 plot.plot("Legendre Polynomials", "legendre_p.svg", "x", "legendre_p(l, x)");
434
435 f2u = boost::math::legendre_q;
436 plot.clear();
437 plot.add(boost::bind(f2u, 1, _1), -0.95, 0.95, "l = 1");
438 plot.add(boost::bind(f2u, 2, _1), -0.95, 0.95, "l = 2");
439 plot.add(boost::bind(f2u, 3, _1), -0.95, 0.95, "l = 3");
440 plot.add(boost::bind(f2u, 4, _1), -0.95, 0.95, "l = 4");
441 plot.add(boost::bind(f2u, 5, _1), -0.95, 0.95, "l = 5");
442 plot.plot("Legendre Polynomials of the Second Kind", "legendre_q.svg", "x", "legendre_q(l, x)");
443
444 f2u = boost::math::laguerre;
445 plot.clear();
446 plot.add(boost::bind(f2u, 0, _1), -5, 10, "n = 0");
447 plot.add(boost::bind(f2u, 1, _1), -5, 10, "n = 1");
448 plot.add(boost::bind(f2u, 2, _1),
449 find_end_point(boost::bind(f2u, 2, _1), -2, 20, false),
450 find_end_point(boost::bind(f2u, 2, _1), 4, 20, true),
451 "n = 2");
452 plot.add(boost::bind(f2u, 3, _1),
453 find_end_point(boost::bind(f2u, 3, _1), -2, 20, false),
454 find_end_point(boost::bind(f2u, 3, _1), 1, 20, false, 8),
455 "n = 3");
456 plot.add(boost::bind(f2u, 4, _1),
457 find_end_point(boost::bind(f2u, 4, _1), -2, 20, false),
458 find_end_point(boost::bind(f2u, 4, _1), 1, 20, true, 8),
459 "n = 4");
460 plot.add(boost::bind(f2u, 5, _1),
461 find_end_point(boost::bind(f2u, 5, _1), -2, 20, false),
462 find_end_point(boost::bind(f2u, 5, _1), 1, 20, true, 8),
463 "n = 5");
464 plot.plot("Laguerre Polynomials", "laguerre.svg", "x", "laguerre(n, x)");
465
466 f2u = boost::math::hermite;
467 plot.clear();
468 plot.add(boost::bind(f2u, 0, _1), -1.8, 1.8, "n = 0");
469 plot.add(boost::bind(f2u, 1, _1), -1.8, 1.8, "n = 1");
470 plot.add(boost::bind(f2u, 2, _1), -1.8, 1.8, "n = 2");
471 plot.add(boost::bind(f2u, 3, _1), -1.8, 1.8, "n = 3");
472 plot.add(boost::bind(f2u, 4, _1), -1.8, 1.8, "n = 4");
473 plot.plot("Hermite Polynomials", "hermite.svg", "x", "hermite(n, x)");
474
475 f2 = boost::math::cyl_bessel_j;
476 plot.clear();
477 plot.add(boost::bind(f2, 0, _1), -20, 20, "v = 0");
478 plot.add(boost::bind(f2, 1, _1), -20, 20, "v = 1");
479 plot.add(boost::bind(f2, 2, _1), -20, 20, "v = 2");
480 plot.add(boost::bind(f2, 3, _1), -20, 20, "v = 3");
481 plot.add(boost::bind(f2, 4, _1), -20, 20, "v = 4");
482 plot.plot("Bessel J", "cyl_bessel_j.svg", "x", "cyl_bessel_j(v, x)");
483
484 f2 = boost::math::cyl_neumann;
485 plot.clear();
486 plot.add(boost::bind(f2, 0, _1), find_end_point(boost::bind(f2, 0, _1), 0.1, -5, true), 20, "v = 0");
487 plot.add(boost::bind(f2, 1, _1), find_end_point(boost::bind(f2, 1, _1), 0.1, -5, true), 20, "v = 1");
488 plot.add(boost::bind(f2, 2, _1), find_end_point(boost::bind(f2, 2, _1), 0.1, -5, true), 20, "v = 2");
489 plot.add(boost::bind(f2, 3, _1), find_end_point(boost::bind(f2, 3, _1), 0.1, -5, true), 20, "v = 3");
490 plot.add(boost::bind(f2, 4, _1), find_end_point(boost::bind(f2, 4, _1), 0.1, -5, true), 20, "v = 4");
491 plot.plot("Bessel Y", "cyl_neumann.svg", "x", "cyl_neumann(v, x)");
492
493 f2 = boost::math::cyl_bessel_i;
494 plot.clear();
495 plot.add(boost::bind(f2, 0, _1), find_end_point(boost::bind(f2, 0, _1), -0.1, 20, false), find_end_point(boost::bind(f2, 0, _1), 0.1, 20, true), "v = 0");
496 plot.add(boost::bind(f2, 2, _1), find_end_point(boost::bind(f2, 2, _1), -0.1, 20, false), find_end_point(boost::bind(f2, 2, _1), 0.1, 20, true), "v = 2");
497 plot.add(boost::bind(f2, 5, _1), find_end_point(boost::bind(f2, 5, _1), -0.1, -20, true), find_end_point(boost::bind(f2, 5, _1), 0.1, 20, true), "v = 5");
498 plot.add(boost::bind(f2, 7, _1), find_end_point(boost::bind(f2, 7, _1), -0.1, -20, true), find_end_point(boost::bind(f2, 7, _1), 0.1, 20, true), "v = 7");
499 plot.add(boost::bind(f2, 10, _1), find_end_point(boost::bind(f2, 10, _1), -0.1, 20, false), find_end_point(boost::bind(f2, 10, _1), 0.1, 20, true), "v = 10");
500 plot.plot("Bessel I", "cyl_bessel_i.svg", "x", "cyl_bessel_i(v, x)");
501
502 f2 = boost::math::cyl_bessel_k;
503 plot.clear();
504 plot.add(boost::bind(f2, 0, _1), find_end_point(boost::bind(f2, 0, _1), 0.1, 10, false), 10, "v = 0");
505 plot.add(boost::bind(f2, 2, _1), find_end_point(boost::bind(f2, 2, _1), 0.1, 10, false), 10, "v = 2");
506 plot.add(boost::bind(f2, 5, _1), find_end_point(boost::bind(f2, 5, _1), 0.1, 10, false), 10, "v = 5");
507 plot.add(boost::bind(f2, 7, _1), find_end_point(boost::bind(f2, 7, _1), 0.1, 10, false), 10, "v = 7");
508 plot.add(boost::bind(f2, 10, _1), find_end_point(boost::bind(f2, 10, _1), 0.1, 10, false), 10, "v = 10");
509 plot.plot("Bessel K", "cyl_bessel_k.svg", "x", "cyl_bessel_k(v, x)");
510
511 f2u = boost::math::sph_bessel;
512 plot.clear();
513 plot.add(boost::bind(f2u, 0, _1), 0, 20, "v = 0");
514 plot.add(boost::bind(f2u, 2, _1), 0, 20, "v = 2");
515 plot.add(boost::bind(f2u, 5, _1), 0, 20, "v = 5");
516 plot.add(boost::bind(f2u, 7, _1), 0, 20, "v = 7");
517 plot.add(boost::bind(f2u, 10, _1), 0, 20, "v = 10");
518 plot.plot("Bessel j", "sph_bessel.svg", "x", "sph_bessel(v, x)");
519
520 f2u = boost::math::sph_neumann;
521 plot.clear();
522 plot.add(boost::bind(f2u, 0, _1), find_end_point(boost::bind(f2u, 0, _1), 0.1, -5, true), 20, "v = 0");
523 plot.add(boost::bind(f2u, 2, _1), find_end_point(boost::bind(f2u, 2, _1), 0.1, -5, true), 20, "v = 2");
524 plot.add(boost::bind(f2u, 5, _1), find_end_point(boost::bind(f2u, 5, _1), 0.1, -5, true), 20, "v = 5");
525 plot.add(boost::bind(f2u, 7, _1), find_end_point(boost::bind(f2u, 7, _1), 0.1, -5, true), 20, "v = 7");
526 plot.add(boost::bind(f2u, 10, _1), find_end_point(boost::bind(f2u, 10, _1), 0.1, -5, true), 20, "v = 10");
527 plot.plot("Bessel y", "sph_neumann.svg", "x", "sph_neumann(v, x)");
528
529 f4 = boost::math::ellint_rj;
530 plot.clear();
531 plot.add(boost::bind(f4, _1, _1, _1, _1), find_end_point(boost::bind(f4, _1, _1, _1, _1), 0.1, 10, false), 4, "RJ");
532 f3 = boost::math::ellint_rf;
533 plot.add(boost::bind(f3, _1, _1, _1), find_end_point(boost::bind(f3, _1, _1, _1), 0.1, 10, false), 4, "RF");
534 plot.plot("Elliptic Integrals", "ellint_carlson.svg", "x", "");
535
536 f2 = boost::math::ellint_1;
537 plot.clear();
538 plot.add(boost::bind(f2, _1, 0.5), -0.9, 0.9, "φ=0.5");
539 plot.add(boost::bind(f2, _1, 0.75), -0.9, 0.9, "φ=0.75");
540 plot.add(boost::bind(f2, _1, 1.25), -0.9, 0.9, "φ=1.25");
541 plot.add(boost::bind(f2, _1, boost::math::constants::pi<double>() / 2), -0.9, 0.9, "φ=π/2");
542 plot.plot("Elliptic Of the First Kind", "ellint_1.svg", "k", "ellint_1(k, phi)");
543
544 f2 = boost::math::ellint_2;
545 plot.clear();
546 plot.add(boost::bind(f2, _1, 0.5), -1, 1, "φ=0.5");
547 plot.add(boost::bind(f2, _1, 0.75), -1, 1, "φ=0.75");
548 plot.add(boost::bind(f2, _1, 1.25), -1, 1, "φ=1.25");
549 plot.add(boost::bind(f2, _1, boost::math::constants::pi<double>() / 2), -1, 1, "φ=π/2");
550 plot.plot("Elliptic Of the Second Kind", "ellint_2.svg", "k", "ellint_2(k, phi)");
551
552 f3 = boost::math::ellint_3;
553 plot.clear();
554 plot.add(boost::bind(f3, _1, 0, 1.25), -1, 1, "n=0 φ=1.25");
555 plot.add(boost::bind(f3, _1, 0.5, 1.25), -1, 1, "n=0.5 φ=1.25");
556 plot.add(boost::bind(f3, _1, 0.25, boost::math::constants::pi<double>() / 2),
557 find_end_point(
558 boost::bind(f3, _1, 0.25, boost::math::constants::pi<double>() / 2),
559 0.5, 4, false, -1),
560 find_end_point(
561 boost::bind(f3, _1, 0.25, boost::math::constants::pi<double>() / 2),
562 -0.5, 4, true, 1), "n=0.25 φ=π/2");
563 plot.add(boost::bind(f3, _1, 0.75, boost::math::constants::pi<double>() / 2),
564 find_end_point(
565 boost::bind(f3, _1, 0.75, boost::math::constants::pi<double>() / 2),
566 0.5, 4, false, -1),
567 find_end_point(
568 boost::bind(f3, _1, 0.75, boost::math::constants::pi<double>() / 2),
569 -0.5, 4, true, 1), "n=0.75 φ=π/2");
570 plot.plot("Elliptic Of the Third Kind", "ellint_3.svg", "k", "ellint_3(k, n, phi)");
571
572 f2 = boost::math::jacobi_sn;
573 plot.clear();
574 plot.add(boost::bind(f2, 0, _1), -10, 10, "k=0");
575 plot.add(boost::bind(f2, 0.5, _1), -10, 10, "k=0.5");
576 plot.add(boost::bind(f2, 0.75, _1), -10, 10, "k=0.75");
577 plot.add(boost::bind(f2, 0.95, _1), -10, 10, "k=0.95");
578 plot.add(boost::bind(f2, 1, _1), -10, 10, "k=1");
579 plot.plot("Jacobi Elliptic sn", "jacobi_sn.svg", "k", "jacobi_sn(k, u)");
580
581 f2 = boost::math::jacobi_cn;
582 plot.clear();
583 plot.add(boost::bind(f2, 0, _1), -10, 10, "k=0");
584 plot.add(boost::bind(f2, 0.5, _1), -10, 10, "k=0.5");
585 plot.add(boost::bind(f2, 0.75, _1), -10, 10, "k=0.75");
586 plot.add(boost::bind(f2, 0.95, _1), -10, 10, "k=0.95");
587 plot.add(boost::bind(f2, 1, _1), -10, 10, "k=1");
588 plot.plot("Jacobi Elliptic cn", "jacobi_cn.svg", "k", "jacobi_cn(k, u)");
589
590 f2 = boost::math::jacobi_dn;
591 plot.clear();
592 plot.add(boost::bind(f2, 0, _1), -10, 10, "k=0");
593 plot.add(boost::bind(f2, 0.5, _1), -10, 10, "k=0.5");
594 plot.add(boost::bind(f2, 0.75, _1), -10, 10, "k=0.75");
595 plot.add(boost::bind(f2, 0.95, _1), -10, 10, "k=0.95");
596 plot.add(boost::bind(f2, 1, _1), -10, 10, "k=1");
597 plot.plot("Jacobi Elliptic dn", "jacobi_dn.svg", "k", "jacobi_dn(k, u)");
598
599 f2 = boost::math::jacobi_cd;
600 plot.clear();
601 plot.add(boost::bind(f2, 0, _1), -10, 10, "k=0");
602 plot.add(boost::bind(f2, 0.5, _1), -10, 10, "k=0.5");
603 plot.add(boost::bind(f2, 0.75, _1), -10, 10, "k=0.75");
604 plot.add(boost::bind(f2, 0.95, _1), -10, 10, "k=0.95");
605 plot.add(boost::bind(f2, 1, _1), -10, 10, "k=1");
606 plot.plot("Jacobi Elliptic cd", "jacobi_cd.svg", "k", "jacobi_cd(k, u)");
607
608 f2 = boost::math::jacobi_cs;
609 plot.clear();
610 plot.add(boost::bind(f2, 0, _1), 0.1, 3, "k=0");
611 plot.add(boost::bind(f2, 0.5, _1), 0.1, 3, "k=0.5");
612 plot.add(boost::bind(f2, 0.75, _1), 0.1, 3, "k=0.75");
613 plot.add(boost::bind(f2, 0.95, _1), 0.1, 3, "k=0.95");
614 plot.add(boost::bind(f2, 1, _1), 0.1, 3, "k=1");
615 plot.plot("Jacobi Elliptic cs", "jacobi_cs.svg", "k", "jacobi_cs(k, u)");
616
617 f2 = boost::math::jacobi_dc;
618 plot.clear();
619 plot.add(boost::bind(f2, 0, _1), -10, 10, "k=0");
620 plot.add(boost::bind(f2, 0.5, _1), -10, 10, "k=0.5");
621 plot.add(boost::bind(f2, 0.75, _1), -10, 10, "k=0.75");
622 plot.add(boost::bind(f2, 0.95, _1), -10, 10, "k=0.95");
623 plot.add(boost::bind(f2, 1, _1), -10, 10, "k=1");
624 plot.plot("Jacobi Elliptic dc", "jacobi_dc.svg", "k", "jacobi_dc(k, u)");
625
626 f2 = boost::math::jacobi_ds;
627 plot.clear();
628 plot.add(boost::bind(f2, 0, _1), 0.1, 3, "k=0");
629 plot.add(boost::bind(f2, 0.5, _1), 0.1, 3, "k=0.5");
630 plot.add(boost::bind(f2, 0.75, _1), 0.1, 3, "k=0.75");
631 plot.add(boost::bind(f2, 0.95, _1), 0.1, 3, "k=0.95");
632 plot.add(boost::bind(f2, 1, _1), 0.1, 3, "k=1");
633 plot.plot("Jacobi Elliptic ds", "jacobi_ds.svg", "k", "jacobi_ds(k, u)");
634
635 f2 = boost::math::jacobi_nc;
636 plot.clear();
637 plot.add(boost::bind(f2, 0, _1), -5, 5, "k=0");
638 plot.add(boost::bind(f2, 0.5, _1), -5, 5, "k=0.5");
639 plot.add(boost::bind(f2, 0.75, _1), -5, 5, "k=0.75");
640 plot.add(boost::bind(f2, 0.95, _1), -5, 5, "k=0.95");
641 plot.add(boost::bind(f2, 1, _1), -5, 5, "k=1");
642 plot.plot("Jacobi Elliptic nc", "jacobi_nc.svg", "k", "jacobi_nc(k, u)");
643
644 f2 = boost::math::jacobi_ns;
645 plot.clear();
646 plot.add(boost::bind(f2, 0, _1), 0.1, 4, "k=0");
647 plot.add(boost::bind(f2, 0.5, _1), 0.1, 4, "k=0.5");
648 plot.add(boost::bind(f2, 0.75, _1), 0.1, 4, "k=0.75");
649 plot.add(boost::bind(f2, 0.95, _1), 0.1, 4, "k=0.95");
650 plot.add(boost::bind(f2, 1, _1), 0.1, 4, "k=1");
651 plot.plot("Jacobi Elliptic ns", "jacobi_ns.svg", "k", "jacobi_ns(k, u)");
652
653 f2 = boost::math::jacobi_nd;
654 plot.clear();
655 plot.add(boost::bind(f2, 0, _1), -2, 2, "k=0");
656 plot.add(boost::bind(f2, 0.5, _1), -2, 2, "k=0.5");
657 plot.add(boost::bind(f2, 0.75, _1), -2, 2, "k=0.75");
658 plot.add(boost::bind(f2, 0.95, _1), -2, 2, "k=0.95");
659 plot.add(boost::bind(f2, 1, _1), -2, 2, "k=1");
660 plot.plot("Jacobi Elliptic nd", "jacobi_nd.svg", "k", "jacobi_nd(k, u)");
661
662 f2 = boost::math::jacobi_sc;
663 plot.clear();
664 plot.add(boost::bind(f2, 0, _1), -5, 5, "k=0");
665 plot.add(boost::bind(f2, 0.5, _1), -5, 5, "k=0.5");
666 plot.add(boost::bind(f2, 0.75, _1), -5, 5, "k=0.75");
667 plot.add(boost::bind(f2, 0.95, _1), -5, 5, "k=0.95");
668 plot.add(boost::bind(f2, 1, _1), -5, 5, "k=1");
669 plot.plot("Jacobi Elliptic sc", "jacobi_sc.svg", "k", "jacobi_sc(k, u)");
670
671 f2 = boost::math::jacobi_sd;
672 plot.clear();
673 plot.add(boost::bind(f2, 0, _1), -2.5, 2.5, "k=0");
674 plot.add(boost::bind(f2, 0.5, _1), -2.5, 2.5, "k=0.5");
675 plot.add(boost::bind(f2, 0.75, _1), -2.5, 2.5, "k=0.75");
676 plot.add(boost::bind(f2, 0.95, _1), -2.5, 2.5, "k=0.95");
677 plot.add(boost::bind(f2, 1, _1), -2.5, 2.5, "k=1");
678 plot.plot("Jacobi Elliptic sd", "jacobi_sd.svg", "k", "jacobi_sd(k, u)");
679
680 f = boost::math::airy_ai;
681 plot.clear();
682 plot.add(f, -20, 20, "");
683 plot.plot("Ai", "airy_ai.svg", "z", "airy_ai(z)");
684
685 f = boost::math::airy_bi;
686 plot.clear();
687 plot.add(f, -20, 3, "");
688 plot.plot("Bi", "airy_bi.svg", "z", "airy_bi(z)");
689
690 f = boost::math::airy_ai_prime;
691 plot.clear();
692 plot.add(f, -20, 20, "");
693 plot.plot("Ai'", "airy_aip.svg", "z", "airy_ai_prime(z)");
694
695 f = boost::math::airy_bi_prime;
696 plot.clear();
697 plot.add(f, -20, 3, "");
698 plot.plot("Bi'", "airy_bip.svg", "z", "airy_bi_prime(z)");
699
700 f = boost::math::trigamma;
701 max_val = 30;
702 plot.clear();
703 plot.add(f, find_end_point(f, 0.1, max_val, false), 5, "");
704 plot.add(f, find_end_point(f, 0.1, max_val, false, -1), find_end_point(f, -0.1, max_val, true), "");
705 plot.add(f, find_end_point(f, 0.1, max_val, false, -2), find_end_point(f, -0.1, max_val, true, -1), "");
706 plot.add(f, find_end_point(f, 0.1, max_val, false, -3), find_end_point(f, -0.1, max_val, true, -2), "");
707 plot.add(f, find_end_point(f, 0.1, max_val, false, -4), find_end_point(f, -0.1, max_val, true, -3), "");
708 plot.add(f, find_end_point(f, 0.1, max_val, false, -5), find_end_point(f, -0.1, max_val, true, -4), "");
709 plot.plot("Trigamma", "trigamma.svg", "x", "trigamma(x)");
710
711 f2i = boost::math::polygamma;
712 max_val = -50;
713 plot.clear();
714 plot.add(boost::bind(f2i, 2, _1), find_end_point(boost::bind(f2i, 2, _1), 0.1, max_val, true), 5, "");
715 plot.add(boost::bind(f2i, 2, _1), find_end_point(boost::bind(f2i, 2, _1), 0.1, max_val, true, -1), find_end_point(boost::bind(f2i, 2, _1), -0.1, -max_val, true), "");
716 plot.add(boost::bind(f2i, 2, _1), find_end_point(boost::bind(f2i, 2, _1), 0.1, max_val, true, -2), find_end_point(boost::bind(f2i, 2, _1), -0.1, -max_val, true, -1), "");
717 plot.add(boost::bind(f2i, 2, _1), find_end_point(boost::bind(f2i, 2, _1), 0.1, max_val, true, -3), find_end_point(boost::bind(f2i, 2, _1), -0.1, -max_val, true, -2), "");
718 plot.add(boost::bind(f2i, 2, _1), find_end_point(boost::bind(f2i, 2, _1), 0.1, max_val, true, -4), find_end_point(boost::bind(f2i, 2, _1), -0.1, -max_val, true, -3), "");
719 plot.add(boost::bind(f2i, 2, _1), find_end_point(boost::bind(f2i, 2, _1), 0.1, max_val, true, -5), find_end_point(boost::bind(f2i, 2, _1), -0.1, -max_val, true, -4), "");
720 plot.add(boost::bind(f2i, 2, _1), find_end_point(boost::bind(f2i, 2, _1), 0.1, max_val, true, -6), find_end_point(boost::bind(f2i, 2, _1), -0.1, -max_val, true, -5), "");
721 plot.plot("Polygamma", "polygamma2.svg", "x", "polygamma(2, x)");
722
723 max_val = 800;
724 plot.clear();
725 plot.add(boost::bind(f2i, 3, _1), find_end_point(boost::bind(f2i, 3, _1), 0.1, max_val, false), 5, "");
726 plot.add(boost::bind(f2i, 3, _1), find_end_point(boost::bind(f2i, 3, _1), 0.1, max_val, false, -1), find_end_point(boost::bind(f2i, 3, _1), -0.1, max_val, true), "");
727 plot.add(boost::bind(f2i, 3, _1), find_end_point(boost::bind(f2i, 3, _1), 0.1, max_val, false, -2), find_end_point(boost::bind(f2i, 3, _1), -0.1, max_val, true, -1), "");
728 plot.add(boost::bind(f2i, 3, _1), find_end_point(boost::bind(f2i, 3, _1), 0.1, max_val, false, -3), find_end_point(boost::bind(f2i, 3, _1), -0.1, max_val, true, -2), "");
729 plot.add(boost::bind(f2i, 3, _1), find_end_point(boost::bind(f2i, 3, _1), 0.1, max_val, false, -4), find_end_point(boost::bind(f2i, 3, _1), -0.1, max_val, true, -3), "");
730 plot.add(boost::bind(f2i, 3, _1), find_end_point(boost::bind(f2i, 3, _1), 0.1, max_val, false, -5), find_end_point(boost::bind(f2i, 3, _1), -0.1, max_val, true, -4), "");
731 plot.add(boost::bind(f2i, 3, _1), find_end_point(boost::bind(f2i, 3, _1), 0.1, max_val, false, -6), find_end_point(boost::bind(f2i, 3, _1), -0.1, max_val, true, -5), "");
732 plot.plot("Polygamma", "polygamma3.svg", "x", "polygamma(3, x)");
733 return 0;
734 }
735
736