chore: import upstream snapshot with attribution

2026-07-13 12:37:28 +08:00
commit 29cfe479ab
432 changed files with 68491 additions and 0 deletions
@@ -0,0 +1,18 @@
+# If necessary, use the RELATIVE flag, otherwise each source file may be listed
+# with full pathname. RELATIVE may makes it easier to extract an executable name
+# automatically.
+file( GLOB APP_SOURCES RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.cpp )
+# file( GLOB APP_SOURCES ${CMAKE_SOURCE_DIR}/*.c )
+# AUX_SOURCE_DIRECTORY(${CMAKE_CURRENT_SOURCE_DIR} APP_SOURCES)
+foreach( testsourcefile ${APP_SOURCES} )
+    # I used a simple string replace, to cut off .cpp.
+    string( REPLACE ".cpp" "" testname ${testsourcefile} )
+    add_executable( ${testname} ${testsourcefile} )
+
+    set_target_properties(${testname} PROPERTIES LINKER_LANGUAGE CXX)
+    if(OpenMP_CXX_FOUND)
+        target_link_libraries(${testname} OpenMP::OpenMP_CXX)
+    endif()
+    install(TARGETS ${testname} DESTINATION "bin/numerical_methods")
+
+endforeach( testsourcefile ${APP_SOURCES} )
@@ -0,0 +1,98 @@
+/**
+ * @file
+ * @brief [A babylonian method
+ * (BM)](https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Babylonian_method)
+ * is an algorithm that computes the square root.
+ * @details
+ * This algorithm has an application in use case scenario where a user wants
+ * find accurate square roots of large numbers
+ * @author [Ameya Chawla](https://github.com/ameyachawlaggsipu)
+ */
+
+#include <cassert>  /// for assert
+#include <cmath>
+#include <iostream>  /// for IO operations
+
+/**
+ * @namespace numerical_methods
+ * @brief Numerical algorithms/methods
+ */
+
+namespace numerical_methods {
+
+/**
+ * @brief Babylonian methods is an iterative function which returns
+ * square root of radicand
+ * @param radicand is the radicand
+ * @returns x1 the square root of radicand
+ */
+
+double babylonian_method(double radicand) {
+    int i = 1;  /// To find initial root or rough approximation
+
+    while (i * i <= radicand) {
+        i++;
+    }
+
+    i--;  /// Real Initial value will be i-1 as loop stops on +1 value
+
+    double x0 = i;  /// Storing previous value for comparison
+    double x1 =
+        (radicand / x0 + x0) / 2;  /// Storing calculated value for comparison
+    double temp = NAN;             /// Temp variable to x0 and x1
+
+    while (std::max(x0, x1) - std::min(x0, x1) < 0.0001) {
+        temp = (radicand / x1 + x1) / 2;  /// Newly calculated root
+        x0 = x1;
+        x1 = temp;
+    }
+
+    return x1;  /// Returning final root
+}
+
+}  // namespace numerical_methods
+
+/**
+ * @brief Self-test implementations
+ * @details
+ * Declaring two test cases and checking for the error
+ * in predicted and true value is less than 0.0001.
+ * @returns void
+ */
+static void test() {
+    /* descriptions of the following test */
+
+    auto testcase1 = 125348;  /// Testcase 1
+    auto testcase2 = 752080;  /// Testcase 2
+
+    auto real_output1 = 354.045194855;  /// Real Output 1
+    auto real_output2 = 867.225460881;  /// Real Output 2
+
+    auto test_result1 = numerical_methods::babylonian_method(testcase1);
+    /// Test result for testcase 1
+    auto test_result2 = numerical_methods::babylonian_method(testcase2);
+    /// Test result for testcase 2
+
+    assert(std::max(test_result1, real_output1) -
+               std::min(test_result1, real_output1) <
+           0.0001);
+    /// Testing for test Case 1
+    assert(std::max(test_result2, real_output2) -
+               std::min(test_result2, real_output2) <
+           0.0001);
+    /// Testing for test Case 2
+
+    std::cout << "All tests have successfully passed!\n";
+}
+
+/**
+ * @brief Main function
+ * calls automated test function to test the working of fast fourier transform.
+ * @returns 0 on exit
+ */
+
+int main() {
+    test();  //  run self-test implementations
+             //  with 2 defined test cases
+    return 0;
+}
@@ -0,0 +1,75 @@
+/**
+ * \file
+ * \brief Solve the equation \f$f(x)=0\f$ using [bisection
+ * method](https://en.wikipedia.org/wiki/Bisection_method)
+ *
+ * Given two points \f$a\f$ and \f$b\f$ such that \f$f(a)<0\f$ and
+ * \f$f(b)>0\f$, then the \f$(i+1)^\text{th}\f$ approximation is given by: \f[
+ * x_{i+1} = \frac{a_i+b_i}{2}
+ * \f]
+ * For the next iteration, the interval is selected
+ * as: \f$[a,x]\f$ if \f$x>0\f$ or \f$[x,b]\f$ if \f$x<0\f$. The Process is
+ * continued till a close enough approximation is achieved.
+ *
+ * \see newton_raphson_method.cpp, false_position.cpp, secant_method.cpp
+ */
+#include <cmath>
+#include <iostream>
+#include <limits>
+
+#define EPSILON \
+    1e-6  // std::numeric_limits<double>::epsilon()  ///< system accuracy limit
+#define MAX_ITERATIONS 50000  ///< Maximum number of iterations to check
+
+/** define \f$f(x)\f$ to find root for
+ */
+static double eq(double i) {
+    return (std::pow(i, 3) - (4 * i) - 9);  // original equation
+}
+
+/** get the sign of any given number */
+template <typename T>
+int sgn(T val) {
+    return (T(0) < val) - (val < T(0));
+}
+
+/** main function */
+int main() {
+    double a = -1, b = 1, x, z;
+    int i;
+
+    // loop to find initial intervals a, b
+    for (int i = 0; i < MAX_ITERATIONS; i++) {
+        z = eq(a);
+        x = eq(b);
+        if (sgn(z) == sgn(x)) {  // same signs, increase interval
+            b++;
+            a--;
+        } else {  // if opposite signs, we got our interval
+            break;
+        }
+    }
+
+    std::cout << "\nFirst initial: " << a;
+    std::cout << "\nSecond initial: " << b;
+
+    // start iterations
+    for (i = 0; i < MAX_ITERATIONS; i++) {
+        x = (a + b) / 2;
+        z = eq(x);
+        std::cout << "\n\nz: " << z << "\t[" << a << " , " << b
+                  << " | Bisect: " << x << "]";
+
+        if (z < 0) {
+            a = x;
+        } else {
+            b = x;
+        }
+
+        if (std::abs(z) < EPSILON)  // stoping criteria
+            break;
+    }
+
+    std::cout << "\n\nRoot: " << x << "\t\tSteps: " << i << std::endl;
+    return 0;
+}
@@ -0,0 +1,216 @@
+/**
+ * \file
+ * \brief Find real extrema of a univariate real function in a given interval
+ * using [Brent's method](https://en.wikipedia.org/wiki/Brent%27s_method).
+ *
+ * Refer the algorithm discoverer's publication
+ * [online](https://maths-people.anu.edu.au/~brent/pd/rpb011i.pdf) and also
+ * associated book:
+ * > R. P. Brent, Algorithms for Minimization without
+ * > Derivatives, Prentice-Hall, Englewood Cliffs, New Jersey, 1973
+ *
+ * \see golden_search_extrema.cpp
+ *
+ * \author [Krishna Vedala](https://github.com/kvedala)
+ */
+#define _USE_MATH_DEFINES  ///< required for MS Visual C++
+#include <cassert>
+#include <cmath>
+#include <cstdint>
+#include <functional>
+#include <iostream>
+#include <limits>
+
+#define EPSILON \
+    std::sqrt(  \
+        std::numeric_limits<double>::epsilon())  ///< system accuracy limit
+
+/**
+ * @brief Get the real root of a function in the given interval.
+ *
+ * @param f function to get root for
+ * @param lim_a lower limit of search window
+ * @param lim_b upper limit of search window
+ * @return root found in the interval
+ */
+double get_minima(const std::function<double(double)> &f, double lim_a,
+                  double lim_b) {
+    uint32_t iters = 0;
+
+    if (lim_a > lim_b) {
+        std::swap(lim_a, lim_b);
+    } else if (std::abs(lim_a - lim_b) <= EPSILON) {
+        std::cerr << "Search range must be greater than " << EPSILON << "\n";
+        return lim_a;
+    }
+
+    // golden ratio value
+    const double M_GOLDEN_RATIO = (3.f - std::sqrt(5.f)) / 2.f;
+
+    double v = lim_a + M_GOLDEN_RATIO * (lim_b - lim_a);
+    double u, w = v, x = v;
+    double fu, fv = f(v);
+    double fw = fv, fx = fv;
+
+    double mid_point = (lim_a + lim_b) / 2.f;
+    double p = 0, q = 0, r = 0;
+
+    double d, e = 0;
+    double tolerance, tolerance2;
+
+    do {
+        mid_point = (lim_a + lim_b) / 2.f;
+        tolerance = EPSILON * std::abs(x);
+        tolerance2 = 2 * tolerance;
+
+        if (std::abs(e) > tolerance2) {
+            // fit parabola
+            r = (x - w) * (fx - fv);
+            q = (x - v) * (fx - fw);
+            p = (x - v) * q - (x - w) * r;
+            q = 2.f * (q - r);
+            if (q > 0)
+                p = -p;
+            else
+                q = -q;
+            r = e;
+            e = d;
+        }
+
+        if (std::abs(p) < std::abs(0.5 * q * r) && p < q * (lim_b - x)) {
+            // parabolic interpolation step
+            d = p / q;
+            u = x + d;
+            if (u - lim_a < tolerance2 || lim_b - u < tolerance2)
+                d = x < mid_point ? tolerance : -tolerance;
+        } else {
+            // golden section interpolation step
+            e = (x < mid_point ? lim_b : lim_a) - x;
+            d = M_GOLDEN_RATIO * e;
+        }
+
+        // evaluate not too close to x
+        if (std::abs(d) >= tolerance)
+            u = d;
+        else if (d > 0)
+            u = tolerance;
+        else
+            u = -tolerance;
+        u += x;
+        fu = f(u);
+
+        // update variables
+        if (fu <= fx) {
+            if (u < x)
+                lim_b = x;
+            else
+                lim_a = x;
+            v = w;
+            fv = fw;
+            w = x;
+            fw = fx;
+            x = u;
+            fx = fu;
+        } else {
+            if (u < x)
+                lim_a = u;
+            else
+                lim_b = u;
+            if (fu <= fw || x == w) {
+                v = w;
+                fv = fw;
+                w = u;
+                fw = fu;
+            } else if (fu <= fv || v == x || v == w) {
+                v = u;
+                fv = fu;
+            }
+        }
+
+        iters++;
+    } while (std::abs(x - mid_point) > (tolerance - (lim_b - lim_a) / 2.f));
+
+    std::cout << " (iters: " << iters << ") ";
+
+    return x;
+}
+
+/**
+ * @brief Test function to find root for the function
+ * \f$f(x)= (x-2)^2\f$
+ * in the interval \f$[1,5]\f$
+ * \n Expected result = 2
+ */
+void test1() {
+    // define the function to minimize as a lambda function
+    std::function<double(double)> f1 = [](double x) {
+        return (x - 2) * (x - 2);
+    };
+
+    std::cout << "Test 1.... ";
+
+    double minima = get_minima(f1, -1, 5);
+
+    std::cout << minima << "...";
+
+    assert(std::abs(minima - 2) < EPSILON);
+    std::cout << "passed\n";
+}
+
+/**
+ * @brief Test function to find root for the function
+ * \f$f(x)= x^{\frac{1}{x}}\f$
+ * in the interval \f$[-2,10]\f$
+ * \n Expected result: \f$e\approx 2.71828182845904509\f$
+ */
+void test2() {
+    // define the function to maximize as a lambda function
+    // since we are maximixing, we negated the function return value
+    std::function<double(double)> func = [](double x) {
+        return -std::pow(x, 1.f / x);
+    };
+
+    std::cout << "Test 2.... ";
+
+    double minima = get_minima(func, -2, 5);
+
+    std::cout << minima << " (" << M_E << ")...";
+
+    assert(std::abs(minima - M_E) < EPSILON);
+    std::cout << "passed\n";
+}
+
+/**
+ * @brief Test function to find *maxima* for the function
+ * \f$f(x)= \cos x\f$
+ * in the interval \f$[0,12]\f$
+ * \n Expected result: \f$\pi\approx 3.14159265358979312\f$
+ */
+void test3() {
+    // define the function to maximize as a lambda function
+    // since we are maximixing, we negated the function return value
+    std::function<double(double)> func = [](double x) { return std::cos(x); };
+
+    std::cout << "Test 3.... ";
+
+    double minima = get_minima(func, -4, 12);
+
+    std::cout << minima << " (" << M_PI << ")...";
+
+    assert(std::abs(minima - M_PI) < EPSILON);
+    std::cout << "passed\n";
+}
+
+/** Main function */
+int main() {
+    std::cout.precision(18);
+
+    std::cout << "Computations performed with machine epsilon: " << EPSILON
+              << "\n";
+
+    test1();
+    test2();
+    test3();
+
+    return 0;
+}
@@ -0,0 +1,206 @@
+/**
+ * @file
+ * @brief Implementation of the Composite Simpson Rule for the approximation
+ *
+ * @details The following is an implementation of the Composite Simpson Rule for
+ * the approximation of definite integrals. More info -> wiki:
+ * https://en.wikipedia.org/wiki/Simpson%27s_rule#Composite_Simpson's_rule
+ *
+ * The idea is to split the interval in an EVEN number N of intervals and use as
+ * interpolation points the xi for which it applies that xi = x0 + i*h, where h
+ * is a step defined as h = (b-a)/N where a and b are the first and last points
+ * of the interval of the integration [a, b].
+ *
+ * We create a table of the xi and their corresponding f(xi) values and we
+ * evaluate the integral by the formula: I = h/3 * {f(x0) + 4*f(x1) + 2*f(x2) +
+ * ... + 2*f(xN-2) + 4*f(xN-1) + f(xN)}
+ *
+ * That means that the first and last indexed i f(xi) are multiplied by 1,
+ * the odd indexed f(xi) by 4 and the even by 2.
+ *
+ * In this program there are 4 sample test functions f, g, k, l that are
+ * evaluated in the same interval.
+ *
+ * Arguments can be passed as parameters from the command line argv[1] = N,
+ * argv[2] = a, argv[3] = b
+ *
+ * N must be even number and a<b.
+ *
+ * In the end of the main() i compare the program's result with the one from
+ * mathematical software with 2 decimal points margin.
+ *
+ * Add sample function by replacing one of the f, g, k, l and the assert
+ *
+ * @author [ggkogkou](https://github.com/ggkogkou)
+ *
+ */
+
+#include <cassert>  /// for assert
+#include <cmath>    /// for math functions
+#include <cmath>
+#include <cstdint>     /// for integer allocation
+#include <cstdlib>     /// for std::atof
+#include <functional>  /// for std::function
+#include <iostream>    /// for IO operations
+#include <map>         /// for std::map container
+
+/**
+ * @namespace numerical_methods
+ * @brief Numerical algorithms/methods
+ */
+namespace numerical_methods {
+/**
+ * @namespace simpson_method
+ * @brief Contains the Simpson's method implementation
+ */
+namespace simpson_method {
+/**
+ * @fn double evaluate_by_simpson(int N, double h, double a,
+ * std::function<double (double)> func)
+ * @brief Calculate integral or assert if integral is not a number (Nan)
+ * @param N number of intervals
+ * @param h step
+ * @param a x0
+ * @param func: choose the function that will be evaluated
+ * @returns the result of the integration
+ */
+double evaluate_by_simpson(std::int32_t N, double h, double a,
+                           const std::function<double(double)>& func) {
+    std::map<std::int32_t, double>
+        data_table;  // Contains the data points. key: i, value: f(xi)
+    double xi = a;   // Initialize xi to the starting point x0 = a
+
+    // Create the data table
+    double temp = NAN;
+    for (std::int32_t i = 0; i <= N; i++) {
+        temp = func(xi);
+        data_table.insert(
+            std::pair<std::int32_t, double>(i, temp));  // add i and f(xi)
+        xi += h;  // Get the next point xi for the next iteration
+    }
+
+    // Evaluate the integral.
+    // Remember: f(x0) + 4*f(x1) + 2*f(x2) + ... + 2*f(xN-2) + 4*f(xN-1) + f(xN)
+    double evaluate_integral = 0;
+    for (std::int32_t i = 0; i <= N; i++) {
+        if (i == 0 || i == N) {
+            evaluate_integral += data_table.at(i);
+        } else if (i % 2 == 1) {
+            evaluate_integral += 4 * data_table.at(i);
+        } else {
+            evaluate_integral += 2 * data_table.at(i);
+        }
+    }
+
+    // Multiply by the coefficient h/3
+    evaluate_integral *= h / 3;
+
+    // If the result calculated is nan, then the user has given wrong input
+    // interval.
+    assert(!std::isnan(evaluate_integral) &&
+           "The definite integral can't be evaluated. Check the validity of "
+           "your input.\n");
+    // Else return
+    return evaluate_integral;
+}
+
+/**
+ * @fn double f(double x)
+ * @brief A function f(x) that will be used to test the method
+ * @param x The independent variable xi
+ * @returns the value of the dependent variable yi = f(xi)
+ */
+double f(double x) { return std::sqrt(x) + std::log(x); }
+/** @brief Another test function */
+double g(double x) { return std::exp(-x) * (4 - std::pow(x, 2)); }
+/** @brief Another test function */
+double k(double x) { return std::sqrt(2 * std::pow(x, 3) + 3); }
+/** @brief Another test function*/
+double l(double x) { return x + std::log(2 * x + 1); }
+}  // namespace simpson_method
+}  // namespace numerical_methods
+
+/**
+ * \brief Self-test implementations
+ * @param N is the number of intervals
+ * @param h is the step
+ * @param a is x0
+ * @param b is the end of the interval
+ * @param used_argv_parameters is 'true' if argv parameters are given and
+ * 'false' if not
+ */
+static void test(std::int32_t N, double h, double a, double b,
+                 bool used_argv_parameters) {
+    // Call the functions and find the integral of each function
+    double result_f = numerical_methods::simpson_method::evaluate_by_simpson(
+        N, h, a, numerical_methods::simpson_method::f);
+    assert((used_argv_parameters || (result_f >= 4.09 && result_f <= 4.10)) &&
+           "The result of f(x) is wrong");
+    std::cout << "The result of integral f(x) on interval [" << a << ", " << b
+              << "] is equal to: " << result_f << std::endl;
+
+    double result_g = numerical_methods::simpson_method::evaluate_by_simpson(
+        N, h, a, numerical_methods::simpson_method::g);
+    assert((used_argv_parameters || (result_g >= 0.27 && result_g <= 0.28)) &&
+           "The result of g(x) is wrong");
+    std::cout << "The result of integral g(x) on interval [" << a << ", " << b
+              << "] is equal to: " << result_g << std::endl;
+
+    double result_k = numerical_methods::simpson_method::evaluate_by_simpson(
+        N, h, a, numerical_methods::simpson_method::k);
+    assert((used_argv_parameters || (result_k >= 9.06 && result_k <= 9.07)) &&
+           "The result of k(x) is wrong");
+    std::cout << "The result of integral k(x) on interval [" << a << ", " << b
+              << "] is equal to: " << result_k << std::endl;
+
+    double result_l = numerical_methods::simpson_method::evaluate_by_simpson(
+        N, h, a, numerical_methods::simpson_method::l);
+    assert((used_argv_parameters || (result_l >= 7.16 && result_l <= 7.17)) &&
+           "The result of l(x) is wrong");
+    std::cout << "The result of integral l(x) on interval [" << a << ", " << b
+              << "] is equal to: " << result_l << std::endl;
+}
+
+/**
+ * @brief Main function
+ * @param argc commandline argument count
+ * @param argv commandline array of arguments
+ * @returns 0 on exit
+ */
+int main(int argc, char** argv) {
+    std::int32_t N = 16;  /// Number of intervals to divide the integration
+                          /// interval. MUST BE EVEN
+    double a = 1, b = 3;  /// Starting and ending point of the integration in
+                          /// the real axis
+    double h = NAN;       /// Step, calculated by a, b and N
+
+    bool used_argv_parameters =
+        false;  // If argv parameters are used then the assert must be omitted
+                // for the tst cases
+
+    // Get user input (by the command line parameters or the console after
+    // displaying messages)
+    if (argc == 4) {
+        N = std::atoi(argv[1]);
+        a = std::atof(argv[2]);
+        b = std::atof(argv[3]);
+        // Check if a<b else abort
+        assert(a < b && "a has to be less than b");
+        assert(N > 0 && "N has to be > 0");
+        if (N < 16 || a != 1 || b != 3) {
+            used_argv_parameters = true;
+        }
+        std::cout << "You selected N=" << N << ", a=" << a << ", b=" << b
+                  << std::endl;
+    } else {
+        std::cout << "Default N=" << N << ", a=" << a << ", b=" << b
+                  << std::endl;
+    }
+
+    // Find the step
+    h = (b - a) / N;
+
+    test(N, h, a, b, used_argv_parameters);  // run self-test implementations
+
+    return 0;
+}
@@ -0,0 +1,340 @@
+/**
+ * @file
+ * \brief Compute all possible approximate roots of any given polynomial using
+ * [Durand Kerner
+ * algorithm](https://en.wikipedia.org/wiki/Durand%E2%80%93Kerner_method)
+ * \author [Krishna Vedala](https://github.com/kvedala)
+ *
+ * Test the algorithm online:
+ * https://gist.github.com/kvedala/27f1b0b6502af935f6917673ec43bcd7
+ *
+ * Try the highly unstable Wilkinson's polynomial:
+ * ```
+ * ./numerical_methods/durand_kerner_roots 1 -210 20615 -1256850 53327946
+ * -1672280820 40171771630 -756111184500 11310276995381 -135585182899530
+ * 1307535010540395 -10142299865511450 63030812099294896 -311333643161390640
+ * 1206647803780373360 -3599979517947607200 8037811822645051776
+ * -12870931245150988800 13803759753640704000 -8752948036761600000
+ * 2432902008176640000
+ * ```
+ * Sample implementation results to compute approximate roots of the equation
+ * \f$x^4-1=0\f$:\n
+ * <img
+ * src="https://raw.githubusercontent.com/TheAlgorithms/C-Plus-Plus/docs/images/numerical_methods/durand_kerner_error.svg"
+ * width="400" alt="Error evolution during root approximations computed every
+ * iteration."/> <img
+ * src="https://raw.githubusercontent.com/TheAlgorithms/C-Plus-Plus/docs/images/numerical_methods/durand_kerner_roots.svg"
+ * width="400" alt="Roots evolution - shows the initial approximation of the
+ * roots and their convergence to a final approximation along with the iterative
+ * approximations" />
+ */
+
+#include <algorithm>
+#include <cassert>
+#include <cmath>
+#include <complex>
+#include <cstdint>
+#include <cstdlib>
+#include <ctime>
+#include <fstream>
+#include <iostream>
+#include <valarray>
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
+#define ACCURACY 1e-10 /**< maximum accuracy limit */
+
+/**
+ * Evaluate the value of a polynomial with given coefficients
+ * \param[in] coeffs coefficients of the polynomial
+ * \param[in] x point at which to evaluate the polynomial
+ * \returns \f$f(x)\f$
+ **/
+std::complex<double> poly_function(const std::valarray<double> &coeffs,
+                                   std::complex<double> x) {
+    double real = 0.f, imag = 0.f;
+    int n;
+
+    // #ifdef _OPENMP
+    // #pragma omp target teams distribute reduction(+ : real, imag)
+    // #endif
+    for (n = 0; n < coeffs.size(); n++) {
+        std::complex<double> tmp =
+            coeffs[n] * std::pow(x, coeffs.size() - n - 1);
+        real += tmp.real();
+        imag += tmp.imag();
+    }
+
+    return std::complex<double>(real, imag);
+}
+
+/**
+ * create a textual form of complex number
+ * \param[in] x point at which to evaluate the polynomial
+ * \returns pointer to converted string
+ */
+const char *complex_str(const std::complex<double> &x) {
+#define MAX_BUFF_SIZE 50
+    static char msg[MAX_BUFF_SIZE];
+
+    std::snprintf(msg, MAX_BUFF_SIZE, "% 7.04g%+7.04gj", x.real(), x.imag());
+
+    return msg;
+}
+
+/**
+ * check for termination condition
+ * \param[in] delta point at which to evaluate the polynomial
+ * \returns `false` if termination not reached
+ * \returns `true` if termination reached
+ */
+bool check_termination(long double delta) {
+    static long double past_delta = INFINITY;
+    if (std::abs(past_delta - delta) <= ACCURACY || delta < ACCURACY)
+        return true;
+    past_delta = delta;
+    return false;
+}
+
+/**
+ * Implements Durand Kerner iterative algorithm to compute all roots of a
+ * polynomial.
+ *
+ * \param[in] coeffs coefficients of the polynomial
+ * \param[out] roots the computed roots of the polynomial
+ * \param[in] write_log flag whether to save the log file (default = `false`)
+ * \returns pair of values - number of iterations taken and final accuracy
+ * achieved
+ */
+std::pair<uint32_t, double> durand_kerner_algo(
+    const std::valarray<double> &coeffs,
+    std::valarray<std::complex<double>> *roots, bool write_log = false) {
+    long double tol_condition = 1;
+    uint32_t iter = 0;
+    int n;
+    std::ofstream log_file;
+
+    if (write_log) {
+        /*
+         * store intermediate values to a CSV file
+         */
+        log_file.open("durand_kerner.log.csv");
+        if (!log_file.is_open()) {
+            perror("Unable to create a storage log file!");
+            std::exit(EXIT_FAILURE);
+        }
+        log_file << "iter#,";
+
+        for (n = 0; n < roots->size(); n++) log_file << "root_" << n << ",";
+
+        log_file << "avg. correction";
+        log_file << "\n0,";
+        for (n = 0; n < roots->size(); n++)
+            log_file << complex_str((*roots)[n]) << ",";
+    }
+
+    bool break_loop = false;
+    while (!check_termination(tol_condition) && iter < INT16_MAX &&
+           !break_loop) {
+        tol_condition = 0;
+        iter++;
+        break_loop = false;
+
+        if (log_file.is_open())
+            log_file << "\n" << iter << ",";
+
+#ifdef _OPENMP
+#pragma omp parallel for shared(break_loop, tol_condition)
+#endif
+        for (n = 0; n < roots->size(); n++) {
+            if (break_loop)
+                continue;
+
+            std::complex<double> numerator, denominator;
+            numerator = poly_function(coeffs, (*roots)[n]);
+            denominator = 1.0;
+            for (int i = 0; i < roots->size(); i++)
+                if (i != n)
+                    denominator *= (*roots)[n] - (*roots)[i];
+
+            std::complex<long double> delta = numerator / denominator;
+
+            if (std::isnan(std::abs(delta)) || std::isinf(std::abs(delta))) {
+                std::cerr << "\n\nOverflow/underrun error - got value = "
+                          << std::abs(delta) << "\n";
+                // return std::pair<uint32_t, double>(iter, tol_condition);
+                break_loop = true;
+            }
+
+            (*roots)[n] -= delta;
+
+#ifdef _OPENMP
+#pragma omp critical
+#endif
+            tol_condition = std::max(tol_condition, std::abs(std::abs(delta)));
+        }
+        // tol_condition /= (degree - 1);
+
+        if (break_loop)
+            break;
+
+        if (log_file.is_open()) {
+            for (n = 0; n < roots->size(); n++)
+                log_file << complex_str((*roots)[n]) << ",";
+        }
+
+#if defined(DEBUG) || !defined(NDEBUG)
+        if (iter % 500 == 0) {
+            std::cout << "Iter: " << iter << "\t";
+            for (n = 0; n < roots->size(); n++)
+                std::cout << "\t" << complex_str((*roots)[n]);
+            std::cout << "\t\tabsolute average change: " << tol_condition
+                      << "\n";
+        }
+#endif
+
+        if (log_file.is_open())
+            log_file << tol_condition;
+    }
+
+    return std::pair<uint32_t, long double>(iter, tol_condition);
+}
+
+/**
+ * Self test the algorithm by checking the roots for \f$x^2+4=0\f$ to which the
+ * roots are \f$0 \pm 2i\f$
+ */
+void test1() {
+    const std::valarray<double> coeffs = {1, 0, 4};  // x^2 - 2 = 0
+    std::valarray<std::complex<double>> roots(2);
+    std::valarray<std::complex<double>> expected = {
+        std::complex<double>(0., 2.),
+        std::complex<double>(0., -2.)  // known expected roots
+    };
+
+    /* initialize root approximations with random values */
+    for (int n = 0; n < roots.size(); n++) {
+        roots[n] = std::complex<double>(std::rand() % 100, std::rand() % 100);
+        roots[n] -= 50.f;
+        roots[n] /= 25.f;
+    }
+
+    auto result = durand_kerner_algo(coeffs, &roots, false);
+
+    for (int i = 0; i < roots.size(); i++) {
+        // check if approximations are have < 0.1% error with one of the
+        // expected roots
+        bool err1 = false;
+        for (int j = 0; j < roots.size(); j++)
+            err1 |= std::abs(std::abs(roots[i] - expected[j])) < 1e-3;
+        assert(err1);
+    }
+
+    std::cout << "Test 1 passed! - " << result.first << " iterations, "
+              << result.second << " accuracy"
+              << "\n";
+}
+
+/**
+ * Self test the algorithm by checking the roots for \f$0.015625x^3-1=0\f$ to
+ * which the roots are \f$(4+0i),\,(-2\pm3.464i)\f$
+ */
+void test2() {
+    const std::valarray<double> coeffs = {// 0.015625 x^3 - 1 = 0
+                                          1. / 64., 0., 0., -1.};
+    std::valarray<std::complex<double>> roots(3);
+    const std::valarray<std::complex<double>> expected = {
+        std::complex<double>(4., 0.), std::complex<double>(-2., 3.46410162),
+        std::complex<double>(-2., -3.46410162)  // known expected roots
+    };
+
+    /* initialize root approximations with random values */
+    for (int n = 0; n < roots.size(); n++) {
+        roots[n] = std::complex<double>(std::rand() % 100, std::rand() % 100);
+        roots[n] -= 50.f;
+        roots[n] /= 25.f;
+    }
+
+    auto result = durand_kerner_algo(coeffs, &roots, false);
+
+    for (int i = 0; i < roots.size(); i++) {
+        // check if approximations are have < 0.1% error with one of the
+        // expected roots
+        bool err1 = false;
+        for (int j = 0; j < roots.size(); j++)
+            err1 |= std::abs(std::abs(roots[i] - expected[j])) < 1e-3;
+        assert(err1);
+    }
+
+    std::cout << "Test 2 passed! - " << result.first << " iterations, "
+              << result.second << " accuracy"
+              << "\n";
+}
+
+/***
+ * Main function.
+ * The comandline input arguments are taken as coeffiecients of a
+ *polynomial. For example, this command
+ * ```sh
+ * ./durand_kerner_roots 1 0 -4
+ * ```
+ * will find roots of the polynomial \f$1\cdot x^2 + 0\cdot x^1 + (-4)=0\f$
+ **/
+int main(int argc, char **argv) {
+    /* initialize random seed: */
+    std::srand(std::time(nullptr));
+
+    if (argc < 2) {
+        test1();  // run tests when no input is provided
+        test2();  // and skip tests when input polynomial is provided
+        std::cout << "Please pass the coefficients of the polynomial as "
+                     "commandline "
+                     "arguments.\n";
+        return 0;
+    }
+
+    int n, degree = argc - 1;              // detected polynomial degree
+    std::valarray<double> coeffs(degree);  // create coefficiencts array
+
+    // number of roots = degree - 1
+    std::valarray<std::complex<double>> s0(degree - 1);
+
+    std::cout << "Computing the roots for:\n\t";
+    for (n = 0; n < degree; n++) {
+        coeffs[n] = strtod(argv[n + 1], nullptr);
+        if (n < degree - 1 && coeffs[n] != 0)
+            std::cout << "(" << coeffs[n] << ") x^" << degree - n - 1 << " + ";
+        else if (coeffs[n] != 0)
+            std::cout << "(" << coeffs[n] << ") x^" << degree - n - 1
+                      << " = 0\n";
+
+        /* initialize root approximations with random values */
+        if (n < degree - 1) {
+            s0[n] = std::complex<double>(std::rand() % 100, std::rand() % 100);
+            s0[n] -= 50.f;
+            s0[n] /= 50.f;
+        }
+    }
+
+    // numerical errors less when the first coefficient is "1"
+    // hence, we normalize the first coefficient
+    {
+        double tmp = coeffs[0];
+        coeffs /= tmp;
+    }
+
+    clock_t end_time, start_time = clock();
+    auto result = durand_kerner_algo(coeffs, &s0, true);
+    end_time = clock();
+
+    std::cout << "\nIterations: " << result.first << "\n";
+    for (n = 0; n < degree - 1; n++)
+        std::cout << "\t" << complex_str(s0[n]) << "\n";
+    std::cout << "absolute average change: " << result.second << "\n";
+    std::cout << "Time taken: "
+              << static_cast<double>(end_time - start_time) / CLOCKS_PER_SEC
+              << " sec\n";
+
+    return 0;
+}
@@ -0,0 +1,128 @@
+/**
+ * \file
+ * \brief Solve the equation \f$f(x)=0\f$ using [false position
+ * method](https://en.wikipedia.org/wiki/Regula_falsi), also known as the Secant
+ * method
+ *
+ * \details
+ * First, multiple intervals are selected with the interval gap provided.
+ * Separate recursive function called for every root.
+ * Roots are printed Separatelt.
+ *
+ * For an interval [a,b] \f$a\f$ and \f$b\f$ such that \f$f(a)<0\f$ and
+ * \f$f(b)>0\f$, then the \f$(i+1)^\text{th}\f$ approximation is given by: \f[
+ * x_{i+1} = \frac{a_i\cdot f(b_i) - b_i\cdot f(a_i)}{f(b_i) - f(a_i)}
+ * \f]
+ * For the next iteration, the interval is selected
+ * as: \f$[a,x]\f$ if \f$x>0\f$ or \f$[x,b]\f$ if \f$x<0\f$. The Process is
+ * continued till a close enough approximation is achieved.
+ *
+ * \see newton_raphson_method.cpp, bisection_method.cpp
+ *
+ * \author Unknown author
+ * \author [Samruddha Patil](https://github.com/sampatil578)
+ */
+#include <cmath>     /// for math operations
+#include <iostream>  /// for io operations
+
+/**
+ * @namespace numerical_methods
+ * @brief Numerical methods
+ */
+namespace numerical_methods {
+/**
+ * @namespace false_position
+ * @brief Functions for [False Position]
+ * (https://en.wikipedia.org/wiki/Regula_falsi) method.
+ */
+namespace false_position {
+/**
+ * @brief This function gives the value of f(x) for given x.
+ * @param x value for which we have to find value of f(x).
+ * @return value of f(x) for given x.
+ */
+static float eq(float x) {
+    return (x * x - x);  // original equation
+}
+
+/**
+* @brief This function finds root of the equation in given interval i.e.
+(x1,x2).
+* @param x1,x2 values for an interval in which root is present.
+  @param y1,y2 values of function at x1, x2 espectively.
+* @return root of the equation in the given interval.
+*/
+static float regula_falsi(float x1, float x2, float y1, float y2) {
+    float diff = x1 - x2;
+    if (diff < 0) {
+        diff = (-1) * diff;
+    }
+    if (diff < 0.00001) {
+        if (y1 < 0) {
+            y1 = -y1;
+        }
+        if (y2 < 0) {
+            y2 = -y2;
+        }
+        if (y1 < y2) {
+            return x1;
+        } else {
+            return x2;
+        }
+    }
+    float x3 = 0, y3 = 0;
+    x3 = x1 - (x1 - x2) * (y1) / (y1 - y2);
+    y3 = eq(x3);
+    return regula_falsi(x2, x3, y2, y3);
+}
+
+/**
+ * @brief This function prints roots of the equation.
+ * @param root which we have to print.
+ * @param count which is count of the root in an interval [-range,range].
+ */
+void printRoot(float root, const int16_t &count) {
+    if (count == 1) {
+        std::cout << "Your 1st root is : " << root << std::endl;
+    } else if (count == 2) {
+        std::cout << "Your 2nd root is : " << root << std::endl;
+    } else if (count == 3) {
+        std::cout << "Your 3rd root is : " << root << std::endl;
+    } else {
+        std::cout << "Your " << count << "th root is : " << root << std::endl;
+    }
+}
+}  // namespace false_position
+}  // namespace numerical_methods
+
+/**
+ * @brief Main function
+ * @returns 0 on exit
+ */
+int main() {
+    float a = 0, b = 0, i = 0, root = 0;
+    int16_t count = 0;
+    float range =
+        100000;       // Range in which we have to find the root. (-range,range)
+    float gap = 0.5;  // interval gap. lesser the gap more the accuracy
+    a = numerical_methods::false_position::eq((-1) * range);
+    i = ((-1) * range + gap);
+    // while loop for selecting proper interval in provided range and with
+    // provided interval gap.
+    while (i <= range) {
+        b = numerical_methods::false_position::eq(i);
+        if (b == 0) {
+            count++;
+            numerical_methods::false_position::printRoot(i, count);
+        }
+        if (a * b < 0) {
+            root = numerical_methods::false_position::regula_falsi(i - gap, i,
+                                                                   a, b);
+            count++;
+            numerical_methods::false_position::printRoot(root, count);
+        }
+        a = b;
+        i += gap;
+    }
+    return 0;
+}
@@ -0,0 +1,165 @@
+/**
+ * @file
+ * @brief [A fast Fourier transform
+ * (FFT)](https://medium.com/@aiswaryamathur/understanding-fast-fouriertransform-from-scratch-to-solve-polynomial-multiplication-8018d511162f)
+ * is an algorithm that computes the
+ * discrete Fourier transform (DFT) of a sequence, or its inverse (IDFT).
+ * @details
+ * This
+ * algorithm has application in use case scenario where a user wants to find
+ points of a
+ * function
+ * in a short time by just using the coefficients of the polynomial
+ * function.
+ * It can be also used to find inverse fourier transform by just switching the
+ value of omega.
+ * Time complexity
+ * this algorithm computes the DFT in O(nlogn) time in comparison to traditional
+ O(n^2).
+ * @author [Ameya Chawla](https://github.com/ameyachawlaggsipu)
+ */
+
+#include <cassert>   /// for assert
+#include <cmath>     /// for mathematical-related functions
+#include <complex>   /// for storing points and coefficents
+#include <cstdint>
+#include <iostream>  /// for IO operations
+#include <vector>    /// for std::vector
+
+/**
+ * @namespace numerical_methods
+ * @brief Numerical algorithms/methods
+ */
+namespace numerical_methods {
+/**
+ * @brief FastFourierTransform is a recursive function which returns list of
+ * complex numbers
+ * @param p List of Coefficents in form of complex numbers
+ * @param n Count of elements in list p
+ * @returns p if n==1
+ * @returns y if n!=1
+ */
+std::complex<double> *FastFourierTransform(std::complex<double> *p, uint8_t n) {
+    if (n == 1) {
+        return p;  /// Base Case To return
+    }
+
+    double pi = 2 * asin(1.0);  /// Declaring value of pi
+
+    std::complex<double> om = std::complex<double>(
+        cos(2 * pi / n), sin(2 * pi / n));  /// Calculating value of omega
+
+    auto *pe = new std::complex<double>[n / 2];  /// Coefficients of even power
+
+    auto *po = new std::complex<double>[n / 2];  /// Coefficients of odd power
+
+    int k1 = 0, k2 = 0;
+    for (int j = 0; j < n; j++) {
+        if (j % 2 == 0) {
+            pe[k1++] = p[j];  /// Assigning values of even Coefficients
+
+        } else {
+            po[k2++] = p[j];  /// Assigning value of odd Coefficients
+        }
+    }
+
+    std::complex<double> *ye =
+        FastFourierTransform(pe, n / 2);  /// Recursive Call
+
+    std::complex<double> *yo =
+        FastFourierTransform(po, n / 2);  /// Recursive Call
+
+    auto *y = new std::complex<double>[n];  /// Final value representation list
+
+    k1 = 0, k2 = 0;
+
+    for (int i = 0; i < n / 2; i++) {
+        y[i] =
+            ye[k1] + pow(om, i) * yo[k2];  /// Updating the first n/2 elements
+        y[i + n / 2] =
+            ye[k1] - pow(om, i) * yo[k2];  /// Updating the last n/2 elements
+
+        k1++;
+        k2++;
+    }
+
+    if (n != 2) {
+        delete[] pe;
+        delete[] po;
+    }
+
+    delete[] ye;  /// Deleting dynamic array ye
+    delete[] yo;  /// Deleting dynamic array yo
+    return y;
+}
+
+}  // namespace numerical_methods
+
+/**
+ * @brief Self-test implementations
+ * @details
+ * Declaring two test cases and checking for the error
+ * in predicted and true value is less than 0.000000000001.
+ * @returns void
+ */
+static void test() {
+    /* descriptions of the following test */
+
+    auto *t1 = new std::complex<double>[2];  /// Test case 1
+    auto *t2 = new std::complex<double>[4];  /// Test case 2
+
+    t1[0] = {1, 0};
+    t1[1] = {2, 0};
+    t2[0] = {1, 0};
+    t2[1] = {2, 0};
+    t2[2] = {3, 0};
+    t2[3] = {4, 0};
+
+    uint8_t n1 = 2;
+    uint8_t n2 = 4;
+    std::vector<std::complex<double>> r1 = {
+        {3, 0}, {-1, 0}};  /// True Answer for test case 1
+
+    std::vector<std::complex<double>> r2 = {
+        {10, 0}, {-2, -2}, {-2, 0}, {-2, 2}};  /// True Answer for test case 2
+
+    std::complex<double> *o1 = numerical_methods::FastFourierTransform(t1, n1);
+    std::complex<double> *t3 =
+        o1;  /// Temporary variable used to delete memory location of o1
+    std::complex<double> *o2 = numerical_methods::FastFourierTransform(t2, n2);
+    std::complex<double> *t4 =
+        o2;  /// Temporary variable used to delete memory location of o2
+    for (uint8_t i = 0; i < n1; i++) {
+        assert((r1[i].real() - o1->real() < 0.000000000001) &&
+               (r1[i].imag() - o1->imag() <
+                0.000000000001));  /// Comparing for both real and imaginary
+                                   /// values for test case 1
+        o1++;
+    }
+
+    for (uint8_t i = 0; i < n2; i++) {
+        assert((r2[i].real() - o2->real() < 0.000000000001) &&
+               (r2[i].imag() - o2->imag() <
+                0.000000000001));  /// Comparing for both real and imaginary
+                                   /// values for test case 2
+        o2++;
+    }
+
+    delete[] t1;
+    delete[] t2;
+    delete[] t3;
+    delete[] t4;
+    std::cout << "All tests have successfully passed!\n";
+}
+
+/**
+ * @brief Main function
+ * calls automated test function to test the working of fast fourier transform.
+ * @returns 0 on exit
+ */
+
+int main() {
+    test();  //  run self-test implementations
+             //  with 2 defined test cases
+    return 0;
+}
@@ -0,0 +1,76 @@
+/**
+ * \file
+ * \brief [Gaussian elimination
+ * method](https://en.wikipedia.org/wiki/Gaussian_elimination)
+ */
+#include <iostream>
+
+/** Main function */
+int main() {
+    int mat_size, i, j, step;
+
+    std::cout << "Matrix size: ";
+    std::cin >> mat_size;
+
+    // create a 2D matrix by dynamic memory allocation
+    double **mat = new double *[mat_size + 1], **x = new double *[mat_size];
+    for (i = 0; i <= mat_size; i++) {
+        mat[i] = new double[mat_size + 1];
+        if (i < mat_size)
+            x[i] = new double[mat_size + 1];
+    }
+
+    // get the matrix elements from user
+    std::cout << std::endl << "Enter value of the matrix: " << std::endl;
+    for (i = 0; i < mat_size; i++) {
+        for (j = 0; j <= mat_size; j++) {
+            std::cin >>
+                mat[i][j];  // Enter (mat_size*mat_size) value of the matrix.
+        }
+    }
+
+    // perform Gaussian elimination
+    for (step = 0; step < mat_size - 1; step++) {
+        for (i = step; i < mat_size - 1; i++) {
+            double a = (mat[i + 1][step] / mat[step][step]);
+
+            for (j = step; j <= mat_size; j++)
+                mat[i + 1][j] = mat[i + 1][j] - (a * mat[step][j]);
+        }
+    }
+
+    std::cout << std::endl
+              << "Matrix using Gaussian Elimination method: " << std::endl;
+    for (i = 0; i < mat_size; i++) {
+        for (j = 0; j <= mat_size; j++) {
+            x[i][j] = mat[i][j];
+            std::cout << mat[i][j] << " ";
+        }
+        std::cout << std::endl;
+    }
+    std::cout << std::endl
+              << "Value of the Gaussian Elimination method: " << std::endl;
+    for (i = mat_size - 1; i >= 0; i--) {
+        double sum = 0;
+        for (j = mat_size - 1; j > i; j--) {
+            x[i][j] = x[j][j] * x[i][j];
+            sum = x[i][j] + sum;
+        }
+        if (x[i][i] == 0)
+            x[i][i] = 0;
+        else
+            x[i][i] = (x[i][mat_size] - sum) / (x[i][i]);
+
+        std::cout << "x" << i << "= " << x[i][i] << std::endl;
+    }
+
+    for (i = 0; i <= mat_size; i++) {
+        delete[] mat[i];
+        if (i < mat_size)
+            delete[] x[i];
+    }
+    delete[] mat;
+    delete[] x;
+
+    return 0;
+}
@@ -0,0 +1,151 @@
+/**
+ * \file
+ * \brief Find extrema of a univariate real function in a given interval using
+ * [golden section search
+ * algorithm](https://en.wikipedia.org/wiki/Golden-section_search).
+ *
+ * \see brent_method_extrema.cpp
+ * \author [Krishna Vedala](https://github.com/kvedala)
+ */
+#define _USE_MATH_DEFINES  //< required for MS Visual C++
+#include <cassert>
+#include <cmath>
+#include <cstdint>
+#include <functional>
+#include <iostream>
+#include <limits>
+
+#define EPSILON 1e-7  ///< solution accuracy limit
+
+/**
+ * @brief Get the minima of a function in the given interval. To get the maxima,
+ * simply negate the function. The golden ratio used here is:\f[
+ * k=\frac{3-\sqrt{5}}{2} \approx 0.381966\ldots\f]
+ *
+ * @param f function to get minima for
+ * @param lim_a lower limit of search window
+ * @param lim_b upper limit of search window
+ * @return local minima found in the interval
+ */
+double get_minima(const std::function<double(double)> &f, double lim_a,
+                  double lim_b) {
+    uint32_t iters = 0;
+    double c, d;
+    double prev_mean, mean = std::numeric_limits<double>::infinity();
+
+    // golden ratio value
+    const double M_GOLDEN_RATIO = (1.f + std::sqrt(5.f)) / 2.f;
+
+    // ensure that lim_a < lim_b
+    if (lim_a > lim_b) {
+        std::swap(lim_a, lim_b);
+    } else if (std::abs(lim_a - lim_b) <= EPSILON) {
+        std::cerr << "Search range must be greater than " << EPSILON << "\n";
+        return lim_a;
+    }
+
+    do {
+        prev_mean = mean;
+
+        // compute the section ratio width
+        double ratio = (lim_b - lim_a) / M_GOLDEN_RATIO;
+        c = lim_b - ratio;  // right-side section start
+        d = lim_a + ratio;  // left-side section end
+
+        if (f(c) < f(d)) {
+            // select left section
+            lim_b = d;
+        } else {
+            // selct right section
+            lim_a = c;
+        }
+
+        mean = (lim_a + lim_b) / 2.f;
+        iters++;
+
+        // continue till the interval width is greater than sqrt(system epsilon)
+    } while (std::abs(lim_a - lim_b) > EPSILON);
+
+    std::cout << " (iters: " << iters << ") ";
+    return prev_mean;
+}
+
+/**
+ * @brief Test function to find minima for the function
+ * \f$f(x)= (x-2)^2\f$
+ * in the interval \f$[1,5]\f$
+ * \n Expected result = 2
+ */
+void test1() {
+    // define the function to minimize as a lambda function
+    std::function<double(double)> f1 = [](double x) {
+        return (x - 2) * (x - 2);
+    };
+
+    std::cout << "Test 1.... ";
+
+    double minima = get_minima(f1, 1, 5);
+
+    std::cout << minima << "...";
+
+    assert(std::abs(minima - 2) < EPSILON);
+    std::cout << "passed\n";
+}
+
+/**
+ * @brief Test function to find *maxima* for the function
+ * \f$f(x)= x^{\frac{1}{x}}\f$
+ * in the interval \f$[-2,10]\f$
+ * \n Expected result: \f$e\approx 2.71828182845904509\f$
+ */
+void test2() {
+    // define the function to maximize as a lambda function
+    // since we are maximixing, we negated the function return value
+    std::function<double(double)> func = [](double x) {
+        return -std::pow(x, 1.f / x);
+    };
+
+    std::cout << "Test 2.... ";
+
+    double minima = get_minima(func, -2, 10);
+
+    std::cout << minima << " (" << M_E << ")...";
+
+    assert(std::abs(minima - M_E) < EPSILON);
+    std::cout << "passed\n";
+}
+
+/**
+ * @brief Test function to find *maxima* for the function
+ * \f$f(x)= \cos x\f$
+ * in the interval \f$[0,12]\f$
+ * \n Expected result: \f$\pi\approx 3.14159265358979312\f$
+ */
+void test3() {
+    // define the function to maximize as a lambda function
+    // since we are maximixing, we negated the function return value
+    std::function<double(double)> func = [](double x) { return std::cos(x); };
+
+    std::cout << "Test 3.... ";
+
+    double minima = get_minima(func, -4, 12);
+
+    std::cout << minima << " (" << M_PI << ")...";
+
+    assert(std::abs(minima - M_PI) < EPSILON);
+    std::cout << "passed\n";
+}
+
+/** Main function */
+int main() {
+    std::cout.precision(9);
+
+    std::cout << "Computations performed with machine epsilon: " << EPSILON
+              << "\n";
+
+    test1();
+    test2();
+    test3();
+
+    return 0;
+}
@@ -0,0 +1,290 @@
+/**
+ * @file
+ * @brief [Gram Schmidt Orthogonalisation
+ * Process](https://en.wikipedia.org/wiki/Gram%E2%80%93Schmidt_process)
+ *
+ * @details
+ * Takes the input of Linearly Independent Vectors,
+ * returns vectors orthogonal to each other.
+ *
+ * ### Algorithm
+ * Take the first vector of given LI vectors as first vector of Orthogonal
+ * vectors. Take projection of second input vector on the first vector of
+ * Orthogonal vector and subtract it from the 2nd LI vector. Take projection of
+ * third vector on the second vector of Othogonal vectors and subtract it from
+ * the 3rd LI vector. Keep repeating the above process until all the vectors in
+ * the given input array are exhausted.
+ *
+ * For Example:
+ * In R2,
+ * Input LI Vectors={(3,1),(2,2)}
+ * then Orthogonal Vectors= {(3, 1),(-0.4, 1.2)}
+ *
+ *  Have defined maximum dimension of vectors to be 10 and number of vectors
+ *  taken is 20.
+ *  Please do not give linearly dependent vectors
+ *
+ *
+ * @author [Akanksha Gupta](https://github.com/Akanksha-Gupta920)
+ */
+
+#include <array>     /// for std::array
+#include <cassert>   /// for assert
+#include <cmath>     /// for fabs
+#include <iostream>  /// for io operations
+
+#include "math.h"
+
+/**
+ * @namespace numerical_methods
+ * @brief Numerical Methods algorithms
+ */
+namespace numerical_methods {
+/**
+ * @namespace gram_schmidt
+ * @brief Functions for [Gram Schmidt Orthogonalisation
+ * Process](https://en.wikipedia.org/wiki/Gram%E2%80%93Schmidt_process)
+ */
+namespace gram_schmidt {
+/**
+ * Dot product function.
+ * Takes 2 vectors along with their dimension as input and returns the dot
+ * product.
+ * @param x vector 1
+ * @param y vector 2
+ * @param c dimension of the vectors
+ *
+ * @returns sum
+ */
+double dot_product(const std::array<double, 10>& x,
+                   const std::array<double, 10>& y, const int& c) {
+    double sum = 0;
+    for (int i = 0; i < c; ++i) {
+        sum += x[i] * y[i];
+    }
+    return sum;
+}
+
+/**
+ * Projection Function
+ * Takes input of 2 vectors along with their dimension and evaluates their
+ * projection in temp
+ *
+ * @param x Vector 1
+ * @param y Vector 2
+ * @param c dimension of each vector
+ *
+ * @returns factor
+ */
+double projection(const std::array<double, 10>& x,
+                  const std::array<double, 10>& y, const int& c) {
+    double dot =
+        dot_product(x, y, c);  /// The dot product of two vectors is taken
+    double anorm =
+        dot_product(y, y, c);  /// The norm of the second vector is taken.
+    double factor =
+        dot /
+        anorm;  /// multiply that factor with every element in a 3rd vector,
+                /// whose initial values are same as the 2nd vector.
+    return factor;
+}
+
+/**
+ * Function to print the orthogonalised vector
+ *
+ * @param r number of vectors
+ * @param c dimenaion of vectors
+ * @param B stores orthogonalised vectors
+ *
+ * @returns void
+ */
+void display(const int& r, const int& c,
+             const std::array<std::array<double, 10>, 20>& B) {
+    for (int i = 0; i < r; ++i) {
+        std::cout << "Vector " << i + 1 << ": ";
+        for (int j = 0; j < c; ++j) {
+            std::cout << B[i][j] << " ";
+        }
+        std::cout << '\n';
+    }
+}
+
+/**
+ * Function for the process of Gram Schimdt Process
+ * @param r number of vectors
+ * @param c dimension of vectors
+ * @param A stores input of given LI vectors
+ * @param B stores orthogonalised vectors
+ *
+ * @returns void
+ */
+void gram_schmidt(int r, const int& c,
+                  const std::array<std::array<double, 10>, 20>& A,
+                  std::array<std::array<double, 10>, 20> B) {
+    if (c < r) {  /// we check whether appropriate dimensions are given or not.
+        std::cout << "Dimension of vector is less than number of vector, hence "
+                     "\n first "
+                  << c << " vectors are orthogonalised\n";
+        r = c;
+    }
+
+    int k = 1;
+
+    while (k <= r) {
+        if (k == 1) {
+            for (int j = 0; j < c; j++)
+                B[0][j] = A[0][j];  /// First vector is copied as it is.
+        }
+
+        else {
+            std::array<double, 10>
+                all_projection{};  /// array to store projections
+            for (int i = 0; i < c; ++i) {
+                all_projection[i] = 0;  /// First initialised to zero
+            }
+
+            int l = 1;
+            while (l < k) {
+                std::array<double, 10>
+                    temp{};           /// to store previous projected array
+                double factor = NAN;  /// to store the factor by which the
+                                      /// previous array will change
+                factor = projection(A[k - 1], B[l - 1], c);
+                for (int i = 0; i < c; ++i) {
+                    temp[i] = B[l - 1][i] * factor;  /// projected array created
+                }
+                for (int j = 0; j < c; ++j) {
+                    all_projection[j] =
+                        all_projection[j] +
+                        temp[j];  /// we take the projection with all the
+                                  /// previous vector and add them.
+                }
+                l++;
+            }
+            for (int i = 0; i < c; ++i) {
+                B[k - 1][i] =
+                    A[k - 1][i] -
+                    all_projection[i];  /// subtract total projection vector
+                                        /// from the input vector
+            }
+        }
+        k++;
+    }
+    display(r, c, B);  // for displaying orthogoanlised vectors
+}
+}  // namespace gram_schmidt
+}  // namespace numerical_methods
+/**
+ * Test Function. Process has been tested for 3 Sample Inputs
+ * @returns void
+ */
+static void test() {
+    std::array<std::array<double, 10>, 20> a1 = {
+        {{1, 0, 1, 0}, {1, 1, 1, 1}, {0, 1, 2, 1}}};
+    std::array<std::array<double, 10>, 20> b1 = {{0}};
+    double dot1 = 0;
+    numerical_methods::gram_schmidt::gram_schmidt(3, 4, a1, b1);
+    int flag = 1;
+    for (int i = 0; i < 2; ++i) {
+        for (int j = i + 1; j < 3; ++j) {
+            dot1 = fabs(
+                numerical_methods::gram_schmidt::dot_product(b1[i], b1[j], 4));
+            if (dot1 > 0.1) {
+                flag = 0;
+                break;
+            }
+        }
+    }
+    if (flag == 0)
+        std::cout << "Vectors are linearly dependent\n";
+    assert(flag == 1);
+    std::cout << "Passed Test Case 1\n ";
+
+    std::array<std::array<double, 10>, 20> a2 = {{{3, 1}, {2, 2}}};
+    std::array<std::array<double, 10>, 20> b2 = {{0}};
+    double dot2 = 0;
+    numerical_methods::gram_schmidt::gram_schmidt(2, 2, a2, b2);
+    flag = 1;
+    for (int i = 0; i < 1; ++i) {
+        for (int j = i + 1; j < 2; ++j) {
+            dot2 = fabs(
+                numerical_methods::gram_schmidt::dot_product(b2[i], b2[j], 2));
+            if (dot2 > 0.1) {
+                flag = 0;
+                break;
+            }
+        }
+    }
+    if (flag == 0)
+        std::cout << "Vectors are linearly dependent\n";
+    assert(flag == 1);
+    std::cout << "Passed Test Case 2\n";
+
+    std::array<std::array<double, 10>, 20> a3 = {{{1, 2, 2}, {-4, 3, 2}}};
+    std::array<std::array<double, 10>, 20> b3 = {{0}};
+    double dot3 = 0;
+    numerical_methods::gram_schmidt::gram_schmidt(2, 3, a3, b3);
+    flag = 1;
+    for (int i = 0; i < 1; ++i) {
+        for (int j = i + 1; j < 2; ++j) {
+            dot3 = fabs(
+                numerical_methods::gram_schmidt::dot_product(b3[i], b3[j], 3));
+            if (dot3 > 0.1) {
+                flag = 0;
+                break;
+            }
+        }
+    }
+    if (flag == 0)
+        std::cout << "Vectors are linearly dependent\n";
+    assert(flag == 1);
+    std::cout << "Passed Test Case 3\n";
+}
+
+/**
+ * @brief Main Function
+ * @return 0 on exit
+ */
+int main() {
+    int r = 0, c = 0;
+    test();  // perform self tests
+    std::cout << "Enter the dimension of your vectors\n";
+    std::cin >> c;
+    std::cout << "Enter the number of vectors you will enter\n";
+    std::cin >> r;
+
+    std::array<std::array<double, 10>, 20>
+        A{};  /// a 2-D array for storing all vectors
+    std::array<std::array<double, 10>, 20> B = {
+        {0}};  /// a 2-D array for storing orthogonalised vectors
+    /// storing vectors in array A
+    for (int i = 0; i < r; ++i) {
+        std::cout << "Enter vector " << i + 1
+                  << '\n';  /// Input of vectors is taken
+        for (int j = 0; j < c; ++j) {
+            std::cout << "Value " << j + 1 << "th of vector: ";
+            std::cin >> A[i][j];
+        }
+        std::cout << '\n';
+    }
+
+    numerical_methods::gram_schmidt::gram_schmidt(r, c, A, B);
+
+    double dot = 0;
+    int flag = 1;  /// To check whether vectors are orthogonal or  not
+    for (int i = 0; i < r - 1; ++i) {
+        for (int j = i + 1; j < r; ++j) {
+            dot = fabs(
+                numerical_methods::gram_schmidt::dot_product(B[i], B[j], c));
+            if (dot > 0.1)  /// take make the process numerically stable, upper
+                            /// bound for the dot product take 0.1
+            {
+                flag = 0;
+                break;
+            }
+        }
+    }
+    if (flag == 0)
+        std::cout << "Vectors are linearly dependent\n";
+    return 0;
+}
@@ -0,0 +1,159 @@
+/**
+ * @file
+ * @brief [An inverse fast Fourier transform
+ * (IFFT)](https://www.geeksforgeeks.org/python-inverse-fast-fourier-transformation/)
+ * is an algorithm that computes the inverse fourier transform.
+ * @details
+ * This algorithm has an application in use case scenario where a user wants
+ * find coefficients of a function in a short time by just using points
+ * generated by DFT. Time complexity this algorithm computes the IDFT in
+ * O(nlogn) time in comparison to traditional O(n^2).
+ * @author [Ameya Chawla](https://github.com/ameyachawlaggsipu)
+ */
+
+#include <cassert>   /// for assert
+#include <cmath>     /// for mathematical-related functions
+#include <complex>   /// for storing points and coefficents
+#include <cstdint>
+#include <iostream>  /// for IO operations
+#include <vector>    /// for std::vector
+
+/**
+ * @namespace numerical_methods
+ * @brief Numerical algorithms/methods
+ */
+namespace numerical_methods {
+/**
+ * @brief InverseFastFourierTransform is a recursive function which returns list
+ * of complex numbers
+ * @param p List of Coefficents in form of complex numbers
+ * @param n Count of elements in list p
+ * @returns p if n==1
+ * @returns y if n!=1
+ */
+std::complex<double> *InverseFastFourierTransform(std::complex<double> *p,
+                                                  uint8_t n) {
+    if (n == 1) {
+        return p;  /// Base Case To return
+    }
+
+    double pi = 2 * asin(1.0);  /// Declaring value of pi
+
+    std::complex<double> om = std::complex<double>(
+        cos(2 * pi / n), sin(2 * pi / n));  /// Calculating value of omega
+
+    om.real(om.real() / n);  /// One change in comparison with DFT
+    om.imag(om.imag() / n);  /// One change in comparison with DFT
+
+    auto *pe = new std::complex<double>[n / 2];  /// Coefficients of even power
+
+    auto *po = new std::complex<double>[n / 2];  /// Coefficients of odd power
+
+    int k1 = 0, k2 = 0;
+    for (int j = 0; j < n; j++) {
+        if (j % 2 == 0) {
+            pe[k1++] = p[j];  /// Assigning values of even Coefficients
+
+        } else {
+            po[k2++] = p[j];  /// Assigning value of odd Coefficients
+        }
+    }
+
+    std::complex<double> *ye =
+        InverseFastFourierTransform(pe, n / 2);  /// Recursive Call
+
+    std::complex<double> *yo =
+        InverseFastFourierTransform(po, n / 2);  /// Recursive Call
+
+    auto *y = new std::complex<double>[n];  /// Final value representation list
+
+    k1 = 0, k2 = 0;
+
+    for (int i = 0; i < n / 2; i++) {
+        y[i] =
+            ye[k1] + pow(om, i) * yo[k2];  /// Updating the first n/2 elements
+        y[i + n / 2] =
+            ye[k1] - pow(om, i) * yo[k2];  /// Updating the last n/2 elements
+
+        k1++;
+        k2++;
+    }
+
+    if (n != 2) {
+        delete[] pe;
+        delete[] po;
+    }
+
+    delete[] ye;  /// Deleting dynamic array ye
+    delete[] yo;  /// Deleting dynamic array yo
+    return y;
+}
+
+}  // namespace numerical_methods
+
+/**
+ * @brief Self-test implementations
+ * @details
+ * Declaring two test cases and checking for the error
+ * in predicted and true value is less than 0.000000000001.
+ * @returns void
+ */
+static void test() {
+    /* descriptions of the following test */
+
+    auto *t1 = new std::complex<double>[2];  /// Test case 1
+    auto *t2 = new std::complex<double>[4];  /// Test case 2
+
+    t1[0] = {3, 0};
+    t1[1] = {-1, 0};
+    t2[0] = {10, 0};
+    t2[1] = {-2, -2};
+    t2[2] = {-2, 0};
+    t2[3] = {-2, 2};
+
+    uint8_t n1 = 2;
+    uint8_t n2 = 4;
+    std::vector<std::complex<double>> r1 = {
+        {1, 0}, {2, 0}};  /// True Answer for test case 1
+
+    std::vector<std::complex<double>> r2 = {
+        {1, 0}, {2, 0}, {3, 0}, {4, 0}};  /// True Answer for test case 2
+
+    std::complex<double> *o1 =
+        numerical_methods::InverseFastFourierTransform(t1, n1);
+
+    std::complex<double> *o2 =
+        numerical_methods::InverseFastFourierTransform(t2, n2);
+
+    for (uint8_t i = 0; i < n1; i++) {
+        assert((r1[i].real() - o1[i].real() < 0.000000000001) &&
+               (r1[i].imag() - o1[i].imag() <
+                0.000000000001));  /// Comparing for both real and imaginary
+                                   /// values for test case 1
+    }
+
+    for (uint8_t i = 0; i < n2; i++) {
+        assert((r2[i].real() - o2[i].real() < 0.000000000001) &&
+               (r2[i].imag() - o2[i].imag() <
+                0.000000000001));  /// Comparing for both real and imaginary
+                                   /// values for test case 2
+    }
+
+    delete[] t1;
+    delete[] t2;
+    delete[] o1;
+    delete[] o2;
+    std::cout << "All tests have successfully passed!\n";
+}
+
+/**
+ * @brief Main function
+ * calls automated test function to test the working of fast fourier transform.
+ * @returns 0 on exit
+ */
+
+int main() {
+    test();  //  run self-test implementations
+             //  with 2 defined test cases
+    return 0;
+}
@@ -0,0 +1,90 @@
+/**
+ * \file
+ * \brief [LU decomposition](https://en.wikipedia.org/wiki/LU_decompositon) of a
+ * square matrix
+ * \author [Krishna Vedala](https://github.com/kvedala)
+ */
+#include <cassert>
+#include <ctime>
+#include <iomanip>
+#include <iostream>
+
+#include "./lu_decomposition.h"
+
+/**
+ * operator to print a matrix
+ */
+template <typename T>
+std::ostream &operator<<(std::ostream &out, matrix<T> const &v) {
+    const int width = 10;
+    const char separator = ' ';
+
+    for (size_t row = 0; row < v.size(); row++) {
+        for (size_t col = 0; col < v[row].size(); col++)
+            out << std::left << std::setw(width) << std::setfill(separator)
+                << v[row][col];
+        out << std::endl;
+    }
+
+    return out;
+}
+
+/**
+ * Test LU decomposition
+ * \todo better ways to self-check a matrix output?
+ */
+void test1() {
+    int mat_size = 3;  // default matrix size
+    const int range = 50;
+    const int range2 = range >> 1;
+
+    /* Create a square matrix with random values */
+    matrix<double> A(mat_size, std::valarray<double>(mat_size));
+    matrix<double> L(mat_size, std::valarray<double>(mat_size));  // output
+    matrix<double> U(mat_size, std::valarray<double>(mat_size));  // output
+    for (int i = 0; i < mat_size; i++) {
+        // calloc so that all valeus are '0' by default
+        for (int j = 0; j < mat_size; j++)
+            /* create random values in the limits [-range2, range-1] */
+            A[i][j] = static_cast<double>(std::rand() % range - range2);
+    }
+
+    std::clock_t start_t = std::clock();
+    lu_decomposition(A, &L, &U);
+    std::clock_t end_t = std::clock();
+    std::cout << "Time taken: "
+              << static_cast<double>(end_t - start_t) / CLOCKS_PER_SEC << "\n";
+
+    std::cout << "A = \n" << A << "\n";
+    std::cout << "L = \n" << L << "\n";
+    std::cout << "U = \n" << U << "\n";
+}
+
+/**
+ * Test determinant computation using LU decomposition
+ */
+void test2() {
+    std::cout << "Determinant test 1...";
+    matrix<int> A1({{1, 2, 3}, {4, 9, 6}, {7, 8, 9}});
+    assert(determinant_lu(A1) == -48);
+    std::cout << "passed\n";
+
+    std::cout << "Determinant test 2...";
+    matrix<int> A2({{1, 2, 3}, {4, 5, 6}, {7, 8, 9}});
+    assert(determinant_lu(A2) == 0);
+    std::cout << "passed\n";
+
+    std::cout << "Determinant test 3...";
+    matrix<float> A3({{1.2, 2.3, 3.4}, {4.5, 5.6, 6.7}, {7.8, 8.9, 9.0}});
+    assert(determinant_lu(A3) == 3.63);
+    std::cout << "passed\n";
+}
+
+/** Main function */
+int main() {
+    std::srand(std::time(NULL));  // random number initializer
+
+    test1();
+    test2();
+    return 0;
+}
@@ -0,0 +1,102 @@
+/**
+ * @file lu_decomposition.h
+ * @author [Krishna Vedala](https://github.com/kvedala)
+ * @brief Functions associated with [LU
+ * Decomposition](https://en.wikipedia.org/wiki/LU_decomposition)
+ * of a square matrix.
+ */
+#pragma once
+
+#include <iostream>
+#include <valarray>
+#include <vector>
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
+/** Define matrix type as a `std::vector` of `std::valarray` */
+template <typename T>
+using matrix = std::vector<std::valarray<T>>;
+
+/** Perform LU decomposition on matrix
+ * \param[in] A matrix to decompose
+ * \param[out] L output L matrix
+ * \param[out] U output U matrix
+ * \returns 0 if no errors
+ * \returns negative if error occurred
+ */
+template <typename T>
+int lu_decomposition(const matrix<T> &A, matrix<double> *L, matrix<double> *U) {
+    int row, col, j;
+    int mat_size = A.size();
+
+    if (mat_size != A[0].size()) {
+        // check matrix is a square matrix
+        std::cerr << "Not a square matrix!\n";
+        return -1;
+    }
+
+    // regularize each row
+    for (row = 0; row < mat_size; row++) {
+        // Upper triangular matrix
+#ifdef _OPENMP
+#pragma omp for
+#endif
+        for (col = row; col < mat_size; col++) {
+            // Summation of L[i,j] * U[j,k]
+            double lu_sum = 0.;
+            for (j = 0; j < row; j++) {
+                lu_sum += L[0][row][j] * U[0][j][col];
+            }
+
+            // Evaluate U[i,k]
+            U[0][row][col] = A[row][col] - lu_sum;
+        }
+
+        // Lower triangular matrix
+#ifdef _OPENMP
+#pragma omp for
+#endif
+        for (col = row; col < mat_size; col++) {
+            if (row == col) {
+                L[0][row][col] = 1.;
+                continue;
+            }
+
+            // Summation of L[i,j] * U[j,k]
+            double lu_sum = 0.;
+            for (j = 0; j < row; j++) {
+                lu_sum += L[0][col][j] * U[0][j][row];
+            }
+
+            // Evaluate U[i,k]
+            L[0][col][row] = (A[col][row] - lu_sum) / U[0][row][row];
+        }
+    }
+
+    return 0;
+}
+
+/**
+ * Compute determinant of an NxN square matrix using LU decomposition.
+ * Using LU decomposition, the determinant is given by the product of diagonal
+ * elements of matrices L and U.
+ *
+ * @tparam T datatype of input matrix - int, unsigned int, double, etc
+ * @param A input square matrix
+ * @return determinant of matrix A
+ */
+template <typename T>
+double determinant_lu(const matrix<T> &A) {
+    matrix<double> L(A.size(), std::valarray<double>(A.size()));
+    matrix<double> U(A.size(), std::valarray<double>(A.size()));
+
+    if (lu_decomposition(A, &L, &U) < 0)
+        return 0;
+
+    double result = 1.f;
+    for (size_t i = 0; i < A.size(); i++) {
+        result *= L[i][i] * U[i][i];
+    }
+    return result;
+}
@@ -0,0 +1,199 @@
+/**
+ * @file
+ * @brief A numerical method for easy [approximation of
+ * integrals](https://en.wikipedia.org/wiki/Midpoint_method)
+ * @details The idea is to split the interval into N of intervals and use as
+ * interpolation points the xi for which it applies that xi = x0 + i*h, where h
+ * is a step defined as h = (b-a)/N where a and b are the first and last points
+ * of the interval of the integration [a, b].
+ *
+ * We create a table of the xi and their corresponding f(xi) values and we
+ * evaluate the integral by the formula: I = h * {f(x0+h/2) + f(x1+h/2) + ... +
+ * f(xN-1+h/2)}
+ *
+ * Arguments can be passed as parameters from the command line argv[1] = N,
+ * argv[2] = a, argv[3] = b. In this case if the default values N=16, a=1, b=3
+ * are changed then the tests/assert are disabled.
+ *
+ *
+ * @author [ggkogkou](https://github.com/ggkogkou)
+ */
+#include <cassert>     /// for assert
+#include <cmath>       /// for math functions
+#include <cstdint>     /// for integer allocation
+#include <cstdlib>     /// for std::atof
+#include <functional>  /// for std::function
+#include <iostream>    /// for IO operations
+#include <map>         /// for std::map container
+
+/**
+ * @namespace numerical_methods
+ * @brief Numerical algorithms/methods
+ */
+namespace numerical_methods {
+/**
+ * @namespace midpoint_rule
+ * @brief Functions for the [Midpoint
+ * Integral](https://en.wikipedia.org/wiki/Midpoint_method) method
+ * implementation
+ */
+namespace midpoint_rule {
+/**
+ * @fn double midpoint(const std::int32_t N, const double h, const double a,
+ * const std::function<double (double)>& func)
+ * @brief Main function for implementing the Midpoint Integral Method
+ * implementation
+ * @param N is the number of intervals
+ * @param h is the step
+ * @param a is x0
+ * @param func is the function that will be integrated
+ * @returns the result of the integration
+ */
+double midpoint(const std::int32_t N, const double h, const double a,
+                const std::function<double(double)>& func) {
+    std::map<int, double>
+        data_table;  // Contains the data points, key: i, value: f(xi)
+    double xi = a;   // Initialize xi to the starting point x0 = a
+
+    // Create the data table
+    // Loop from x0 to xN-1
+    double temp = NAN;
+    for (std::int32_t i = 0; i < N; i++) {
+        temp = func(xi + h / 2);  // find f(xi+h/2)
+        data_table.insert(
+            std::pair<std::int32_t, double>(i, temp));  // add i and f(xi)
+        xi += h;  // Get the next point xi for the next iteration
+    }
+
+    // Evaluate the integral.
+    // Remember: {f(x0+h/2) + f(x1+h/2) + ... + f(xN-1+h/2)}
+    double evaluate_integral = 0;
+    for (std::int32_t i = 0; i < N; i++) evaluate_integral += data_table.at(i);
+
+    // Multiply by the coefficient h
+    evaluate_integral *= h;
+
+    // If the result calculated is nan, then the user has given wrong input
+    // interval.
+    assert(!std::isnan(evaluate_integral) &&
+           "The definite integral can't be evaluated. Check the validity of "
+           "your input.\n");
+    // Else return
+    return evaluate_integral;
+}
+
+/**
+ * @brief A function f(x) that will be used to test the method
+ * @param x The independent variable xi
+ * @returns the value of the dependent variable yi = f(xi) = sqrt(xi) + ln(xi)
+ */
+double f(double x) { return std::sqrt(x) + std::log(x); }
+/**
+ * @brief A function g(x) that will be used to test the method
+ * @param x The independent variable xi
+ * @returns the value of the dependent variable yi = g(xi) = e^(-xi) * (4 -
+ * xi^2)
+ */
+double g(double x) { return std::exp(-x) * (4 - std::pow(x, 2)); }
+/**
+ * @brief A function k(x) that will be used to test the method
+ * @param x The independent variable xi
+ * @returns the value of the dependent variable yi = k(xi) = sqrt(2*xi^3 + 3)
+ */
+double k(double x) { return std::sqrt(2 * std::pow(x, 3) + 3); }
+/**
+ * @brief A function l(x) that will be used to test the method
+ * @param x The independent variable xi
+ * @returns the value of the dependent variable yi = l(xi) = xi + ln(2*xi + 1)
+ */
+double l(double x) { return x + std::log(2 * x + 1); }
+
+}  // namespace midpoint_rule
+}  // namespace numerical_methods
+
+/**
+ * @brief Self-test implementations
+ * @param N is the number of intervals
+ * @param h is the step
+ * @param a is x0
+ * @param b is the end of the interval
+ * @param used_argv_parameters is 'true' if argv parameters are given and
+ * 'false' if not
+ */
+static void test(std::int32_t N, double h, double a, double b,
+                 bool used_argv_parameters) {
+    // Call midpoint() for each of the test functions f, g, k, l
+    // Assert with two decimal point precision
+    double result_f = numerical_methods::midpoint_rule::midpoint(
+        N, h, a, numerical_methods::midpoint_rule::f);
+    assert((used_argv_parameters || (result_f >= 4.09 && result_f <= 4.10)) &&
+           "The result of f(x) is wrong");
+    std::cout << "The result of integral f(x) on interval [" << a << ", " << b
+              << "] is equal to: " << result_f << std::endl;
+
+    double result_g = numerical_methods::midpoint_rule::midpoint(
+        N, h, a, numerical_methods::midpoint_rule::g);
+    assert((used_argv_parameters || (result_g >= 0.27 && result_g <= 0.28)) &&
+           "The result of g(x) is wrong");
+    std::cout << "The result of integral g(x) on interval [" << a << ", " << b
+              << "] is equal to: " << result_g << std::endl;
+
+    double result_k = numerical_methods::midpoint_rule::midpoint(
+        N, h, a, numerical_methods::midpoint_rule::k);
+    assert((used_argv_parameters || (result_k >= 9.06 && result_k <= 9.07)) &&
+           "The result of k(x) is wrong");
+    std::cout << "The result of integral k(x) on interval [" << a << ", " << b
+              << "] is equal to: " << result_k << std::endl;
+
+    double result_l = numerical_methods::midpoint_rule::midpoint(
+        N, h, a, numerical_methods::midpoint_rule::l);
+    assert((used_argv_parameters || (result_l >= 7.16 && result_l <= 7.17)) &&
+           "The result of l(x) is wrong");
+    std::cout << "The result of integral l(x) on interval [" << a << ", " << b
+              << "] is equal to: " << result_l << std::endl;
+}
+
+/**
+ * @brief Main function
+ * @param argc commandline argument count
+ * @param argv commandline array of arguments
+ * @returns 0 on exit
+ */
+int main(int argc, char** argv) {
+    std::int32_t N =
+        16;  /// Number of intervals to divide the integration interval.
+    /// MUST BE EVEN
+    double a = 1, b = 3;  /// Starting and ending point of the integration in
+    /// the real axis
+    double h = NAN;  /// Step, calculated by a, b and N
+
+    bool used_argv_parameters =
+        false;  // If argv parameters are used then the assert must be omitted
+    // for the test cases
+
+    // Get user input (by the command line parameters or the console after
+    // displaying messages)
+    if (argc == 4) {
+        N = std::atoi(argv[1]);
+        a = std::atof(argv[2]);
+        b = std::atof(argv[3]);
+        // Check if a<b else abort
+        assert(a < b && "a has to be less than b");
+        assert(N > 0 && "N has to be > 0");
+        if (N < 4 || a != 1 || b != 3) {
+            used_argv_parameters = true;
+        }
+        std::cout << "You selected N=" << N << ", a=" << a << ", b=" << b
+                  << std::endl;
+    } else {
+        std::cout << "Default N=" << N << ", a=" << a << ", b=" << b
+                  << std::endl;
+    }
+
+    // Find the step
+    h = (b - a) / N;
+
+    test(N, h, a, b, used_argv_parameters);  // run self-test implementations
+
+    return 0;
+}
@@ -0,0 +1,68 @@
+/**
+ * \file
+ * \brief Solve the equation \f$f(x)=0\f$ using [Newton-Raphson
+ * method](https://en.wikipedia.org/wiki/Newton%27s_method) for both real and
+ * complex solutions
+ *
+ * The \f$(i+1)^\text{th}\f$ approximation is given by:
+ * \f[
+ * x_{i+1} = x_i - \frac{f(x_i)}{f'(x_i)}
+ * \f]
+ *
+ * \author [Krishna Vedala](https://github.com/kvedala)
+ * \see bisection_method.cpp, false_position.cpp
+ */
+#include <cmath>
+#include <cstdint>
+#include <ctime>
+#include <iostream>
+#include <limits>
+
+constexpr double EPSILON = 1e-10;              ///< system accuracy limit
+constexpr int16_t MAX_ITERATIONS = INT16_MAX;  ///< Maximum number of iterations
+
+/** define \f$f(x)\f$ to find root for.
+ * Currently defined as:
+ * \f[
+ * f(x) = x^3 - 4x - 9
+ * \f]
+ */
+static double eq(double i) {
+    return (std::pow(i, 3) - (4 * i) - 9);  // original equation
+}
+
+/** define the derivative function \f$f'(x)\f$
+ * For the current problem, it is:
+ * \f[
+ * f'(x) = 3x^2 - 4
+ * \f]
+ */
+static double eq_der(double i) {
+    return ((3 * std::pow(i, 2)) - 4);  // derivative of equation
+}
+
+/** Main function */
+int main() {
+    std::srand(std::time(nullptr));  // initialize randomizer
+
+    double z = NAN, c = std::rand() % 100, m = NAN, n = NAN;
+    int i = 0;
+
+    std::cout << "\nInitial approximation: " << c;
+
+    // start iterations
+    for (i = 0; i < MAX_ITERATIONS; i++) {
+        m = eq(c);
+        n = eq_der(c);
+
+        z = c - (m / n);
+        c = z;
+
+        if (std::abs(m) < EPSILON) {  // stoping criteria
+            break;
+        }
+    }
+
+    std::cout << "\n\nRoot: " << z << "\t\tSteps: " << i << std::endl;
+    return 0;
+}
@@ -0,0 +1,211 @@
+/**
+ * \file
+ * \authors [Krishna Vedala](https://github.com/kvedala)
+ * \brief Solve a multivariable first order [ordinary differential equation
+ * (ODEs)](https://en.wikipedia.org/wiki/Ordinary_differential_equation) using
+ * [forward Euler
+ * method](https://en.wikipedia.org/wiki/Numerical_methods_for_ordinary_differential_equations#Euler_method)
+ *
+ * \details
+ * The ODE being solved is:
+ * \f{eqnarray*}{
+ * \dot{u} &=& v\\
+ * \dot{v} &=& -\omega^2 u\\
+ * \omega &=& 1\\
+ * [x_0, u_0, v_0] &=& [0,1,0]\qquad\ldots\text{(initial values)}
+ * \f}
+ * The exact solution for the above problem is:
+ * \f{eqnarray*}{
+ * u(x) &=& \cos(x)\\
+ * v(x) &=& -\sin(x)\\
+ * \f}
+ * The computation results are stored to a text file `forward_euler.csv` and the
+ * exact soltuion results in `exact.csv` for comparison.
+ * <img
+ * src="https://raw.githubusercontent.com/TheAlgorithms/C-Plus-Plus/docs/images/numerical_methods/ode_forward_euler.svg"
+ * alt="Implementation solution"/>
+ *
+ * To implement [Van der Pol
+ * oscillator](https://en.wikipedia.org/wiki/Van_der_Pol_oscillator), change the
+ * ::problem function to:
+ * ```cpp
+ * const double mu = 2.0;
+ * dy[0] = y[1];
+ * dy[1] = mu * (1.f - y[0] * y[0]) * y[1] - y[0];
+ * ```
+ * \see ode_midpoint_euler.cpp, ode_semi_implicit_euler.cpp
+ */
+
+#include <cmath>
+#include <ctime>
+#include <fstream>
+#include <iostream>
+#include <valarray>
+
+/**
+ * @brief Problem statement for a system with first-order differential
+ * equations. Updates the system differential variables.
+ * \note This function can be updated to and ode of any order.
+ *
+ * @param[in] 		x 		independent variable(s)
+ * @param[in,out]	y		dependent variable(s)
+ * @param[in,out]	dy	    first-derivative of dependent variable(s)
+ */
+void problem(const double &x, std::valarray<double> *y,
+             std::valarray<double> *dy) {
+    const double omega = 1.F;             // some const for the problem
+    (*dy)[0] = (*y)[1];                   // x dot // NOLINT
+    (*dy)[1] = -omega * omega * (*y)[0];  // y dot // NOLINT
+}
+
+/**
+ * @brief Exact solution of the problem. Used for solution comparison.
+ *
+ * @param[in] 		x 		independent variable
+ * @param[in,out]	y		dependent variable
+ */
+void exact_solution(const double &x, std::valarray<double> *y) {
+    y[0][0] = std::cos(x);
+    y[0][1] = -std::sin(x);
+}
+
+/** \addtogroup ode Ordinary Differential Equations
+ * Integration functions for implementations with solving [ordinary differential
+ * equations](https://en.wikipedia.org/wiki/Ordinary_differential_equation)
+ * (ODEs) of any order and and any number of independent variables.
+ * @{
+ */
+/**
+ * @brief Compute next step approximation using the forward-Euler
+ * method. @f[y_{n+1}=y_n + dx\cdot f\left(x_n,y_n\right)@f]
+ * @param[in] 		dx	step size
+ * @param[in] 	    x	take \f$x_n\f$ and compute \f$x_{n+1}\f$
+ * @param[in,out] 	y	take \f$y_n\f$ and compute \f$y_{n+1}\f$
+ * @param[in,out]	dy	compute \f$f\left(x_n,y_n\right)\f$
+ */
+void forward_euler_step(const double dx, const double x,
+                        std::valarray<double> *y, std::valarray<double> *dy) {
+    problem(x, y, dy);
+    *y += *dy * dx;
+}
+
+/**
+ * @brief Compute approximation using the forward-Euler
+ * method in the given limits.
+ * @param[in] 		dx  	step size
+ * @param[in]   	x0  	initial value of independent variable
+ * @param[in] 	    x_max	final value of independent variable
+ * @param[in,out] 	y	    take \f$y_n\f$ and compute \f$y_{n+1}\f$
+ * @param[in] save_to_file	flag to save results to a CSV file (1) or not (0)
+ * @returns time taken for computation in seconds
+ */
+double forward_euler(double dx, double x0, double x_max,
+                     std::valarray<double> *y, bool save_to_file = false) {
+    std::valarray<double> dy = *y;
+
+    std::ofstream fp;
+    if (save_to_file) {
+        fp.open("forward_euler.csv", std::ofstream::out);
+        if (!fp.is_open()) {
+            std::perror("Error! ");
+        }
+    }
+
+    std::size_t L = y->size();
+
+    /* start integration */
+    std::clock_t t1 = std::clock();
+    double x = x0;
+
+    do {  // iterate for each step of independent variable
+        if (save_to_file && fp.is_open()) {
+            // write to file
+            fp << x << ",";
+            for (int i = 0; i < L - 1; i++) {
+                fp << y[0][i] << ",";  // NOLINT
+            }
+            fp << y[0][L - 1] << "\n";  // NOLINT
+        }
+
+        forward_euler_step(dx, x, y, &dy);  // perform integration
+        x += dx;                            // update step
+    } while (x <= x_max);  // till upper limit of independent variable
+    /* end of integration */
+    std::clock_t t2 = std::clock();
+
+    if (fp.is_open()) {
+        fp.close();
+    }
+
+    return static_cast<double>(t2 - t1) / CLOCKS_PER_SEC;
+}
+
+/** @} */
+
+/**
+ * Function to compute and save exact solution for comparison
+ *
+ * \param [in]    X0  	    initial value of independent variable
+ * \param [in] 	  X_MAX	    final value of independent variable
+ * \param [in] 	  step_size	independent variable step size
+ * \param [in]    Y0	    initial values of dependent variables
+ */
+void save_exact_solution(const double &X0, const double &X_MAX,
+                         const double &step_size,
+                         const std::valarray<double> &Y0) {
+    double x = X0;
+    std::valarray<double> y(Y0);
+
+    std::ofstream fp("exact.csv", std::ostream::out);
+    if (!fp.is_open()) {
+        std::perror("Error! ");
+        return;
+    }
+    std::cout << "Finding exact solution\n";
+
+    std::clock_t t1 = std::clock();
+    do {
+        fp << x << ",";
+        for (int i = 0; i < y.size() - 1; i++) {
+            fp << y[i] << ",";  // NOLINT
+        }
+        fp << y[y.size() - 1] << "\n";  // NOLINT
+
+        exact_solution(x, &y);
+
+        x += step_size;
+    } while (x <= X_MAX);
+
+    std::clock_t t2 = std::clock();
+    double total_time = static_cast<double>(t2 - t1) / CLOCKS_PER_SEC;
+    std::cout << "\tTime = " << total_time << " ms\n";
+
+    fp.close();
+}
+
+/**
+ * Main Function
+ */
+int main(int argc, char *argv[]) {
+    double X0 = 0.f;                    /* initial value of x0 */
+    double X_MAX = 10.F;                /* upper limit of integration */
+    std::valarray<double> Y0{1.f, 0.f}; /* initial value Y = y(x = x_0) */
+    double step_size = NAN;
+
+    if (argc == 1) {
+        std::cout << "\nEnter the step size: ";
+        std::cin >> step_size;
+    } else {
+        // use commandline argument as independent variable step size
+        step_size = std::atof(argv[1]);
+    }
+
+    // get approximate solution
+    double total_time = forward_euler(step_size, X0, X_MAX, &Y0, true);
+    std::cout << "\tTime = " << total_time << " ms\n";
+
+    /* compute exact solution for comparion */
+    save_exact_solution(X0, X_MAX, step_size, Y0);
+
+    return 0;
+}
@@ -0,0 +1,214 @@
+/**
+ * \file
+ * \authors [Krishna Vedala](https://github.com/kvedala)
+ * \brief Solve a multivariable first order [ordinary differential equation
+ * (ODEs)](https://en.wikipedia.org/wiki/Ordinary_differential_equation) using
+ * [midpoint Euler
+ * method](https://en.wikipedia.org/wiki/Midpoint_method)
+ *
+ * \details
+ * The ODE being solved is:
+ * \f{eqnarray*}{
+ * \dot{u} &=& v\\
+ * \dot{v} &=& -\omega^2 u\\
+ * \omega &=& 1\\
+ * [x_0, u_0, v_0] &=& [0,1,0]\qquad\ldots\text{(initial values)}
+ * \f}
+ * The exact solution for the above problem is:
+ * \f{eqnarray*}{
+ * u(x) &=& \cos(x)\\
+ * v(x) &=& -\sin(x)\\
+ * \f}
+ * The computation results are stored to a text file `midpoint_euler.csv` and
+ * the exact soltuion results in `exact.csv` for comparison. <img
+ * src="https://raw.githubusercontent.com/TheAlgorithms/C-Plus-Plus/docs/images/numerical_methods/ode_midpoint_euler.svg"
+ * alt="Implementation solution"/>
+ *
+ * To implement [Van der Pol
+ * oscillator](https://en.wikipedia.org/wiki/Van_der_Pol_oscillator), change the
+ * ::problem function to:
+ * ```cpp
+ * const double mu = 2.0;
+ * dy[0] = y[1];
+ * dy[1] = mu * (1.f - y[0] * y[0]) * y[1] - y[0];
+ * ```
+ * \see ode_forward_euler.cpp, ode_semi_implicit_euler.cpp
+ */
+
+#include <cmath>
+#include <ctime>
+#include <fstream>
+#include <iostream>
+#include <valarray>
+
+/**
+ * @brief Problem statement for a system with first-order differential
+ * equations. Updates the system differential variables.
+ * \note This function can be updated to and ode of any order.
+ *
+ * @param[in] 		x 		independent variable(s)
+ * @param[in,out]	y		dependent variable(s)
+ * @param[in,out]	dy	    first-derivative of dependent variable(s)
+ */
+void problem(const double &x, std::valarray<double> *y,
+             std::valarray<double> *dy) {
+    const double omega = 1.F;             // some const for the problem
+    dy[0][0] = y[0][1];                   // x dot
+    dy[0][1] = -omega * omega * y[0][0];  // y dot
+}
+
+/**
+ * @brief Exact solution of the problem. Used for solution comparison.
+ *
+ * @param[in] 		x 		independent variable
+ * @param[in,out]	y		dependent variable
+ */
+void exact_solution(const double &x, std::valarray<double> *y) {
+    y[0][0] = std::cos(x);
+    y[0][1] = -std::sin(x);
+}
+
+/** \addtogroup ode Ordinary Differential Equations
+ * @{
+ */
+/**
+ * @brief Compute next step approximation using the midpoint-Euler
+ * method.
+ * @f[y_{n+1} = y_n + dx\, f\left(x_n+\frac{1}{2}dx,
+ * y_n + \frac{1}{2}dx\,f\left(x_n,y_n\right)\right)@f]
+ *
+ * @param[in] 		dx	step size
+ * @param[in] 	    x	take \f$x_n\f$ and compute \f$x_{n+1}\f$
+ * @param[in,out] 	y	take \f$y_n\f$ and compute \f$y_{n+1}\f$
+ * @param[in,out]	dy	compute \f$f\left(x_n,y_n\right)\f$
+ */
+void midpoint_euler_step(const double dx, const double &x,
+                         std::valarray<double> *y, std::valarray<double> *dy) {
+    problem(x, y, dy);
+    double tmp_x = x + 0.5 * dx;
+
+    std::valarray<double> tmp_y = y[0] + dy[0] * (0.5 * dx);
+
+    problem(tmp_x, &tmp_y, dy);
+
+    y[0] += dy[0] * dx;
+}
+
+/**
+ * @brief Compute approximation using the midpoint-Euler
+ * method in the given limits.
+ * @param[in] 		dx  	step size
+ * @param[in]   	x0  	initial value of independent variable
+ * @param[in] 	    x_max	final value of independent variable
+ * @param[in,out] 	y	    take \f$y_n\f$ and compute \f$y_{n+1}\f$
+ * @param[in] save_to_file	flag to save results to a CSV file (1) or not (0)
+ * @returns time taken for computation in seconds
+ */
+double midpoint_euler(double dx, double x0, double x_max,
+                      std::valarray<double> *y, bool save_to_file = false) {
+    std::valarray<double> dy = y[0];
+
+    std::ofstream fp;
+    if (save_to_file) {
+        fp.open("midpoint_euler.csv", std::ofstream::out);
+        if (!fp.is_open()) {
+            std::perror("Error! ");
+        }
+    }
+
+    std::size_t L = y->size();
+
+    /* start integration */
+    std::clock_t t1 = std::clock();
+    double x = x0;
+    do {  // iterate for each step of independent variable
+        if (save_to_file && fp.is_open()) {
+            // write to file
+            fp << x << ",";
+            for (int i = 0; i < L - 1; i++) {
+                fp << y[0][i] << ",";
+            }
+            fp << y[0][L - 1] << "\n";
+        }
+
+        midpoint_euler_step(dx, x, y, &dy);  // perform integration
+        x += dx;                             // update step
+    } while (x <= x_max);  // till upper limit of independent variable
+    /* end of integration */
+    std::clock_t t2 = std::clock();
+
+    if (fp.is_open())
+        fp.close();
+
+    return static_cast<double>(t2 - t1) / CLOCKS_PER_SEC;
+}
+
+/** @} */
+
+/**
+ * Function to compute and save exact solution for comparison
+ *
+ * \param [in]    X0  	    initial value of independent variable
+ * \param [in] 	  X_MAX	    final value of independent variable
+ * \param [in] 	  step_size	independent variable step size
+ * \param [in]    Y0	    initial values of dependent variables
+ */
+void save_exact_solution(const double &X0, const double &X_MAX,
+                         const double &step_size,
+                         const std::valarray<double> &Y0) {
+    double x = X0;
+    std::valarray<double> y = Y0;
+
+    std::ofstream fp("exact.csv", std::ostream::out);
+    if (!fp.is_open()) {
+        std::perror("Error! ");
+        return;
+    }
+    std::cout << "Finding exact solution\n";
+
+    std::clock_t t1 = std::clock();
+    do {
+        fp << x << ",";
+        for (int i = 0; i < y.size() - 1; i++) {
+            fp << y[i] << ",";
+        }
+        fp << y[y.size() - 1] << "\n";
+
+        exact_solution(x, &y);
+
+        x += step_size;
+    } while (x <= X_MAX);
+
+    std::clock_t t2 = std::clock();
+    double total_time = static_cast<double>(t2 - t1) / CLOCKS_PER_SEC;
+    std::cout << "\tTime = " << total_time << " ms\n";
+
+    fp.close();
+}
+
+/**
+ * Main Function
+ */
+int main(int argc, char *argv[]) {
+    double X0 = 0.f;                       /* initial value of x0 */
+    double X_MAX = 10.F;                   /* upper limit of integration */
+    std::valarray<double> Y0 = {1.f, 0.f}; /* initial value Y = y(x = x_0) */
+    double step_size;
+
+    if (argc == 1) {
+        std::cout << "\nEnter the step size: ";
+        std::cin >> step_size;
+    } else {
+        // use commandline argument as independent variable step size
+        step_size = std::atof(argv[1]);
+    }
+
+    // get approximate solution
+    double total_time = midpoint_euler(step_size, X0, X_MAX, &Y0, true);
+    std::cout << "\tTime = " << total_time << " ms\n";
+
+    /* compute exact solution for comparion */
+    save_exact_solution(X0, X_MAX, step_size, Y0);
+
+    return 0;
+}
@@ -0,0 +1,211 @@
+/**
+ * \file
+ * \authors [Krishna Vedala](https://github.com/kvedala)
+ * \brief Solve a multivariable first order [ordinary differential equation
+ * (ODEs)](https://en.wikipedia.org/wiki/Ordinary_differential_equation) using
+ * [semi implicit Euler
+ * method](https://en.wikipedia.org/wiki/Semi-implicit_Euler_method)
+ *
+ * \details
+ * The ODE being solved is:
+ * \f{eqnarray*}{
+ * \dot{u} &=& v\\
+ * \dot{v} &=& -\omega^2 u\\
+ * \omega &=& 1\\
+ * [x_0, u_0, v_0] &=& [0,1,0]\qquad\ldots\text{(initial values)}
+ * \f}
+ * The exact solution for the above problem is:
+ * \f{eqnarray*}{
+ * u(x) &=& \cos(x)\\
+ * v(x) &=& -\sin(x)\\
+ * \f}
+ * The computation results are stored to a text file `semi_implicit_euler.csv`
+ * and the exact soltuion results in `exact.csv` for comparison. <img
+ * src="https://raw.githubusercontent.com/TheAlgorithms/C-Plus-Plus/docs/images/numerical_methods/ode_semi_implicit_euler.svg"
+ * alt="Implementation solution"/>
+ *
+ * To implement [Van der Pol
+ * oscillator](https://en.wikipedia.org/wiki/Van_der_Pol_oscillator), change the
+ * ::problem function to:
+ * ```cpp
+ * const double mu = 2.0;
+ * dy[0] = y[1];
+ * dy[1] = mu * (1.f - y[0] * y[0]) * y[1] - y[0];
+ * ```
+ * \see ode_midpoint_euler.cpp, ode_forward_euler.cpp
+ */
+
+#include <cmath>
+#include <ctime>
+#include <fstream>
+#include <iostream>
+#include <valarray>
+
+/**
+ * @brief Problem statement for a system with first-order differential
+ * equations. Updates the system differential variables.
+ * \note This function can be updated to and ode of any order.
+ *
+ * @param[in] 		x 		independent variable(s)
+ * @param[in,out]	y		dependent variable(s)
+ * @param[in,out]	dy	    first-derivative of dependent variable(s)
+ */
+void problem(const double &x, std::valarray<double> *y,
+             std::valarray<double> *dy) {
+    const double omega = 1.F;             // some const for the problem
+    dy[0][0] = y[0][1];                   // x dot
+    dy[0][1] = -omega * omega * y[0][0];  // y dot
+}
+
+/**
+ * @brief Exact solution of the problem. Used for solution comparison.
+ *
+ * @param[in] 		x 		independent variable
+ * @param[in,out]	y		dependent variable
+ */
+void exact_solution(const double &x, std::valarray<double> *y) {
+    y[0][0] = std::cos(x);
+    y[0][1] = -std::sin(x);
+}
+
+/** \addtogroup ode Ordinary Differential Equations
+ * @{
+ */
+/**
+ * @brief Compute next step approximation using the semi-implicit-Euler
+ * method. @f[y_{n+1}=y_n + dx\cdot f\left(x_n,y_n\right)@f]
+ * @param[in] 		dx	step size
+ * @param[in] 	    x	take \f$x_n\f$ and compute \f$x_{n+1}\f$
+ * @param[in,out] 	y	take \f$y_n\f$ and compute \f$y_{n+1}\f$
+ * @param[in,out]	dy	compute \f$f\left(x_n,y_n\right)\f$
+ */
+void semi_implicit_euler_step(const double dx, const double &x,
+                              std::valarray<double> *y,
+                              std::valarray<double> *dy) {
+    problem(x, y, dy);         // update dy once
+    y[0][0] += dx * dy[0][0];  // update y0
+    problem(x, y, dy);         // update dy once more
+
+    dy[0][0] = 0.f;      // ignore y0
+    y[0] += dy[0] * dx;  // update remaining using new dy
+}
+
+/**
+ * @brief Compute approximation using the semi-implicit-Euler
+ * method in the given limits.
+ * @param[in] 		dx  	step size
+ * @param[in]   	x0  	initial value of independent variable
+ * @param[in] 	    x_max	final value of independent variable
+ * @param[in,out] 	y	    take \f$y_n\f$ and compute \f$y_{n+1}\f$
+ * @param[in] save_to_file	flag to save results to a CSV file (1) or not (0)
+ * @returns time taken for computation in seconds
+ */
+double semi_implicit_euler(double dx, double x0, double x_max,
+                           std::valarray<double> *y,
+                           bool save_to_file = false) {
+    std::valarray<double> dy = y[0];
+
+    std::ofstream fp;
+    if (save_to_file) {
+        fp.open("semi_implicit_euler.csv", std::ofstream::out);
+        if (!fp.is_open()) {
+            std::perror("Error! ");
+        }
+    }
+
+    std::size_t L = y->size();
+
+    /* start integration */
+    std::clock_t t1 = std::clock();
+    double x = x0;
+    do {  // iterate for each step of independent variable
+        if (save_to_file && fp.is_open()) {
+            // write to file
+            fp << x << ",";
+            for (int i = 0; i < L - 1; i++) {
+                fp << y[0][i] << ",";
+            }
+            fp << y[0][L - 1] << "\n";
+        }
+
+        semi_implicit_euler_step(dx, x, y, &dy);  // perform integration
+        x += dx;                                  // update step
+    } while (x <= x_max);  // till upper limit of independent variable
+    /* end of integration */
+    std::clock_t t2 = std::clock();
+
+    if (fp.is_open())
+        fp.close();
+
+    return static_cast<double>(t2 - t1) / CLOCKS_PER_SEC;
+}
+
+/** @} */
+
+/**
+ * Function to compute and save exact solution for comparison
+ *
+ * \param [in]    X0  	    initial value of independent variable
+ * \param [in] 	  X_MAX	    final value of independent variable
+ * \param [in] 	  step_size	independent variable step size
+ * \param [in]    Y0	    initial values of dependent variables
+ */
+void save_exact_solution(const double &X0, const double &X_MAX,
+                         const double &step_size,
+                         const std::valarray<double> &Y0) {
+    double x = X0;
+    std::valarray<double> y = Y0;
+
+    std::ofstream fp("exact.csv", std::ostream::out);
+    if (!fp.is_open()) {
+        std::perror("Error! ");
+        return;
+    }
+    std::cout << "Finding exact solution\n";
+
+    std::clock_t t1 = std::clock();
+    do {
+        fp << x << ",";
+        for (int i = 0; i < y.size() - 1; i++) {
+            fp << y[i] << ",";
+        }
+        fp << y[y.size() - 1] << "\n";
+
+        exact_solution(x, &y);
+
+        x += step_size;
+    } while (x <= X_MAX);
+
+    std::clock_t t2 = std::clock();
+    double total_time = static_cast<double>(t2 - t1) / CLOCKS_PER_SEC;
+    std::cout << "\tTime = " << total_time << " ms\n";
+
+    fp.close();
+}
+
+/**
+ * Main Function
+ */
+int main(int argc, char *argv[]) {
+    double X0 = 0.f;                       /* initial value of x0 */
+    double X_MAX = 10.F;                   /* upper limit of integration */
+    std::valarray<double> Y0 = {1.f, 0.f}; /* initial value Y = y(x = x_0) */
+    double step_size;
+
+    if (argc == 1) {
+        std::cout << "\nEnter the step size: ";
+        std::cin >> step_size;
+    } else {
+        // use commandline argument as independent variable step size
+        step_size = std::atof(argv[1]);
+    }
+
+    // get approximate solution
+    double total_time = semi_implicit_euler(step_size, X0, X_MAX, &Y0, true);
+    std::cout << "\tTime = " << total_time << " ms\n";
+
+    /* compute exact solution for comparion */
+    save_exact_solution(X0, X_MAX, step_size, Y0);
+
+    return 0;
+}
@@ -0,0 +1,210 @@
+/**
+ * @file
+ * \brief Library functions to compute [QR
+ * decomposition](https://en.wikipedia.org/wiki/QR_decomposition) of a given
+ * matrix.
+ * \author [Krishna Vedala](https://github.com/kvedala)
+ */
+
+#ifndef NUMERICAL_METHODS_QR_DECOMPOSE_H_
+#define NUMERICAL_METHODS_QR_DECOMPOSE_H_
+
+#include <cmath>
+#include <cstdlib>
+#include <iomanip>
+#include <iostream>
+#include <limits>
+#include <numeric>
+#include <valarray>
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
+/** \namespace qr_algorithm
+ * \brief Functions to compute [QR
+ * decomposition](https://en.wikipedia.org/wiki/QR_decomposition) of any
+ * rectangular matrix
+ */
+namespace qr_algorithm {
+/**
+ * operator to print a matrix
+ */
+template <typename T>
+std::ostream &operator<<(std::ostream &out,
+                         std::valarray<std::valarray<T>> const &v) {
+    const int width = 12;
+    const char separator = ' ';
+
+    out.precision(4);
+    for (size_t row = 0; row < v.size(); row++) {
+        for (size_t col = 0; col < v[row].size(); col++)
+            out << std::right << std::setw(width) << std::setfill(separator)
+                << v[row][col];
+        out << std::endl;
+    }
+
+    return out;
+}
+
+/**
+ * operator to print a vector
+ */
+template <typename T>
+std::ostream &operator<<(std::ostream &out, std::valarray<T> const &v) {
+    const int width = 10;
+    const char separator = ' ';
+
+    out.precision(4);
+    for (size_t row = 0; row < v.size(); row++) {
+        out << std::right << std::setw(width) << std::setfill(separator)
+            << v[row];
+    }
+
+    return out;
+}
+
+/**
+ * Compute dot product of two vectors of equal lengths
+ *
+ * If \f$\vec{a}=\left[a_0,a_1,a_2,...,a_L\right]\f$ and
+ * \f$\vec{b}=\left[b_0,b_1,b_1,...,b_L\right]\f$ then
+ * \f$\vec{a}\cdot\vec{b}=\displaystyle\sum_{i=0}^L a_i\times b_i\f$
+ *
+ * \returns \f$\vec{a}\cdot\vec{b}\f$
+ */
+template <typename T>
+inline double vector_dot(const std::valarray<T> &a, const std::valarray<T> &b) {
+    return (a * b).sum();
+    // could also use following
+    // return std::inner_product(std::begin(a), std::end(a), std::begin(b),
+    // 0.f);
+}
+
+/**
+ * Compute magnitude of vector.
+ *
+ * If \f$\vec{a}=\left[a_0,a_1,a_2,...,a_L\right]\f$ then
+ * \f$\left|\vec{a}\right|=\sqrt{\displaystyle\sum_{i=0}^L a_i^2}\f$
+ *
+ * \returns \f$\left|\vec{a}\right|\f$
+ */
+template <typename T>
+inline double vector_mag(const std::valarray<T> &a) {
+    double dot = vector_dot(a, a);
+    return std::sqrt(dot);
+}
+
+/**
+ * Compute projection of vector \f$\vec{a}\f$ on \f$\vec{b}\f$ defined as
+ * \f[\text{proj}_\vec{b}\vec{a}=\frac{\vec{a}\cdot\vec{b}}{\left|\vec{b}\right|^2}\vec{b}\f]
+ *
+ * \returns NULL if error, otherwise pointer to output
+ */
+template <typename T>
+std::valarray<T> vector_proj(const std::valarray<T> &a,
+                             const std::valarray<T> &b) {
+    double num = vector_dot(a, b);
+    double deno = vector_dot(b, b);
+
+    /*! check for division by zero using machine epsilon */
+    if (deno <= std::numeric_limits<double>::epsilon()) {
+        std::cerr << "[" << __func__ << "] Possible division by zero\n";
+        return a;  // return vector a back
+    }
+
+    double scalar = num / deno;
+
+    return b * scalar;
+}
+
+/**
+ * Decompose matrix \f$A\f$ using [Gram-Schmidt
+ *process](https://en.wikipedia.org/wiki/QR_decomposition).
+ *
+ * \f{eqnarray*}{
+ * \text{given that}\quad A &=&
+ *\left[\mathbf{a}_1,\mathbf{a}_2,\ldots,\mathbf{a}_{N-1},\right]\\
+ * \text{where}\quad\mathbf{a}_i &=&
+ * \left[a_{0i},a_{1i},a_{2i},\ldots,a_{(M-1)i}\right]^T\quad\ldots\mbox{(column
+ * vectors)}\\
+ * \text{then}\quad\mathbf{u}_i &=& \mathbf{a}_i
+ *-\sum_{j=0}^{i-1}\text{proj}_{\mathbf{u}_j}\mathbf{a}_i\\
+ * \mathbf{e}_i &=&\frac{\mathbf{u}_i}{\left|\mathbf{u}_i\right|}\\
+ * Q &=& \begin{bmatrix}\mathbf{e}_0 & \mathbf{e}_1 & \mathbf{e}_2 & \dots &
+ * \mathbf{e}_{N-1}\end{bmatrix}\\
+ * R &=& \begin{bmatrix}\langle\mathbf{e}_0\,,\mathbf{a}_0\rangle &
+ * \langle\mathbf{e}_1\,,\mathbf{a}_1\rangle &
+ * \langle\mathbf{e}_2\,,\mathbf{a}_2\rangle & \dots \\
+ *                  0 & \langle\mathbf{e}_1\,,\mathbf{a}_1\rangle &
+ * \langle\mathbf{e}_2\,,\mathbf{a}_2\rangle & \dots\\
+ *                  0 & 0 & \langle\mathbf{e}_2\,,\mathbf{a}_2\rangle &
+ * \dots\\ \vdots & \vdots & \vdots & \ddots
+ *      \end{bmatrix}\\
+ * \f}
+ */
+template <typename T>
+void qr_decompose(
+    const std::valarray<std::valarray<T>> &A, /**< input matrix to decompose */
+    std::valarray<std::valarray<T>> *Q,       /**< output decomposed matrix */
+    std::valarray<std::valarray<T>> *R        /**< output decomposed matrix */
+) {
+    std::size_t ROWS = A.size();        // number of rows of A
+    std::size_t COLUMNS = A[0].size();  // number of columns of A
+    std::valarray<T> col_vector(ROWS);
+    std::valarray<T> col_vector2(ROWS);
+    std::valarray<T> tmp_vector(ROWS);
+
+    for (int i = 0; i < COLUMNS; i++) {
+        /* for each column => R is a square matrix of NxN */
+        int j;
+        R[0][i] = 0.; /* make R upper triangular */
+
+        /* get corresponding Q vector */
+#ifdef _OPENMP
+// parallelize on threads
+#pragma omp for
+#endif
+        for (j = 0; j < ROWS; j++) {
+            tmp_vector[j] = A[j][i]; /* accumulator for uk */
+            col_vector[j] = A[j][i];
+        }
+        for (j = 0; j < i; j++) {
+            for (int k = 0; k < ROWS; k++) {
+                col_vector2[k] = Q[0][k][j];
+            }
+            col_vector2 = vector_proj(col_vector, col_vector2);
+            tmp_vector -= col_vector2;
+        }
+
+        double mag = vector_mag(tmp_vector);
+
+#ifdef _OPENMP
+// parallelize on threads
+#pragma omp for
+#endif
+        for (j = 0; j < ROWS; j++) Q[0][j][i] = tmp_vector[j] / mag;
+
+            /* compute upper triangular values of R */
+#ifdef _OPENMP
+// parallelize on threads
+#pragma omp for
+#endif
+        for (int kk = 0; kk < ROWS; kk++) {
+            col_vector[kk] = Q[0][kk][i];
+        }
+
+#ifdef _OPENMP
+// parallelize on threads
+#pragma omp for
+#endif
+        for (int k = i; k < COLUMNS; k++) {
+            for (int kk = 0; kk < ROWS; kk++) {
+                col_vector2[kk] = A[kk][k];
+            }
+            R[0][i][k] = (col_vector * col_vector2).sum();
+        }
+    }
+}
+}  // namespace qr_algorithm
+
+#endif  // NUMERICAL_METHODS_QR_DECOMPOSE_H_
@@ -0,0 +1,58 @@
+/**
+ * @file
+ * \brief Program to compute the [QR
+ * decomposition](https://en.wikipedia.org/wiki/QR_decomposition) of a given
+ * matrix.
+ * \author [Krishna Vedala](https://github.com/kvedala)
+ */
+
+#include <array>
+#include <cmath>
+#include <cstdlib>
+#include <ctime>
+#include <iostream>
+
+#include "./qr_decompose.h"
+
+using qr_algorithm::qr_decompose;
+using qr_algorithm::operator<<;
+
+/**
+ * main function
+ */
+int main(void) {
+    unsigned int ROWS, COLUMNS;
+
+    std::cout << "Enter the number of rows and columns: ";
+    std::cin >> ROWS >> COLUMNS;
+
+    std::cout << "Enter matrix elements row-wise:\n";
+
+    std::valarray<std::valarray<double>> A(ROWS);
+    std::valarray<std::valarray<double>> Q(ROWS);
+    std::valarray<std::valarray<double>> R(COLUMNS);
+    for (int i = 0; i < std::max(ROWS, COLUMNS); i++) {
+        if (i < ROWS) {
+            A[i] = std::valarray<double>(COLUMNS);
+            Q[i] = std::valarray<double>(COLUMNS);
+        }
+        if (i < COLUMNS) {
+            R[i] = std::valarray<double>(COLUMNS);
+        }
+    }
+
+    for (int i = 0; i < ROWS; i++)
+        for (int j = 0; j < COLUMNS; j++) std::cin >> A[i][j];
+
+    std::cout << A << "\n";
+
+    clock_t t1 = clock();
+    qr_decompose(A, &Q, &R);
+    double dtime = static_cast<double>(clock() - t1) / CLOCKS_PER_SEC;
+
+    std::cout << Q << "\n";
+    std::cout << R << "\n";
+    std::cout << "Time taken to compute: " << dtime << " sec\n ";
+
+    return 0;
+}
@@ -0,0 +1,284 @@
+/**
+ * @file
+ * \brief Compute real eigen values and eigen vectors of a symmetric matrix
+ * using [QR decomposition](https://en.wikipedia.org/wiki/QR_decomposition)
+ * method.
+ * \author [Krishna Vedala](https://github.com/kvedala)
+ */
+#include <cassert>
+#include <cmath>
+#include <cstdlib>
+#include <ctime>
+#include <iostream>
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
+#include "./qr_decompose.h"
+using qr_algorithm::operator<<;
+
+#define LIMS 9 /**< limit of range of matrix values */
+
+/**
+ * create a symmetric square matrix of given size with random elements. A
+ * symmetric square matrix will *always* have real eigen values.
+ *
+ * \param[out] A matrix to create (must be pre-allocated in memory)
+ */
+void create_matrix(std::valarray<std::valarray<double>> *A) {
+    int i, j, tmp, lim2 = LIMS >> 1;
+    int N = A->size();
+
+#ifdef _OPENMP
+#pragma omp for
+#endif
+    for (i = 0; i < N; i++) {
+        A[0][i][i] = (std::rand() % LIMS) - lim2;
+        for (j = i + 1; j < N; j++) {
+            tmp = (std::rand() % LIMS) - lim2;
+            A[0][i][j] = tmp;  // summetrically distribute random values
+            A[0][j][i] = tmp;
+        }
+    }
+}
+
+/**
+ * Perform multiplication of two matrices.
+ * * R2 must be equal to C1
+ * * Resultant matrix size should be R1xC2
+ * \param[in] A first matrix to multiply
+ * \param[in] B second matrix to multiply
+ * \param[out] OUT output matrix (must be pre-allocated)
+ * \returns pointer to resultant matrix
+ */
+void mat_mul(const std::valarray<std::valarray<double>> &A,
+             const std::valarray<std::valarray<double>> &B,
+             std::valarray<std::valarray<double>> *OUT) {
+    int R1 = A.size();
+    int C1 = A[0].size();
+    int R2 = B.size();
+    int C2 = B[0].size();
+    if (C1 != R2) {
+        perror("Matrix dimensions mismatch!");
+        return;
+    }
+
+    for (int i = 0; i < R1; i++) {
+        for (int j = 0; j < C2; j++) {
+            OUT[0][i][j] = 0.f;
+            for (int k = 0; k < C1; k++) {
+                OUT[0][i][j] += A[i][k] * B[k][j];
+            }
+        }
+    }
+}
+
+namespace qr_algorithm {
+/** Compute eigen values using iterative shifted QR decomposition algorithm as
+ * follows:
+ * 1. Use last diagonal element of A as eigen value approximation \f$c\f$
+ * 2. Shift diagonals of matrix \f$A' = A - cI\f$
+ * 3. Decompose matrix \f$A'=QR\f$
+ * 4. Compute next approximation \f$A'_1 = RQ \f$
+ * 5. Shift diagonals back \f$A_1 = A'_1 + cI\f$
+ * 6. Termination condition check: last element below diagonal is almost 0
+ *   1. If not 0, go back to step 1 with the new approximation \f$A_1\f$
+ *   2. If 0, continue to step 7
+ * 7. Save last known \f$c\f$ as the eigen value.
+ * 8. Are all eigen values found?
+ *   1. If not, remove last row and column of \f$A_1\f$ and go back to step 1.
+ *   2. If yes, stop.
+ *
+ * \note The matrix \f$A\f$ gets modified
+ *
+ * \param[in,out] A matrix to compute eigen values for
+ * \param[in] print_intermediates (optional) whether to print intermediate A, Q
+ * and R matrices (default = `false`)
+ */
+std::valarray<double> eigen_values(std::valarray<std::valarray<double>> *A,
+                                   bool print_intermediates = false) {
+    int rows = A->size();
+    int columns = rows;
+    int counter = 0, num_eigs = rows - 1;
+    double last_eig = 0;
+
+    std::valarray<std::valarray<double>> Q(rows);
+    std::valarray<std::valarray<double>> R(columns);
+
+    /* number of eigen values = matrix size */
+    std::valarray<double> eigen_vals(rows);
+    for (int i = 0; i < rows; i++) {
+        Q[i] = std::valarray<double>(columns);
+        R[i] = std::valarray<double>(columns);
+    }
+
+    /* continue till all eigen values are found */
+    while (num_eigs > 0) {
+        /* iterate with QR decomposition */
+        while (std::abs(A[0][num_eigs][num_eigs - 1]) >
+               std::numeric_limits<double>::epsilon()) {
+            // initial approximation = last diagonal element
+            last_eig = A[0][num_eigs][num_eigs];
+            for (int i = 0; i < rows; i++) {
+                A[0][i][i] -= last_eig; /* A - cI */
+            }
+
+            qr_decompose(*A, &Q, &R);
+
+            if (print_intermediates) {
+                std::cout << *A << "\n";
+                std::cout << Q << "\n";
+                std::cout << R << "\n";
+                printf("-------------------- %d ---------------------\n",
+                       ++counter);
+            }
+
+            // new approximation A' = R * Q
+            mat_mul(R, Q, A);
+
+            for (int i = 0; i < rows; i++) {
+                A[0][i][i] += last_eig; /* A + cI */
+            }
+        }
+
+        /* store the converged eigen value */
+        eigen_vals[num_eigs] = last_eig;
+        // A[0][num_eigs][num_eigs];
+        if (print_intermediates) {
+            std::cout << "========================\n";
+            std::cout << "Eigen value: " << last_eig << ",\n";
+            std::cout << "========================\n";
+        }
+
+        num_eigs--;
+        rows--;
+        columns--;
+    }
+    eigen_vals[0] = A[0][0][0];
+
+    if (print_intermediates) {
+        std::cout << Q << "\n";
+        std::cout << R << "\n";
+    }
+
+    return eigen_vals;
+}
+
+}  // namespace qr_algorithm
+
+/**
+ * test function to compute eigen values of a 2x2 matrix
+ * \f[\begin{bmatrix}
+ * 5 & 7\\
+ * 7 & 11
+ * \end{bmatrix}\f]
+ * which are approximately, {15.56158, 0.384227}
+ */
+void test1() {
+    std::valarray<std::valarray<double>> X = {{5, 7}, {7, 11}};
+    double y[] = {15.56158, 0.384227};  // corresponding y-values
+
+    std::cout << "------- Test 1 -------" << std::endl;
+    std::valarray<double> eig_vals = qr_algorithm::eigen_values(&X);
+
+    for (int i = 0; i < 2; i++) {
+        std::cout << i + 1 << "/2 Checking for " << y[i] << " --> ";
+        bool result = false;
+        for (int j = 0; j < 2 && !result; j++) {
+            if (std::abs(y[i] - eig_vals[j]) < 0.1) {
+                result = true;
+                std::cout << "(" << eig_vals[j] << ") ";
+            }
+        }
+        assert(result);  // ensure that i^th expected eigen value was computed
+        std::cout << "found\n";
+    }
+    std::cout << "Test 1 Passed\n\n";
+}
+
+/**
+ * test function to compute eigen values of a 2x2 matrix
+ * \f[\begin{bmatrix}
+ * -4& 4& 2& 0& -3\\
+ * 4& -4& 4& -3& -1\\
+ * 2& 4& 4& 3& -3\\
+ * 0& -3& 3& -1&-1\\
+ * -3& -1& -3& -3& 0
+ * \end{bmatrix}\f]
+ * which are approximately, {9.27648, -9.26948, 2.0181, -1.03516, -5.98994}
+ */
+void test2() {
+    std::valarray<std::valarray<double>> X = {{-4, 4, 2, 0, -3},
+                                              {4, -4, 4, -3, -1},
+                                              {2, 4, 4, 3, -3},
+                                              {0, -3, 3, -1, -3},
+                                              {-3, -1, -3, -3, 0}};
+    double y[] = {9.27648, -9.26948, 2.0181, -1.03516,
+                  -5.98994};  // corresponding y-values
+
+    std::cout << "------- Test 2 -------" << std::endl;
+    std::valarray<double> eig_vals = qr_algorithm::eigen_values(&X);
+
+    std::cout << X << "\n"
+              << "Eigen values: " << eig_vals << "\n";
+
+    for (int i = 0; i < 5; i++) {
+        std::cout << i + 1 << "/5 Checking for " << y[i] << " --> ";
+        bool result = false;
+        for (int j = 0; j < 5 && !result; j++) {
+            if (std::abs(y[i] - eig_vals[j]) < 0.1) {
+                result = true;
+                std::cout << "(" << eig_vals[j] << ") ";
+            }
+        }
+        assert(result);  // ensure that i^th expected eigen value was computed
+        std::cout << "found\n";
+    }
+    std::cout << "Test 2 Passed\n\n";
+}
+
+/**
+ * main function
+ */
+int main(int argc, char **argv) {
+    int mat_size = 5;
+    if (argc == 2) {
+        mat_size = atoi(argv[1]);
+    } else {  // if invalid input argument is given run tests
+        test1();
+        test2();
+        std::cout << "Usage: ./qr_eigen_values [mat_size]\n";
+        return 0;
+    }
+
+    if (mat_size < 2) {
+        fprintf(stderr, "Matrix size should be > 2\n");
+        return -1;
+    }
+
+    // initialize random number generator
+    std::srand(std::time(nullptr));
+
+    int i, rows = mat_size, columns = mat_size;
+
+    std::valarray<std::valarray<double>> A(rows);
+
+    for (int i = 0; i < rows; i++) {
+        A[i] = std::valarray<double>(columns);
+    }
+
+    /* create a random matrix */
+    create_matrix(&A);
+
+    std::cout << A << "\n";
+
+    clock_t t1 = clock();
+    std::valarray<double> eigen_vals = qr_algorithm::eigen_values(&A);
+    double dtime = static_cast<double>(clock() - t1) / CLOCKS_PER_SEC;
+
+    std::cout << "Eigen vals: ";
+    for (i = 0; i < mat_size; i++) std::cout << eigen_vals[i] << "\t";
+    std::cout << "\nTime taken to compute: " << dtime << " sec\n";
+
+    return 0;
+}
@@ -0,0 +1,133 @@
+/**
+ * @{
+ * \file
+ * \brief [Runge Kutta fourth
+ * order](https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods) method
+ * implementation
+ *
+ * \author [Rudra Prasad Das](http://github.com/rudra697)
+ *
+ * \details
+ * It solves the unknown value of y
+ * for a given value of x
+ * only first order differential equations
+ * can be solved
+ * \example
+ * it solves \frac{\mathrm{d} y}{\mathrm{d} x}= \frac{\left ( x-y \right )}{2}
+ * given x for given initial
+ * conditions
+ * There can be many such equations
+ */
+#include <cassert>   /// asserting the test functions
+#include <cstdint>
+#include <iostream>  /// for io operations
+#include <vector>    /// for using the vector container
+/**
+ * @brief The change() function is used
+ * to return the updated iterative value corresponding
+ * to the given function
+ * @param x is the value corresponding to the x coordinate
+ * @param y is the value corresponding to the y coordinate
+ * @returns the computed function value at that call
+ */
+static double change(double x, double y) { return ((x - y) / 2.0); }
+
+/**
+ * @namespace numerical_methods
+ * @brief Numerical Methods
+ */
+namespace numerical_methods {
+/**
+ * @namespace runge_kutta
+ * @brief Functions for [Runge Kutta fourth
+ * order](https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods) method
+ */
+namespace runge_kutta {
+/**
+ * @brief the Runge Kutta method finds the value of integration of a function in
+ * the given limits. the lower limit of integration as the initial value and the
+ * upper limit is the given x
+ * @param init_x is the value of initial x and is updated after each call
+ * @param init_y is the value of initial x and is updated after each call
+ * @param x is current iteration at which the function needs to be evaluated
+ * @param h is the step value
+ * @returns the value of y at thr required value of x from the initial
+ * conditions
+ */
+double rungeKutta(double init_x, const double &init_y, const double &x,
+                  const double &h) {
+    // Count number of iterations
+    // using step size or
+    // step height h
+
+    // n calucates the number of iterations
+    // k1, k2, k3, k4 are the Runge Kutta variables
+    // used for calculation of y at each iteration
+
+    auto n = static_cast<uint64_t>((x - init_x) / h);
+    // used a vector container for the variables
+    std::vector<double> k(4, 0.0);
+
+    // Iterate for number of iterations
+
+    double y = init_y;
+    for (int i = 1; i <= n; ++i) {
+        // Apply Runge Kutta Formulas
+        // to find next value of y
+        k[0] = h * change(init_x, y);
+        k[1] = h * change(init_x + 0.5 * h, y + 0.5 * k[0]);
+        k[2] = h * change(init_x + 0.5 * h, y + 0.5 * k[1]);
+        k[3] = h * change(init_x + h, y + k[2]);
+
+        // Update next value of y
+
+        y += (1.0 / 6.0) * (k[0] + 2 * k[1] + 2 * k[2] + k[3]);
+
+        // Update next value of x
+
+        init_x += h;
+    }
+
+    return y;
+}
+}  // namespace runge_kutta
+}  // namespace numerical_methods
+
+/**
+ * @brief Tests to check algorithm implementation.
+ * @returns void
+ */
+static void test() {
+    std::cout << "The Runge Kutta function will be tested on the basis of "
+                 "precomputed values\n";
+
+    std::cout << "Test 1...."
+              << "\n";
+    double valfirst = numerical_methods::runge_kutta::rungeKutta(
+        2, 3, 4, 0.2);  // Tests the function with pre calculated values
+    assert(valfirst == 3.10363932323749570);
+    std::cout << "Passed Test 1\n";
+
+    std::cout << "Test 2...."
+              << "\n";
+    double valsec = numerical_methods::runge_kutta::rungeKutta(
+        1, 2, 5, 0.1);  // The value of step changed
+    assert(valsec == 3.40600589380261409);
+    std::cout << "Passed Test 2\n";
+
+    std::cout << "Test 3...."
+              << "\n";
+    double valthird = numerical_methods::runge_kutta::rungeKutta(
+        -1, 3, 4, 0.1);  // Tested with negative value
+    assert(valthird == 2.49251005860244268);
+    std::cout << "Passed Test 3\n";
+}
+
+/**
+ * @brief Main function
+ * @returns 0 on exit
+ */
+int main() {
+    test();  // Execute the tests
+    return 0;
+}
@@ -0,0 +1,40 @@
+/**
+ * \file
+ * \brief Method of successive approximations using [fixed-point
+ * iteration](https://en.wikipedia.org/wiki/Fixed-point_iteration) method
+ */
+#include <cmath>
+#include <iostream>
+
+/** equation 1
+ * \f[f(y) = 3y - \cos y -2\f]
+ */
+static float eq(float y) { return (3 * y) - cos(y) - 2; }
+
+/** equation 2
+ * \f[f(y) = \frac{\cos y+2}{2}\f]
+ */
+static float eqd(float y) { return 0.5 * (cos(y) + 2); }
+
+/** Main function */
+int main() {
+    float y, x1, x2, sum;
+    int i, n;
+
+    for (i = 0; i < 10; i++) {
+        sum = eq(y);
+        std::cout << "value of equation at " << i << " " << sum << "\n";
+        y++;
+    }
+    std::cout << "enter the x1->";
+    std::cin >> x1;
+    std::cout << "enter the no iteration to perform->\n";
+    std::cin >> n;
+
+    for (i = 0; i <= n; i++) {
+        x2 = eqd(x1);
+        std::cout << "\nenter the x2->" << x2;
+        x1 = x2;
+    }
+    return 0;
+}