00001 #include "SundanceGauss1D.hpp"
00002 #ifdef _MSC_VER
00003 # include "winmath.h"
00004 #endif
00005
00006 using namespace Sundance;
00007 using namespace Sundance;
00008 using namespace Sundance;
00009 using namespace Sundance;
00010 using namespace Sundance;
00011 using namespace Teuchos;
00012
00013 Gauss1D::Gauss1D(int n)
00014 : nodes_(n), weights_(n)
00015 {
00016 computeWeights(n, -1.0, 1.0);
00017 }
00018
00019 Gauss1D::Gauss1D(int n, double a, double b)
00020 : nodes_(n), weights_(n)
00021 {
00022 computeWeights(n, a, b);
00023 }
00024
00025
00026
00027 void Gauss1D::computeWeights(int n, double a, double b)
00028 {
00029 int m = (n+1)/2;
00030
00031 double xMid = (b+a)/2.0;
00032 double halfWidth = (b-a)/2.0;
00033
00034 for (int i=0; i<m; i++)
00035 {
00036
00037 double z = cos(M_PI*(i+0.75)/(n+0.5));
00038 double dP;
00039 double zOld;
00040 double tol = 1.0e-14;
00041
00042 do
00043 {
00044 double p1 = 1.0;
00045 double p2 = 0.0;
00046 for (int j=1; j<=n; j++)
00047 {
00048 double p3 = p2;
00049 p2 = p1;
00050 p1 = ((2.0*j-1.0)*z*p2 - (j-1.0)*p3)/j;
00051 }
00052 dP = n*(z*p1-p2)/(z*z-1.0);
00053 zOld =z;
00054 z = zOld - p1/dP;
00055 }
00056 while ( fabs(z-zOld) > tol );
00057 nodes_[i] = xMid - halfWidth*z;
00058 nodes_[n-i-1] = xMid + halfWidth*z;
00059 weights_[i] = 2.0*halfWidth/((1.0-z*z)*dP*dP);
00060 weights_[n-i-1] = weights_[i];
00061 }
00062 }
00063
00064
00065 bool Gauss1D::unitTest()
00066 {
00067 std::cerr << "------------------ Gauss1D unit test ----------------------"
00068 << std::endl;
00069
00070 Gauss1D q(20, 0.0, M_PI);
00071
00072 double sum = 0.0;
00073 for (int i=0; i<q.nPoints(); i++)
00074 {
00075 sum += q.weights()[i]*sin(q.nodes()[i]);
00076 }
00077 std::cerr << "integral of sin(x) over [0, pi] = " << sum << std::endl;
00078 double sumErr = fabs(sum - 2.0);
00079 bool sumPass = sumErr < 1.0e-10;
00080 std::cerr << "error = " << sumErr << std::endl;
00081 if (sumPass) std::cerr << "Gauss1D sine test PASSED" << std::endl;
00082 else std::cerr << "Gauss1D sine test FAILED" << std::endl;
00083 std::cerr << "------------------ End Gauss1D unit test ----------------------"
00084 << std::endl;
00085 return sumPass;
00086 }
00087
00088
