Add tests for Chebyshev and Laguerre polynomials

Orthogonal Polynomials/CMakeLists.txt

@@ -0,0 +1,17 @@
+cmake_minimum_required(VERSION 3.12)
+
+# set( CMAKE_CXX_COMPILER "C:/msys64/mingw64/bin/g++.exe" )
+# set( CMAKE_C_COMPILER "C:/msys64/mingw64/bin/gcc.exe" )
+
+project(
+    OrthPol
+    LANGUAGES CXX
+)
+
+set(CMAKE_CXX_STANDARD 20)
+
+add_subdirectory(modules)
+add_subdirectory(bin)
+
+enable_testing()
+add_subdirectory(TEST)

Orthogonal Polynomials/TEST/CMakeLists.txt

@@ -0,0 +1,54 @@
+include(FetchContent)
+
+FetchContent_Declare(
+  googletest
+  GIT_REPOSITORY https://github.com/google/googletest.git
+  GIT_TAG release-1.12.1
+)
+
+# For Windows: Prevent overriding the parent project's compiler/linker settings
+set(gtest_force_shared_crt ON CACHE BOOL "" FORCE)
+FetchContent_MakeAvailable(googletest)
+
+enable_testing()
+
+add_executable(
+  OrthPolTEST
+  OrthPolTEST.cpp
+)
+
+target_link_libraries(
+  OrthPolTEST
+  GTest::gtest_main
+)
+
+target_include_directories(OrthPolTEST PUBLIC ${PROJECT_SOURCE_DIR})
+
+include(GoogleTest)
+
+gtest_discover_tests(OrthPolTEST)
+
+
+# cmake_minimum_required(VERSION 3.0)
+
+# project(OrthPolTEST)
+
+# find_package(GTest REQUIRED)
+# find_package(Threads REQUIRED)
+
+# set(CMAKE_CXX_STANDARD 11)
+# set(CMAKE_CXX_STANDARD_REQUIRED on)
+
+# include_directories(
+#     ${GTEST_INCLUDE_DIRS}
+# )
+
+# add_executable(
+#   OrthPolTEST ./OrthPolTEST.cpp
+# )
+
+# target_link_libraries(
+#   OrthPolTEST ${GTEST_LIBRARIES} Threads::Threads)
+
+# enable_testing()
+# add_test(OrthPolTEST "./OrthPolTEST")

Orthogonal Polynomials/TEST/OrthPolTEST.cpp

@@ -0,0 +1,192 @@
+#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();
+// }
+
+
+

Orthogonal Polynomials/bin/CMakeLists.txt

@@ -0,0 +1,5 @@
+add_executable(${PROJECT_NAME} main.cpp)
+
+
+target_link_libraries(${PROJECT_NAME} PRIVATE Calc)
+target_include_directories(${PROJECT_NAME} PUBLIC ${PROJECT_SOURCE_DIR})

Orthogonal Polynomials/bin/main.cpp

@@ -0,0 +1,42 @@
+#include "modules/OrthPol.h"
+#include <iomanip>
+#include <string>
+
+int main()
+{   
+    // std::cout << std::setw(2) << "n    ";
+    // for (int j=1; j <=5; j++)
+    // {
+    //     if (j == 5)
+    //     {
+    //         std::cout << std::setw(3) << std::right << std::setprecision(8) << j * 0.2 << std::endl;
+    //     }
+    //     else
+    //     {
+    //         std::cout << std::setw(3) << std::right << std::setprecision(8) << j * 0.2 << "    ";
+    //     }
+    // }
+    // std::cout << "------------------------------------------------------------------------" << std::endl;
+    // for (int i=0; i<=20; ++i)
+    // {   
+    //     std::cout << std::setw(2) << i << "    ";
+    //     for (int j=1; j <=5; ++j)
+    //     {
+
+    //         std::vector<double> res = OrthPol(1, i, j*0.2);
+    //         if (j == 5)
+    //         {    
+    //             std::cout << std::to_string(res[0]).substr(0,8) << std::endl;        
+    //             //std::cout << std::setw(10) << std::right << std::setprecision(8) << res[0] << std::endl;
+    //         }
+    //         else
+    //         {
+    //             std::cout << std::to_string(res[0]).substr(0,8) << " ";  
+    //             //std::cout << std::setw(10) << std::right << std::setprecision(8) << res[0] << " " ;
+    //         }
+    //     }
+    // }
+    for (int i = 0; i < 51; ++i) {
+        OrthPol(1, i, 0);
+    }
+}

Orthogonal Polynomials/modules/CMakeLists.txt

@@ -0,0 +1 @@
+add_library(Calc OrthPol.cpp OrthPol.h)

Orthogonal Polynomials/modules/OrthPol.cpp

@@ -0,0 +1,79 @@
+// #include <iostream> 
+// #include <string>
+// #include <vector>
+#include  "OrthPol.h"
+
+/**
+ * Purpose: Compute orthogonal polynomials: Chebyshev, Laguerre, Hermite polynomials, and their derivatives
+ * 
+ * @param KF
+ *      KF = 1 for Chebyshev polynomial of the first kind
+ *      KF = 2 for Chebyshev polynomial of the second kind
+ *      KF = 3 for Laguerre polynomial
+ *      KF = 4 for Hermite polynomial
+ * @param n - Integer order of orthogonal polynomial.
+ * @param x - Argument of orthogonal polynomial.
+ * @return Vector with values of the orthogonal polynomial (the first place) of order n with argument x, and its first derivatives (the second place).
+ * @throws None
+ */
+std::vector<double> OrthPol(const int& KF, const int& n, const double& x)
+{
+    std::vector<double> result;
+
+    double* PL = new double[n + 1];
+    double* DPL = new double[n + 1];
+    
+    double A = 2.0;
+    double B = 0.0;
+    double C = 1.0;
+    double Y0 = 1.0;
+    double Y1 = 2.0 * x;
+    double DY0 = 0.0;
+    double DY1 = 2.0;
+    PL[0] = 1.0;
+    PL[1] = 2.0 * x;
+    DPL[0] = 0.0;
+    DPL[1] = 2.0;
+
+    if (KF == 1)
+    {
+        Y1 = x;
+        DY1 = 1.0;
+        PL[1] = x;
+        DPL[1] = 1.0;
+    }
+    else if (KF == 3)
+    {
+        Y1 = 1.0 - x;
+        DY1 = -1.0;
+        PL[1] = 1.0 - x;
+        DPL[1] = -1.0;
+    }
+    
+    for (int k = 2; k <= n; k++)
+    {
+        if (KF == 3)
+        {   
+            A = -1.0 / k;
+            B = 2.0 + A;
+            C = 1.0 + A;
+        }
+        else if (KF == 4)
+        {
+            C = 2.0 * (k - 1.0);
+        }
+        double YN = (A * x + B) * Y1 - C * Y0;
+        double DYN = A * Y1 + (A * x + B) * DY1 - C * DY0;
+        PL[k] = YN;
+        DPL[k] = DYN;
+        Y0 = Y1;
+        Y1 = YN;
+        DY0 = DY1;
+        DY1 = DYN;        
+    }
+    
+    result.push_back(PL[n]);
+    result.push_back(DPL[n]);
+
+    return result;
+}

Orthogonal Polynomials/modules/OrthPol.h

@@ -0,0 +1,6 @@
+#pragma once
+#include <iostream> 
+#include <string>
+#include <vector>  
+
+std::vector<double> OrthPol(const int& KF, const int& n, const double& x);