Page MenuHomeHEPForge
Files F7878278

No OneTemporary
Actions

View Options

Index: trunk/tests/test_ClassicalOrthoPolys1D.cc
===================================================================
--- trunk/tests/test_ClassicalOrthoPolys1D.cc (revision 573)
+++ trunk/tests/test_ClassicalOrthoPolys1D.cc (revision 574)
@@ -1,170 +1,171 @@
#include <cmath>
#include "UnitTest++.h"
#include "test_utils.hh"
#include "npstat/nm/ClassicalOrthoPolys1D.hh"
#include "npstat/nm/GaussLegendreQuadrature.hh"
+#include "npstat/nm/FejerQuadrature.hh"
#include "npstat/rng/MersenneTwister.hh"
using namespace npstat;
using namespace std;
namespace {
TEST(LegendreOrthoPoly1D_values)
{
const double eps = 1.0e-15;
LegendreOrthoPoly1D leg;
CHECK(leg.weight(0.123) == 0.5);
CHECK(leg.weight(10.0) == 0.0);
unsigned deg = 0;
CHECK_CLOSE(1.0, leg.poly(deg, 0.3), eps);
deg = 3;
CHECK_CLOSE(-71/343.0L*sqrtl(2*deg+1.0L), leg.poly(deg, 1/7.0L), eps);
deg = 7;
CHECK_CLOSE(-326569/1771561.0L*sqrtl(2*deg+1.0L), leg.poly(deg, 1/11.0L), eps);
GaussLegendreQuadrature glq(10);
for (unsigned i=0; i<5; ++i)
for (unsigned j=0; j<5; ++j)
CHECK_CLOSE((i == j ? 1.0 : 0.0), leg.empiricalKroneckerDelta(glq, i, j), eps);
CHECK_CLOSE(3374/715.0L, leg.integratePoly(3, 4), eps);
unsigned degs[3] = {3, 4, 5};
CHECK_CLOSE(1.051943782704350297, leg.jointIntegral(degs, 3), eps);
}
TEST(JacobiOrthoPoly1D_values)
{
const double eps = 1.0e-15;
JacobiOrthoPoly1D jac(3, 4);
long double x = 1/3.0L;
CHECK_CLOSE(35/32.0L*powl(1-x,3)*powl(1+x,4), jac.weight(x), eps);
GaussLegendreQuadrature glq(10);
for (unsigned i=0; i<5; ++i)
for (unsigned j=0; j<5; ++j)
{
CHECK_CLOSE((i == j ? 1.0 : 0.0), jac.empiricalKroneckerDelta(glq, i, j), eps);
unsigned degs[2];
degs[0] = i;
degs[1] = j;
CHECK_CLOSE((i == j ? 1.0 : 0.0), jac.jointAverage(glq, degs, 2), eps);
}
JacobiOrthoPoly1D jac2(2, 5);
for (unsigned i=0; i<5; ++i)
for (unsigned j=0; j<5; ++j)
{
CHECK_CLOSE((i == j ? 1.0 : 0.0), jac2.empiricalKroneckerDelta(glq, i, j), eps);
unsigned degs[2];
degs[0] = i;
degs[1] = j;
CHECK_CLOSE((i == j ? 1.0 : 0.0), jac2.jointAverage(glq, degs, 2), eps);
}
JacobiOrthoPoly1D jac3(0, 0);
CHECK_CLOSE(0.5, jac3.weight(0.123), eps);
CHECK(jac3.weight(10.0) == 0.0);
unsigned deg = 7;
CHECK_CLOSE(-326569/1771561.0L*sqrtl(2*deg+1.0L), jac3.poly(deg, 1/11.0L), eps);
MersenneTwister rng;
const unsigned npoints = 64;
const unsigned maxdeg = 7;
const unsigned ntries = 5;
std::vector<double> points(npoints);
for (unsigned itry=0; itry<ntries; ++itry)
{
rng.run(&points[0], npoints, npoints);
for (unsigned ipt=0; ipt<npoints; ++ipt)
{
points[ipt] *= 2.0;
points[ipt] -= 1.0;
}
const Matrix<double> avemat(jac.sampleProductAverages(&points[0], npoints, maxdeg));
std::vector<double> averages(maxdeg + 1);
jac.sampleAverages(&points[0], npoints, &averages[0], maxdeg);
CHECK_EQUAL(maxdeg + 1, avemat.nRows());
CHECK_EQUAL(maxdeg + 1, avemat.nColumns());
for (unsigned i=0; i<=maxdeg; ++i)
{
CHECK_CLOSE(averages[i], avemat[i][0], eps);
CHECK_CLOSE(averages[i], avemat[0][i], eps);
for (unsigned j=0; j<=maxdeg; ++j)
{
unsigned degs[2];
degs[0] = i;
degs[1] = j;
const double ave = jac.jointSampleAverage(&points[0], npoints, degs, 2);
CHECK_CLOSE(ave, avemat[i][j], eps);
CHECK_CLOSE(ave, avemat[j][i], eps);
}
}
// Test "twoJointSampleAverages" method
const double degmax = maxdeg + 0.99;
unsigned degs0[3];
for (unsigned i=0; i<3; ++i)
degs0[i] = rng()*degmax;
unsigned degs1[4];
for (unsigned i=0; i<4; ++i)
degs1[i] = rng()*degmax;
const double ave0 = jac.jointSampleAverage(&points[0], npoints, degs0, 3);
const double ave1 = jac.jointSampleAverage(&points[0], npoints, degs1, 4);
const std::pair<double, double>& ave2 = jac.twoJointSampleAverages(
&points[0], npoints, degs0, 3, degs1, 4);
CHECK_CLOSE(ave0, ave2.first, eps);
CHECK_CLOSE(ave1, ave2.second, eps);
}
}
TEST(ChebyshevOrthoPoly2nd_values)
{
const double epsVal = 1.0e-15;
const double epsKron = 1.0e-7;
ChebyshevOrthoPoly2nd cheb;
unsigned deg = 3;
CHECK_CLOSE(-0.54810495626822157434L, cheb.poly(deg, 1/7.0L), epsVal);
deg = 7;
CHECK_CLOSE(-0.66835314371696127673L, cheb.poly(deg, 1/11.0L), epsVal);
GaussLegendreQuadrature glq(1024);
for (unsigned i=0; i<5; ++i)
for (unsigned j=0; j<5; ++j)
CHECK_CLOSE((i == j ? 1.0 : 0.0), cheb.empiricalKroneckerDelta(glq, i, j), epsKron);
}
TEST(ChebyshevOrthoPoly1st_values)
{
const double epsVal = 1.0e-15;
- const double epsKron = 0.005;
+ const double epsKron = 1.0e-3;
const double norm = 1.0/sqrt(2.0);
ChebyshevOrthoPoly1st cheb;
unsigned deg = 0;
CHECK_CLOSE(1.0, cheb.poly(deg, 0.3), epsVal);
deg = 3;
CHECK_CLOSE(-0.41690962099125364431L/norm, cheb.poly(deg, 1/7.0L), epsVal);
deg = 7;
CHECK_CLOSE(-0.59498215518301758629L/norm, cheb.poly(deg, 1/11.0L), epsVal);
- GaussLegendreQuadrature glq(1024);
+ FejerQuadrature fej(1000);
for (unsigned i=0; i<5; ++i)
for (unsigned j=0; j<5; ++j)
- CHECK_CLOSE((i == j ? 1.0 : 0.0), cheb.empiricalKroneckerDelta(glq, i, j), epsKron);
+ CHECK_CLOSE((i == j ? 1.0 : 0.0), cheb.empiricalKroneckerDelta(fej, i, j), epsKron);
}
}

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 5:36 PM (1 d, 12 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805419
Default Alt Text
(6 KB)

Event Timeline