Algorythms/Orthogonal Polynomials/TEST/OrthPolTEST.cpp

#include <gtest/gtest.h>
#include "modules/OrthPol.cpp"
#include <cmath>

#define EPS_CHEB pow(10, -7)
#define EPS_LAG pow(10, -5)

class OrthPolTest : public testing::Test {};

// Test that OrthPol returns correct values if x=0
// TEST_F(OrthPolTest, ChebPolZeroCase) {
//    double zero = 0.0;
//    for (int i=0; i<=50; ++i) 
//    {

//       std::vector<double> res = OrthPol(1, i, 0);
//       if (i % 2 == 1)
//       {
//         EXPECT_DOUBLE_EQ(res[0], zero) << "Chebyshev polinom of the first kind " << i << "th-order are not equal 0";
//         EXPECT_DOUBLE_EQ(res[1], pow(-1.0, (i-1)/2)*i) << "Derivative of Chebyshev polinom of the first kind " << i << "th-order are not equal "<< i << "*(-1)^" << (i-1)/2;
//       }
//       else
//       {
//         EXPECT_DOUBLE_EQ(res[0], pow(-1.0, i/2)) << "Chebyshev polinom of the first kind " << i << "th-order are not equal (-1)^" << i/2;
//         EXPECT_DOUBLE_EQ(res[1], zero) << "Derivative of Chebyshev polinom of the first kind " << i << "th-order are not equal 0";
//       }
//    }
//    for (int i=0; i<=50; ++i) 
//    {
//       std::vector<double> res = OrthPol(2, i, 0);
//       if (i % 2 == 1)
//       {
//         EXPECT_DOUBLE_EQ(res[0], zero) << "Chebyshev polinom of the second kind " << i << "th-order are not equal 0";
//         EXPECT_DOUBLE_EQ(res[1], pow(-1.0, (i-1)/2)*(i+1)) << "Derivative of Chebyshev polinom of the second kind " << i << "th-order are not equal "<< (i+1) << "*(-1)^" << (i-1)/2;
//       }
//       else
//       {
//         EXPECT_DOUBLE_EQ(res[0], pow(-1.0, i/2)) << "Chebyshev polinom of the second kind " << i << "th-order are not equal (-1)^" << i/2;
//         EXPECT_DOUBLE_EQ(res[1], zero) << "Derivative of Chebyshev polinom of the second kind " << i << "th-order are not equal 0";
//       }
//    }
// }

// Test that OrthPol returns correct values if x=+-1
// TEST_F(OrthPolTest, ChebPolOneCase) {
//    for (int i=0; i<=50; ++i) 
//    {
//     std::vector<double> res = OrthPol(1, i, 1);
//     EXPECT_DOUBLE_EQ(res[0], 1) << "Chebyshev polinom of the first kind " << i << "th-order are not equal " << 1;
//     EXPECT_DOUBLE_EQ(res[1], pow(i, 2)) << "Derivative of Chebyshev polinom of the first kind " << i << "th-order are not equal "<< pow(i, 2);

//     res = OrthPol(2, i, 1);
//     EXPECT_DOUBLE_EQ(res[0], i+1) << "Chebyshev polinom of the second kind " << i << "th-order are not equal " << i+1;
//     EXPECT_DOUBLE_EQ(res[1], i*(i+1)*(i+2)/3) << "Derivative of Chebyshev polinom of the first kind " << i << "th-order are not equal "<< i*(i+1)*(i+2)/3;

//     res = OrthPol(1, i, -1);
//     EXPECT_DOUBLE_EQ(res[0], pow(-1, i)) << "Chebyshev polinom of the first kind " << i << "th-order are not equal " << pow(-1, i);
//     EXPECT_DOUBLE_EQ(res[1], pow(-1, i-1) * pow(i, 2)) << "Derivative of Chebyshev polinom of the first kind " << i << "th-order are not equal "<< pow(-1, i-1) * pow(i, 2);

//     res = OrthPol(2, i, -1);
//     EXPECT_DOUBLE_EQ(res[0], pow(-1, i)*(i+1)) << "Chebyshev polinom of the first kind " << i << "th-order are not equal " << pow(-1, i)*(i+1);
//     EXPECT_DOUBLE_EQ(res[1], pow(-1, i+1) * (i*(i+1)*(i+2)/3)) << "Derivative of Chebyshev polinom of the first kind " << i << "th-order are not equal "<< pow(-1, i+1) * (i*(i+1)*(i+2)/3);
//    }
// }

double ExplicitExpressionCheb(const int& type, const int& n, const double& x)
{
    if (n == 0)
    {
        return 1;
    }
    if (type == 1)
    {
        if (n == 1)
        {
            return x;
        }
        int N = n / 2;
        double d = static_cast<double>(n) / 2;
        double c = pow(-1, 0) * std::tgamma(n) / (std::tgamma(1) * std::tgamma(n + 1));
        double g = pow(2*x, n);
        double res = 0;
        for (int m=0; m<=N; ++m)
        {
            c = pow(-1, m) * std::tgamma(n - m) / (std::tgamma(m + 1) * std::tgamma(n - 2*m + 1));
            g = pow(2*x, n - 2*m);
            res += c * g;
        }
    return d * res;
    }
    else if (type == 2)
    {
        if (n == 1)
        {
            return 2*x;
        }
        int N = n / 2;
        double d = 1;
        double c = 1;
        double g = pow(2*x, n);
        double res = 0;
        for (int m=0; m<=N; ++m)
        {
            c = pow(-1, m) * std::tgamma(n - m + 1) / (std::tgamma(m + 1) * std::tgamma(n - 2*m + 1));
            g = pow(2*x, n - 2*m);
            res += c * g;
        }
    return d * res;
    }
    else
    {
        printf("The order must be 1 or 2!");
        std::exit(9);
    }
}

double ExplicitExpressionLaguerre(const int& n, const double& x)
{
    if (n == 0)
    {
        return 1;
    }
    else if (n == 1)
    {
        return 1 - x;
    }
    int N = n;
    double d = 1;
    double c = pow(-1, 0) * std::tgamma(n) / (std::tgamma(1) * std::tgamma(n + 1));
    double g = pow(2*x, n);
    double res = 0;
    for (int m=0; m<=N; ++m)
    {
        c = pow(-1, m) * std::tgamma(n + 1) / (pow(std::tgamma(m + 1), 2) * std::tgamma(n - m + 1));
        g = pow(x, m);
        res += c * g;
    }
    return d * res;
}


TEST_F(OrthPolTest, ExplicitExpressionCase){
    double x = -1; // the left boundary
    double length = 2; //the length of the segmment
    double step = 0.2;
    std::vector<double> res_rec;
    double res_expl;
    for (int i=0; i <= static_cast<int>(length / step); ++i)
    {
        for (int j=0; j <= 20; ++j)
        {
            res_rec = OrthPol(1, j, x);
            res_expl = ExplicitExpressionCheb(1, j, x);
            // std::cout << j << " " << x << " " << res_rec[0] << "\n";
            EXPECT_NEAR(res_rec[0], res_expl, EPS_CHEB);

            res_rec = OrthPol(2, j, x);
            res_expl = ExplicitExpressionCheb(2, j, x);
            // std::cout << j << " " << x << " " << res_rec[0] << "\n";
            EXPECT_NEAR(res_rec[0], res_expl, EPS_CHEB);

        }
        x += step;
    }

    x = -10;
    length = 20;
    step = 0.5;
    for (int i=0; i <= static_cast<int>(length / step); ++i)
    {
        for (int j=0; j <= 20; ++j)
        {
            res_rec = OrthPol(3, j, x);
            res_expl = ExplicitExpressionLaguerre(j, x);
            // std::cout << j << " " << x << " " << res_rec[0] << "\n";
            EXPECT_NEAR(res_rec[0], res_expl, EPS_LAG);
        }
        x += step;
    }


}



// int main(int argc, char **argv) {
//   testing::InitGoogleTest(&argc, argv);
//   return RUN_ALL_TESTS();
// }