diff options
| -rw-r--r-- | stokes.hpp | 37 |
1 files changed, 34 insertions, 3 deletions
@@ -5,6 +5,38 @@ #include <iostream> +template <class Real> +Real jacobi(unsigned α, unsigned β, unsigned n, Real x) { + Real f = 0; + + Real y = (x - 1) / 2; + Real z = (x + 1) / 2; + + for (unsigned s = 0; s <= n; s++) { + f += pow(y, n - s) * pow(z, s) / (tgamma(s + 1) * tgamma(n + α - s + 1) * tgamma(β + s + 1) * tgamma(n - s + 1)); + } + + return tgamma(n + α + 1) * tgamma(n + β + 1) * f; +} + +template <class Real> +Real dJacobi(unsigned α, unsigned β, unsigned n, Real x) { + Real f = 0; + + Real y = (x - 1) / 2; + Real z = (x + 1) / 2; + + for (unsigned s = 0; s < n; s++) { + f += pow(y, n - s - 1) * pow(z, s) / (tgamma(s + 1) * tgamma(n + α - s + 1) * tgamma(β + s + 1) * tgamma(n - s - 1 + 1)); + } + + for (unsigned s = 1; s <= n; s++) { + f += pow(y, n - s) * pow(z, s - 1) / (tgamma(s - 1 + 1) * tgamma(n + α - s + 1) * tgamma(β + s + 1) * tgamma(n - s + 1)); + } + + return 0.5 * tgamma(n + α + 1) * tgamma(n + β + 1) * f; +} + template <class Scalar, int ...ps> class Cord { public: @@ -28,12 +60,11 @@ public: } Real gCoeff(unsigned i, Real t) const { - return (1 - t) * t * std::legendre(i, 2 * t - 1); + return (1 - t) * t * jacobi(1, 1, i, 2 * t - 1); } Real dgCoeff(unsigned i, Real t) const { - Real x = 2 * t - 1; - return ((i - 1.0) * x * std::legendre(i, x) - (i + 1) * std::legendre(i + 1, x)) / 2.0; + return (1 - 2 * t) * jacobi(1, 1, i, 2 * t - 1) + 2 * (1 - t) * t * dJacobi(1, 1, i, 2 * t - 1); } Vector<Scalar> f(Real t) const { |