76 Array< RCP<const OneDOrthogPolyBasis<int,double> > > qBases(1);
77 qBases[0] = rcp(
new LegendreBasis<int,double>(quadOrder));
79 RCP<const CompletePolynomialBasis<int,double> > qBasis =
80 rcp(
new CompletePolynomialBasis<int,double>(qBases));
82 quad_ = rcp(
new TensorProductQuadrature<int, double>(qBasis, quadOrder));
97 JacobiBasis<int, double> basis(nMax, alpha, beta,
false);
99 Array<Array<double> > qp =
quad_->getQuadPoints();
100 Array<double> qw =
quad_->getQuadWeights();
103 for (
int n=0; n<=nMax; n++)
105 double nFact = tgamma(n+1.0);
106 cout <<
"n=" << n << endl;
107 for (
double x=-1.0; x<=1.0; x+=0.25)
109 cout << setw(20) << x << setw(20) << basis.evaluate(x, n)
112 for (
int m=0; m<=nMax; m++)
115 for (
int q=0; q<qw.size(); q++)
118 double w = qw[q] * pow(1-x,alpha)*pow(1+x,beta);
119 double Pn = basis.evaluate(x, n);
120 double Pm = basis.evaluate(x, m);
125 exact = pow(2.0, alpha+beta+1.0)/(2.0*n+alpha+beta+1.0)
126 * tgamma(n+alpha+1.0)*tgamma(n+beta+1.0)
127 /tgamma(n+alpha+beta+1.0)/nFact;
128 double err =
fabs(exact - sum);
129 cout << setw(4) << n << setw(4) << m
130 << setw(20) << exact << setw(20) << sum << setw(20) << err
138 cout <<
"***** FAIL ******" << endl;