Integration with quadrature and multiprecision boost libraries in C++

Question

after searching i had found amazing code for integration by

quadrature boost library.

rather than

log(x)/(1+x)

want to integrate

(poly[0]+poly[1]*x+poly[2]*x^2+...+poly[n]*x^n)*log(x)/(1+x). But, i do not

know how to insert the vector

poly

to

struct f

or even how to call these operators from main function. The code :

#include<iostream>
#include<boost/math/constnats/constants.hpp>
#include<boost/multiprecision/cpp_dec_float.hpp>
#include <boost/numeric/quadrature/adaptive.hpp>
#include <boost/numeric/quadrature/kronrodgauss.hpp>
#include <boost/numeric/quadrature/epsilon.hpp>
using namespace std;
using boost::multiprecision::cpp_dec_float_50;
namespace quadrature=boost::numeric::quadrature;
struct f
{
 double operator()(double x) const {
 return (log(x)/(1+x); }
};
int main()
{ 
vector<cpp_dec_float_50> poly(0);
cpp_dec_float_50 p = 0;
for (int i=0;i<=n;i++)
{
   p=polynomial(i,n);
   poly.push_back(p); 
}

double answer,error_estimate;
quadrature::adaptive().relative_accuracy(1e-5).absolute_accuracy(1e-7)
(f(),0.,1.,answer,error_estimate);
cout<<"ans"<<answer<<endl;
return 0;
}
cpp_dec_float_50 polynomial(int k ,int n)
{
.
.
.

}

Also, when changing the double operator, to cpp_dec_float_50 operator in

struct f

many problems arise. and the later type is necessary in my project. Any one can fix that ?

EDIT

i tried this, but i do sth wrong

#include<iostream>
#include <boost/numeric/quadrature/adaptive.hpp>
#include <boost/numeric/quadrature/kronrodgauss.hpp>
#include <boost/numeric/quadrature/epsilon.hpp>
#include<boost/math/constants/constants.hpp>
#include<boost/multiprecision/cpp_dec_float.hpp>
using namespace std;
using boost::multiprecision::cpp_dec_float_50;
namespace quadrature=boost::numeric::quadrature;
double polynomial(int k ,int n);
struct f
{ const cpp_dec_float_50 s=0;
 vector<cpp_dec_float_50> poly;
 cpp_dec_float_50 sum()const{
 for(int i=0;i<=poly.size();i++)
  s+=poly[i];
 return s

 }
 double operator()(double x) const {
 return
 s*log(x)/(1+x); }
 };
  int main()
  {
int n=2;
 f fun;
 cpp_dec_float_50 p = 0;
 for (int i=0;i<=n;i++)
 {
 p=polynomial(i,n);
 fun.poly.push_back(p);
 }

 double answer,error_estimate;
 quadrature::adaptive().relative_accuracy(1e-5).absolute_accuracy(1e-7)
 (fun,0.,1.,answer,error_estimate);
 cout<<"ans"<<answer<<endl;
 return 0;
 }
 double polynomial(int k ,int n)
 {
  return k;

  }

Edit when using Patstew suggestion Two errors occur

mmm firstly, how to put the above expression in the struct in order to integrate it; — Sarah, May 18 '16 at 10:19

patstew · Answer 1 · 2016-05-20T15:40:39.000

0

Try something along the lines of:

struct f
{
 vector<cpp_dec_float_50> poly;
 double operator()(double x) const {
 return (poly[0]+poly[1]*x+poly[2]*x^2+...+poly[n]*x^n)*log(x)/(1+x); }
};
int main()
{ 
f fun;
cpp_dec_float_50 p = 0;
for (int i=0;i<=n;i++)
{
   p=polynomial(i,n);
   fun.poly.push_back(p); 
}

double answer,error_estimate;
quadrature::adaptive().relative_accuracy(1e-5).absolute_accuracy(1e-7)
(fun,0.,1.,answer,error_estimate);
cout<<"ans"<<answer<<endl;
return 0;
}

EDIT: RE you own answer, you never call sum (and s is const so you couldn't change it if you did) so s is always 0 and you will always get 0 as your answer. Also you are iterating all the way up to poly.size() in sum(), but poly[poly.size()-1] is the last element. I think you really want your sum function to calculate a polynomial? Try this:

#include<iostream>
#include <boost/numeric/quadrature/adaptive.hpp>
#include <boost/numeric/quadrature/kronrodgauss.hpp>
#include <boost/numeric/quadrature/epsilon.hpp>
#include<boost/math/constants/constants.hpp>
#include<boost/multiprecision/cpp_dec_float.hpp>
using namespace std;
using boost::multiprecision::cpp_dec_float_50;
namespace quadrature=boost::numeric::quadrature;
double polynomial(int k ,int n);
struct f
{
 vector<double> poly;
 double polysum(double x) {
 double s = poly[0];
 double p = 1;
 for(int i=1;i<poly.size();i++) {
  p = p*x;
  s+= p*poly[i];
 }
 return s

 }
 double operator()(double x) {
 return polysum(x)*log(x)/(1+x); }
 };
  int main()
  {
int n=2;
 f fun;
 double p = 0;
 for (int i=0;i<=n;i++)
 {
 p=polynomial(i,n);
 fun.poly.push_back(p);
 }

 double answer,error_estimate;
 quadrature::adaptive().relative_accuracy(1e-5).absolute_accuracy(1e-7)
 (fun,0.,1.,answer,error_estimate);
 cout<<"ans"<<answer<<endl;
 return 0;
 }
 double polynomial(int k ,int n)
 {
  return k;

  }

edited May 20 '16 at 15:40

answered May 19 '16 at 21:05

patstew

1,806
17
21

it seems god, see plz my trial – Sarah May 19 '16 at 21:31
This is exactly what i want. But, there still two errors. one of them related to use double operator, were the polysum is cpp_dec_float_50 – Sarah May 20 '16 at 13:44
It look like your quadrature function isn't expecting a const function, so remove the const from operator(). Also, do you really need multiprecision numbers? double is usually sufficient. – patstew May 20 '16 at 15:43
it really works! just add const like this: double polysum(double x)const – Sarah May 20 '16 at 17:01
i need to represent entries with no. of digits exceeds 200 digit, does double support this no. of digits ? i mean to convert cpp_dec_float_50 to double. not o use double from the scratch – Sarah May 20 '16 at 17:04
What on earth are you doing that needs that level of precision? You can count the number of protons it takes stacked end to end to stretch across the universe with 45 digits... Anyway, a double is accurate to ~16 decimal digits. You need to be very careful about using high precision numbers, if you mix in a lower precision number like a double (e.g. some of your original functions returned double) you can ruin the accuracy of your answer. – patstew May 20 '16 at 17:13
these are the coefficients of chebyshev polynomial which are in terms of factorial function which increases in large way as the degree of polynomial increase. No solution to this problem ? – Sarah May 20 '16 at 17:24
Doubles are fine for storing large numbers - it can represent 10^200. It's just that only the first ~16 figures are stored. So problems mainly arise when you add small numbers and big numbers. e.g. (10^200 + 1) - 10^200 will give you 0 rather than 1. But 10^200 * 10^-200 should give you 1 as expected, multiplication is fine. Normally you can either decide that the small number is irrelevant and ignore it, or rejig your equations so that the big numbers cancel out first. – patstew May 20 '16 at 17:48
No, this will not be fine to my issue.Thank u very much Patstew i voted up for u but they will not counter it. For his problem i will ask again about this – Sarah May 20 '16 at 17:59

Integration with quadrature and multiprecision boost libraries in C++

EDIT

1 Answers1