Skip to content

Commit

Permalink
added gauss_lobattor tests
Browse files Browse the repository at this point in the history
  • Loading branch information
@BryanFlynt
BryanFlynt committed Sep 4, 2021
1 parent 2605554 commit 4313517
Show file tree
Hide file tree
Showing 8 changed files with 268 additions and 114 deletions.
4 changes: 2 additions & 2 deletions include/polycalc/quadrature/gauss_jacobi.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,8 +9,8 @@
#include <cassert>
#include <vector>

#include "parameters.hpp"
#include "polynomial/jacobi.hpp"
#include "polycalc/parameters.hpp"
#include "polycalc/polynomial/jacobi.hpp"

namespace polycalc {
namespace quadrature {
Expand Down
153 changes: 110 additions & 43 deletions include/polycalc/quadrature/gauss_lobatto.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,8 +9,8 @@
#include <cassert>
#include <vector>

#include "parameters.hpp"
#include "polynomial/jacobi.hpp"
#include "polycalc/parameters.hpp"
#include "polycalc/polynomial/jacobi.hpp"

namespace polycalc {
namespace quadrature {
Expand Down Expand Up @@ -53,21 +53,56 @@ template <typename T, typename P>
std::vector<typename GaussLobatto<T, P>::value_type> GaussLobatto<T, P>::zeros(const unsigned n) const {
assert(n > 0);

// Good Decimal Calculator found at following site
// https://keisan.casio.com/exec/system/1280801905

// Zeros to return
std::vector<value_type> x(n);

if (1 == n) {
x[0] = 0.0;
} else if (2 == n) {
x[0] = -1.0;
x[1] = +1.0;
} else {
polynomial jac(alpha_ + 1, beta_ + 1);
auto zeros = jac.zeros(n - 2);

x[0] = -1.0;
std::copy(zeros.begin(), zeros.end(), x.begin() + 1);
x[n - 1] = +1.0;
switch (n) {
case 1:
x[0] = 0.0;
break;
case 2:
x[0] = -1.0;
x[1] = +1.0;
break;
case 3:
x[0] = -1.0;
x[1] = 0.0;
x[2] = +1.0;
break;
case 4:
x[0] = -1.0;
x[1] = -std::sqrt(5.0L) / 5.0L;
x[2] = +std::sqrt(5.0L) / 5.0L;
x[3] = +1.0;
break;
case 5:
x[0] = -1.0;
x[1] = -std::sqrt(21.0L) / 7.0L;
x[2] = 0.0;
x[3] = +std::sqrt(21.0L) / 7.0L;
x[4] = +1.0;
break;
case 6:
x[0] = -1.0;
x[1] = -std::sqrt((7.0L+2.0L*std::sqrt(7.0L))/21.0L);
x[2] = -std::sqrt((7.0L-2.0L*std::sqrt(7.0L))/21.0L);
x[3] = +std::sqrt((7.0L-2.0L*std::sqrt(7.0L))/21.0L);
x[4] = +std::sqrt((7.0L+2.0L*std::sqrt(7.0L))/21.0L);
x[5] = +1.0;
break;
default:
polynomial jac(alpha_ + 1, beta_ + 1);
auto zeros = jac.zeros(n - 2);

x.front() = -1.0;
std::copy(zeros.begin(), zeros.end(), x.begin() + 1);
x.back() = +1.0;
if (!(n % 2 == 0)) {
x[std::ptrdiff_t(n / 2)] = 0; // Correct 10E-16 error at zero
}
}
return x;
}
Expand All @@ -76,38 +111,70 @@ template <typename T, typename P>
std::vector<typename GaussLobatto<T, P>::value_type> GaussLobatto<T, P>::weights(const unsigned n) const {
assert(n > 0);

// Good Decimal Calculator found at following site
// https://keisan.casio.com/exec/system/1280801905

// Weights to return
std::vector<value_type> w(n);

if (1 == n) {
w[0] = 2.0;
} else if (2 == n) {
w[0] = 1.0;
w[1] = 1.0;
} else {
// Get location of zeros
auto z = this->zeros(n);

// Evaluate Jacobi n-1 polynomial at each zero
polynomial jac(alpha_, beta_);
std::vector<value_type> w(n);
for (size_type i = 0; i < n; ++i) {
w[i] = jac.eval(n - 1, z[i]);
}

const value_type one = 1;
const value_type two = 2;
const value_type apb = alpha_ + beta_;

value_type fac;
fac = std::pow(two, apb + one) * std::tgamma(alpha_ + n) * std::tgamma(beta_ + n);
fac /= (n - 1) * std::tgamma(n) * std::tgamma(alpha_ + beta_ + n + one);

for (size_type i = 0; i < n; ++i) {
w[i] = fac / (w[i] * w[i]);
}
w[0] *= (beta_ + one);
w[n - 1] *= (alpha_ + one);
switch (n) {
case 1:
w[0] = +2.0;
break;
case 2:
w[0] = +1.0;
w[1] = +1.0;
break;
case 3:
w[0] = 1.0L / 3.0L;
w[1] = 4.0L / 3.0L;
w[2] = 1.0L / 3.0L;
break;
case 4:
w[0] = 1.0L / 6.0L;
w[1] = 5.0L / 6.0L;
w[2] = 5.0L / 6.0L;
w[3] = 1.0L / 6.0L;
break;
case 5:
w[0] = 1.0L / 10.0L;
w[1] = 49.0L / 90.0L;
w[2] = 32.0L / 45.0L;
w[3] = 49.0L / 90.0L;
w[4] = 1.0L / 10.0L;
break;
case 6:
w[0] = 1.0L / 15.0L;
w[1] = (14.0L-std::sqrt(7.0L))/30.0L;
w[2] = (14.0L+std::sqrt(7.0L))/30.0L;
w[3] = (14.0L+std::sqrt(7.0L))/30.0L;
w[4] = (14.0L-std::sqrt(7.0L))/30.0L;
w[5] = 1.0L / 15.0L;
break;
default:

// Get location of zeros
auto z = this->zeros(n);

// Evaluate Jacobi n-1 polynomial at each zero
polynomial jac(alpha_, beta_);
for (size_type i = 0; i < n; ++i) {
w[i] = jac.eval(n - 1, z[i]);
}

const value_type one = 1;
const value_type two = 2;
const value_type apb = alpha_ + beta_;

value_type fac;
fac = std::pow(two, apb + one) * std::tgamma(alpha_ + n) * std::tgamma(beta_ + n);
fac /= (n - 1) * std::tgamma(n) * std::tgamma(alpha_ + beta_ + n + one);

for (size_type i = 0; i < n; ++i) {
w[i] = fac / (w[i] * w[i]);
}
w[0] *= (beta_ + one);
w[n - 1] *= (alpha_ + one);
}
return w;
}
Expand Down
4 changes: 2 additions & 2 deletions include/polycalc/quadrature/gauss_radau.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,8 +9,8 @@
#include <cassert>
#include <vector>

#include "parameters.hpp"
#include "polynomial/jacobi.hpp"
#include "polycalc/parameters.hpp"
#include "polycalc/polynomial/jacobi.hpp"

namespace polycalc {
namespace quadrature {
Expand Down
3 changes: 2 additions & 1 deletion test/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,4 +44,5 @@ endfunction()


add_catch_test(jacobi)
add_catch_test(lagrange)
add_catch_test(lagrange)
add_catch_test(gauss_lobatto)
108 changes: 108 additions & 0 deletions test/gauss_lobatto.cpp
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,108 @@


#include <cmath>
#include <iomanip>
#include <iostream>

#include "polycalc/quadrature/gauss_lobatto.hpp"


#include "catch.hpp"
#include "utilities.hpp"


TEST_CASE("Jacobi Polynomial", "[default]") {
using namespace ::polycalc::quadrature;
using value_type = double;

constexpr value_type zero(0);
constexpr value_type one(1);
constexpr value_type two(2);
constexpr value_type three(3);
constexpr value_type four(4);
constexpr value_type five(5);
constexpr value_type six(6);
constexpr value_type seven(7);
constexpr value_type eight(8);
constexpr value_type nine(9);
constexpr value_type ten(10);

SECTION("Legendre Expansion N = 3") {
GaussLobatto<value_type> gll(0,0);

const int N = 3;
auto z = gll.zeros(N);
auto w = gll.weights(N);

REQUIRE(almost_equal(z[0], -one));
REQUIRE(almost_equal(z[1], zero));
REQUIRE(almost_equal(z[2], +one));

REQUIRE(almost_equal(w[0], one/three));
REQUIRE(almost_equal(w[1], four/three));
REQUIRE(almost_equal(w[2], one/three));
}

SECTION("Legendre Expansion N = 4") {
GaussLobatto<value_type> gll(0,0);

const int N = 4;
auto z = gll.zeros(N);
auto w = gll.weights(N);

REQUIRE(almost_equal(z[0], -one));
REQUIRE(almost_equal(z[1], -std::sqrt(five)/five));
REQUIRE(almost_equal(z[2], +std::sqrt(five)/five));
REQUIRE(almost_equal(z[3], +one));

REQUIRE(almost_equal(w[0], one/six));
REQUIRE(almost_equal(w[1], five/six));
REQUIRE(almost_equal(w[2], five/six));
REQUIRE(almost_equal(w[3], one/six));
}

SECTION("Legendre Expansion N = 5") {
GaussLobatto<value_type> gll(0,0);

const int N = 5;
auto z = gll.zeros(N);
auto w = gll.weights(N);

REQUIRE(almost_equal(z[0], -one));
REQUIRE(almost_equal(z[1], -std::sqrt(21.0)/seven));
REQUIRE(almost_equal(z[2], zero));
REQUIRE(almost_equal(z[3], +std::sqrt(21.0)/seven));
REQUIRE(almost_equal(z[4], +one));

REQUIRE(almost_equal(w[0], one/ten));
REQUIRE(almost_equal(w[1], 49./90.));
REQUIRE(almost_equal(w[2], 32./45.));
REQUIRE(almost_equal(w[3], 49./90.));
REQUIRE(almost_equal(w[4], one/ten));
}

SECTION("Legendre Expansion N = 7") {
GaussLobatto<value_type> gll(0,0);

const int N = 7;
auto z = gll.zeros(N);
auto w = gll.weights(N);

REQUIRE(almost_equal(z[0], -1.0));
REQUIRE(almost_equal(z[1], -0.830223896278566929872));
REQUIRE(almost_equal(z[2], -0.4688487934707142138038));
REQUIRE(almost_equal(z[3], 0.0));
REQUIRE(almost_equal(z[4], +0.4688487934707142138038));
REQUIRE(almost_equal(z[5], +0.830223896278566929872));
REQUIRE(almost_equal(z[6], +1.0));

REQUIRE(almost_equal(w[0], 0.04761904761904761904762));
REQUIRE(almost_equal(w[1], 0.2768260473615657)); // Not exact
REQUIRE(almost_equal(w[2], 0.4317453812098626234179));
REQUIRE(almost_equal(w[3], 0.487619047619047619048));
REQUIRE(almost_equal(w[4], 0.4317453812098626234179));
REQUIRE(almost_equal(w[5], 0.2768260473615657)); // Not exact
REQUIRE(almost_equal(w[6], 0.04761904761904761904762));
}

}
34 changes: 1 addition & 33 deletions test/jacobi.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,40 +7,8 @@
#include <iostream>

#include "catch.hpp"
#include "utilities.hpp"

template <class T>
typename std::enable_if<!std::numeric_limits<T>::is_integer, bool>::type almost_equal(T x, T y) {
constexpr T eps = std::numeric_limits<T>::epsilon();
constexpr T zero = 0;
const T abs_diff = std::abs(x - y);

// std::cout << std::setprecision(16);
// std::cout << "x = " << x << std::endl;
// std::cout << "y = " << y << std::endl;
// std::cout << "std::abs(x - y) = " << abs_diff << std::endl;
// std::cout << "eps = " << eps << std::endl;
// std::cout << "eps * abs(x) = " << eps * std::abs(x) << std::endl;
// std::cout << "eps * abs(y) = " << eps * std::abs(y) << std::endl;

if ((x == zero) || (y == zero)) {
if (abs_diff <= eps) {
return true;
}
} else {
if ((abs_diff <= (eps * std::abs(x))) || (abs_diff <= (eps * std::abs(y)))) {
return true;
}
}
return false;
}

// template <typename T>
// void display_vector(const std::vector<T>& vec) {
// std::cout << std::setprecision(16);
// for (auto i = 0; i < vec.size(); ++i) {
// std::cout << vec[i] << std::endl;
// }
// }

TEST_CASE("Jacobi Polynomial", "[default]") {
using namespace ::polycalc::polynomial;
Expand Down
Loading

0 comments on commit 4313517

Please sign in to comment.