00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031 #include "PlayaExceptions.hpp"
00032 #include "SundanceGaussLobatto1D.hpp"
00033 #ifdef _MSC_VER
00034 # include "winmath.h"
00035 #endif
00036
00037 using namespace Sundance;
00038 using namespace Teuchos;
00039
00040 GaussLobatto1D::GaussLobatto1D(int n) :
00041 nodes_(n), weights_(n)
00042 {
00043 computeWeights(n, -1.0, 1.0);
00044 }
00045
00046 GaussLobatto1D::GaussLobatto1D(int n, double a, double b) :
00047 nodes_(n), weights_(n)
00048 {
00049 computeWeights(n, a, b);
00050 }
00051
00052 void GaussLobatto1D::computeWeights(int n, double a, double b)
00053 {
00054
00055 TEUCHOS_TEST_FOR_EXCEPTION(n < 2, std::runtime_error, "number of points=" << n
00056 << " must be at least 2 for Gauss-Lobatto-Legendre quadrature!");
00057
00058 int m = (n + 1) / 2;
00059
00060 double xMid = (b + a) / 2.0;
00061 double halfWidth = (b - a) / 2.0;
00062
00063 for (int i = 0; i < m; i++)
00064 {
00065
00066 double z = cos(M_PI * i / (n - 1));
00067 double p1;
00068 double zOld;
00069 double tol = 1.0e-14;
00070
00071 do
00072 {
00073 p1 = 1.0;
00074 double p2 = 0.0;
00075 for (int j = 1; j < n; j++)
00076 {
00077 double p3 = p2;
00078 p2 = p1;
00079 p1 = ((2.0 * j - 1.0) * z * p2 - (j - 1.0) * p3) / j;
00080 }
00081
00082 zOld = z;
00083
00084
00085
00086
00087
00088 z = zOld - (z * p1 - p2) / (n * p1);
00089 } while (fabs(z - zOld) > tol);
00090
00091 if (i == 0)
00092 {
00093 nodes_[0] = a;
00094 nodes_[n - 1] = b;
00095 }
00096 else
00097 {
00098 nodes_[i] = xMid - halfWidth * z;
00099 nodes_[n - i - 1] = xMid + halfWidth * z;
00100 }
00101 weights_[i] = 2.0 * halfWidth / ((n - 1) * n * p1 * p1);
00102 weights_[n - i - 1] = weights_[i];
00103 }
00104 }
00105
00106 bool GaussLobatto1D::unitTest()
00107 {
00108 std::cerr
00109 << "------------------ GaussLobatto1D unit test ----------------------"
00110 << std::endl;
00111
00112 GaussLobatto1D q(20, 0.0, M_PI);
00113
00114 double sum = 0.0;
00115 for (int i = 0; i < q.nPoints(); i++)
00116 {
00117 sum += q.weights()[i] * sin(q.nodes()[i]);
00118 }
00119 std::cerr << "integral of sin(x) over [0, pi] = " << sum << std::endl;
00120 double sumErr = fabs(sum - 2.0);
00121 bool sumPass = sumErr < 1.0e-10;
00122 std::cerr << "error = " << sumErr << std::endl;
00123 if (sumPass)
00124 std::cerr << "GaussLobatto1D sine test PASSED" << std::endl;
00125 else
00126 std::cerr << "GaussLobatto1D sine test FAILED" << std::endl;
00127 std::cerr
00128 << "------------------ End GaussLobatto1D unit test ----------------------"
00129 << std::endl;
00130 return sumPass;
00131 }
00132