Numerical Math for Engineering Students

Numerical Math for Engineering Students
Explore algorithms for solving linear and non-linear systems based on "Numerical Methods for Engineers" by S.K. Gupta.
View Code Samples
Linear & Non-linear Systems
Linear Systems
Form Ax = b, where A is a coefficient matrix, x is unknown, and b is constant.
Non-linear Systems
Variables not related by simple linear equations. Common in physics, engineering, and biology.
Solution Methods
Require numerical approximation. Can exhibit complex behavior like multiple solutions.
Newton-Raphson Method
Initial Guess
Choose starting point x0
Compute Jacobian
Calculate Jacobian matrix J at x0
Solve Linear System
Solve J(x1 - x0) = f(x0) for x1
Iterate
Set x0 = x1 and repeat until convergence
The Solver Implementation
Here's an OOP implementation of the Newton-Raphson method using modern C++ syntax.
#include  #include  #include  #include  class NewtonRaphsonSolver { public: // Constructor taking function, its derivative, and optional parameters NewtonRaphsonSolver(std::function func, std::function deriv, float tolerance = 1e-6, int max_iterations = 100) : m_func(func), m_deriv(deriv), m_tolerance(tolerance), m_max_iterations(max_iterations) {} // Solve method with initial guess float solve(float initial_guess) { float x_old = initial_guess; int iterations = 0; while (iterations < m_max_iterations) { float N = m_func(x_old); float D = m_deriv(x_old); // Check for zero function value (exact root found) if (std::abs(N) < m_tolerance) { return x_old; } // Check for zero derivative (method fails) if (std::abs(D) < m_tolerance) { throw std::runtime_error("Newton-Raphson method failed: derivative is zero"); } // Newton-Raphson update float x_new = x_old - (N / D); // Check for convergence if (std::abs(x_new - x_old) < m_tolerance) { return x_new; } // Update for next iteration x_old = x_new; iterations++; } throw std::runtime_error("Newton-Raphson method failed to converge within maximum iterations"); } // Setters for tolerance and max iterations void setTolerance(double tolerance) { m_tolerance = tolerance; } void setMaxIterations(int max_iter) { m_max_iterations = max_iter; } private: std::function m_func; // Function to find root of std::function m_deriv; // Derivative of function double m_tolerance; // Convergence tolerance int m_max_iterations; // Maximum allowed iterations }; // Example usage: int main() { // Define the function and its derivative (example: find root of x^2 - 2) auto func = [](float x) { return x*x - 2.0f; }; auto deriv = [](float x) { return 2.0f*x; }; // Create solver instance NewtonRaphsonSolver solver(func, deriv); try { // Solve with initial guess float root = solver.solve(1.0f); std::cout << "Approximate root: " << root << std::endl; std::cout << "Function value at root: " << func(root) << std::endl; } catch (const std::exception& e) { std::cerr << "Error: " << e.what() << std::endl; } return 0; }
#include <iostream>
#include <cmath>
#include <functional>
#include <stdexcept>

class NewtonRaphsonSolver {
public:
    // Constructor taking function, its derivative, and optional parameters
    NewtonRaphsonSolver(std::function<float(float)> func, 
                       std::function<float(float)> deriv,
                       float tolerance = 1e-6,
                       int max_iterations = 100)
        : m_func(func), m_deriv(deriv), 
          m_tolerance(tolerance), m_max_iterations(max_iterations) {}

    // Solve method with initial guess
    float solve(float initial_guess) {
        float x_old = initial_guess;
        int iterations = 0;
        
        while (iterations < m_max_iterations) {
            float N = m_func(x_old);
            float D = m_deriv(x_old);
            
            // Check for zero function value (exact root found)
            if (std::abs(N) < m_tolerance) {
                return x_old;
            }
            
            // Check for zero derivative (method fails)
            if (std::abs(D) < m_tolerance) {
                throw std::runtime_error("Newton-Raphson method failed: derivative is zero");
            }
            
            // Newton-Raphson update
            float x_new = x_old - (N / D);
            
            // Check for convergence
            if (std::abs(x_new - x_old) < m_tolerance) {
                return x_new;
            }
            
            // Update for next iteration
            x_old = x_new;
            iterations++;
        }
        
        throw std::runtime_error("Newton-Raphson method failed to converge within maximum iterations");
    }

    // Setters for tolerance and max iterations
    void setTolerance(double tolerance) { m_tolerance = tolerance; }
    void setMaxIterations(int max_iter) { m_max_iterations = max_iter; }

private:
    std::function<float(float)> m_func;    // Function to find root of
    std::function<float(float)> m_deriv;   // Derivative of function
    double m_tolerance;                    // Convergence tolerance
    int m_max_iterations;                  // Maximum allowed iterations
};

// Example usage:
int main() {
    // Define the function and its derivative (example: find root of x^2 - 2)
    auto func = [](float x) { return x*x - 2.0f; };
    auto deriv = [](float x) { return 2.0f*x; };
    
    // Create solver instance
    NewtonRaphsonSolver solver(func, deriv);
    
    try {
        // Solve with initial guess
        float root = solver.solve(1.0f);
        std::cout << "Approximate root: " << root << std::endl;
        std::cout << "Function value at root: " << func(root) << std::endl;
    } catch (const std::exception& e) {
        std::cerr << "Error: " << e.what() << std::endl;
    }
    
    return 0;
}
Tridiagonal System Solvers
ODM Algorithm
Orthogonal decomposition method for tridiagonal systems.
Decomposes into independent subproblems
Solves subproblems in parallel
Combines solutions for final result
O(n) complexity, faster than other methods
Thomas Algorithm
Direct method for tridiagonal systems.
Transforms system into solvable form
Uses forward and backward substitution
Reduces complexity from O(n³) to O(n)
Fast but may be expensive for large systems
Interpolation & Differential Equations
Lagrange Interpolation
Constructs polynomial passing through n+1 data points. Simple but prone to oscillation with high-degree polynomials.
Runge-Kutta 4(2)
Fourth-order method for solving ODEs. Uses four evaluations per step for higher accuracy than lower-order methods.
Least Squares
Finds best-fit line by minimizing sum of squared differences between observed and predicted values.
Boundary Value Problems
Guess shooting parameter
Start with initial value
Solve resulting IVP
Use Runge-Kutta method
Compare with boundary condition
Check discrepancy
Adjust parameter and repeat
Until convergence achieved
Three-Point Method
Select Three Points
Choose x1, x0, x2 where x1 < x0 < x2
Evaluate Function
Calculate y1 = f(x1), y0 = f(x0), y2 = f(x2)
Construct Polynomial
Fit polynomial p(x) through the three points
Calculate Derivative
Take derivative of p(x) at x0
Matrix Design & Implementation
Legacy Version
A classic matrix implementation with proven reliability and stability.
New AI-Developed Version
An innovative matrix implementation created with advanced AI technology, tested with MathLAB.
Matrix Routines
Explore the matrix operations and functionalities provided by the AI-composed matrix routines.
Feel free to explore the detailed C++ code snippets for matrix operations below:
#ifndef __MATRIX_H # define __MATRIX_H # include  # include  # include  # include  # include  # include  using namespace std; template  class CMatrix { public: CMatrix(size_t rows, size_t cols) : rows_(rows), cols_(cols), data_(rows * cols) { std::fill(data_.begin(), data_.end(), (T)0); } CMatrix(size_t rows, size_t cols, const std::vector& data) : rows_(rows), cols_(cols), data_(data) { if (rows * cols != data.size()) { throw std::invalid_argument("invalid data size"); } } CMatrix(const CMatrix& other) : rows_(other.rows_), cols_(other.cols_), data_(other.data_) {} T& operator()(size_t row, size_t col) { return data_[row * cols_ + col]; } const T& operator()(size_t row, size_t col) const { return data_[row * cols_ + col]; } // Matrix operations and functions... friend std::ostream& operator<<(std::ostream& out, const CMatrix& m) { out.precision(std::numeric_limits< double >::max_digits10); for (size_t i = 0; i < m.rows_; i++) { for (size_t j = 0; j < m.cols_; j++) { out << m(i, j) << "\t"; } out << std::endl; } return out; } size_t rows() const { return rows_; } size_t cols() const { return cols_; } protected: // Helper functions for matrix manipulation... private: size_t rows_; size_t cols_; std::vector data_; }; #endif // __MATRIX_H
#ifndef __MATRIX_H
# define __MATRIX_H

# include <vector>
# include <algorithm>
# include <cstddef>
# include <stdexcept>
# include <iomanip>
# include <limits>

using namespace std;

template <typename T>
class CMatrix
{
public:
    CMatrix(size_t rows, size_t cols)
      : rows_(rows), cols_(cols), data_(rows * cols)
    {
        std::fill(data_.begin(), data_.end(), (T)0);
    }

    CMatrix(size_t rows, size_t cols, const std::vector<T>& data)
      : rows_(rows), cols_(cols), data_(data)
    {
        if (rows * cols != data.size())
        {
            throw std::invalid_argument("invalid data size");
        }
    }

    CMatrix(const CMatrix<T>& other)
      : rows_(other.rows_), cols_(other.cols_), data_(other.data_) {}

    T& operator()(size_t row, size_t col)
    {
        return data_[row * cols_ + col];
    }

    const T& operator()(size_t row, size_t col) const
    {
        return data_[row * cols_ + col];
    }

    // Matrix operations and functions...

    friend std::ostream& operator<<(std::ostream& out, const CMatrix<T>& m)
    {
        out.precision(std::numeric_limits< double >::max_digits10);
        for (size_t i = 0; i < m.rows_; i++)
        {
            for (size_t j = 0; j < m.cols_; j++)
            {
                out << m(i, j) << "\t";
            }
            out << std::endl;
        }
        return out;
    }

    size_t rows() const { return rows_; }
    size_t cols() const { return cols_; }

protected:
    // Helper functions for matrix manipulation...

private:
    size_t rows_;
    size_t cols_;
    std::vector<T> data_;
};

#endif // __MATRIX_H
Lagrangian Interpolation
Here's an object-oriented implementation of Lagrangian interpolation using modern C++ features.
#include  #include  #include  #include  class LagrangianInterpolator { public: // Constructor taking data points LagrangianInterpolator(const std::vector>& data_points) { if (data_points.empty()) { throw std::invalid_argument("Data points cannot be empty"); } m_data_points = data_points; } // Add a single data point void addDataPoint(float x, float y) { m_data_points.emplace_back(x, y); } // Add multiple data points void addDataPoints(const std::vector>& points) { m_data_points.insert(m_data_points.end(), points.begin(), points.end()); } // Clear all data points void clearDataPoints() { m_data_points.clear(); } // Perform interpolation at point x float interpolate(float x) const { if (m_data_points.empty()) { throw std::runtime_error("No data points available for interpolation"); } float result = 0.0f; const size_t n = m_data_points.size(); for (size_t j = 0; j < n; ++j) { float term = m_data_points[j].second; // y_j for (size_t i = 0; i < n; ++i) { if (i == j) continue; term *= (x - m_data_points[i].first) / (m_data_points[j].first - m_data_points[i].first); } result += term; } return result; } // Get the number of data points size_t size() const { return m_data_points.size(); } // Print all data points (for debugging/visualization) void printDataPoints() const { std::cout << "Data Points:\n"; std::cout << std::setw(10) << "x" << std::setw(10) << "y" << "\n"; std::cout << "---------------------" << "\n"; for (const auto& point : m_data_points) { std::cout << std::setw(10) << point.first << std::setw(10) << point.second << "\n"; } } private: std::vector> m_data_points; // (x, y) pairs }; // Example usage int main() { try { // Create interpolator with initial data points LagrangianInterpolator interpolator({ {1.0f, 1.0f}, {2.0f, 4.0f}, {3.0f, 9.0f}, // Example: f(x) = x^2 {4.0f, 16.0f} }); // Add another point interpolator.addDataPoint(5.0f, 25.0f); // Print all data points interpolator.printDataPoints(); // Perform interpolation at x = 2.5 float x = 2.5f; float y = interpolator.interpolate(x); std::cout << "\nInterpolated value at x = " << x << " is y = " << y << std::endl; } catch (const std::exception& e) { std::cerr << "Error: " << e.what() << std::endl; return 1; } return 0; }
#include <iostream>
#include <vector>
#include <stdexcept>
#include <iomanip>

class LagrangianInterpolator {
public:
    // Constructor taking data points
    LagrangianInterpolator(const std::vector<std::pair<float, float>>& data_points) {
        if (data_points.empty()) {
            throw std::invalid_argument("Data points cannot be empty");
        }
        m_data_points = data_points;
    }

    // Add a single data point
    void addDataPoint(float x, float y) {
        m_data_points.emplace_back(x, y);
    }

    // Add multiple data points
    void addDataPoints(const std::vector<std::pair<float, float>>& points) {
        m_data_points.insert(m_data_points.end(), points.begin(), points.end());
    }

    // Clear all data points
    void clearDataPoints() {
        m_data_points.clear();
    }

    // Perform interpolation at point x
    float interpolate(float x) const {
        if (m_data_points.empty()) {
            throw std::runtime_error("No data points available for interpolation");
        }

        float result = 0.0f;
        const size_t n = m_data_points.size();

        for (size_t j = 0; j < n; ++j) {
            float term = m_data_points[j].second; // y_j
            for (size_t i = 0; i < n; ++i) {
                if (i == j) continue;
                term *= (x - m_data_points[i].first) / 
                       (m_data_points[j].first - m_data_points[i].first);
            }
            result += term;
        }

        return result;
    }

    // Get the number of data points
    size_t size() const {
        return m_data_points.size();
    }

    // Print all data points (for debugging/visualization)
    void printDataPoints() const {
        std::cout << "Data Points:\n";
        std::cout << std::setw(10) << "x" << std::setw(10) << "y" << "\n";
        std::cout << "---------------------" << "\n";
        for (const auto& point : m_data_points) {
            std::cout << std::setw(10) << point.first 
                      << std::setw(10) << point.second << "\n";
        }
    }

private:
    std::vector<std::pair<float, float>> m_data_points; // (x, y) pairs
};

// Example usage
int main() {
    try {
        // Create interpolator with initial data points
        LagrangianInterpolator interpolator({
            {1.0f, 1.0f},
            {2.0f, 4.0f},
            {3.0f, 9.0f},  // Example: f(x) = x^2
            {4.0f, 16.0f}
        });

        // Add another point
        interpolator.addDataPoint(5.0f, 25.0f);

        // Print all data points
        interpolator.printDataPoints();

        // Perform interpolation at x = 2.5
        float x = 2.5f;
        float y = interpolator.interpolate(x);
        
        std::cout << "\nInterpolated value at x = " << x 
                  << " is y = " << y << std::endl;

    } catch (const std::exception& e) {
        std::cerr << "Error: " << e.what() << std::endl;
        return 1;
    }

    return 0;
}