HeatEq1DNeumannBoundary.hh

Go to the documentation of this file.

#ifndef NPSTAT_HEATEQ1DNEUMANNBOUNDARY_HH_
#define NPSTAT_HEATEQ1DNEUMANNBOUNDARY_HH_
 
/*!
// \file HeatEq1DNeumannBoundary.hh
//
// \brief Solution of 1-d heat equation with Neumann boundary conditions
//        (no heat transfer through the boundaries). Useful mainly for
//        generating doubly stochastic matrices.
//
// Author: I. Volobouev
//
// January 2020
*/
 
#include <vector>
 
#include "npstat/nm/Matrix.hh"
 
namespace npstat {
    class HeatEq1DNeumannBoundary
    {
    public:
        /**
        // The number of spatial discretization points, nx, must be
        // at least 3
        */
        HeatEq1DNeumannBoundary(const double* thermalDiffusivity, unsigned nx,
                                double spatialStep, double dt);
 
        //@{
        /** A simple inspector of object properties */
        inline double dt() const {return dt_;}
        inline double time() const {return dt_*cnt_;}
        inline unsigned nCoords() const {return nx_;}
        inline double intervalWidth() const {return intervalWidth_;}
        //@}
 
        /** Evolve the system by time dt */
        void step();
 
        /** Number of time steps made so far */
        inline unsigned long nSteps() const {return cnt_;} 
 
        /**
        // GF standard deviation for the given coordinate index.
        // The index should be smaller than the "nx" argument
        // provided in the constructor.
        */
        double gfWidth(unsigned index) const;
 
        /**
        // Note that the matrix representing the Green's Function
        // changes after each step. To get the solution of the
        // heat equation, multiply the GF matrix by the column of
        // intitial conditions.
        */
        Matrix<double> GF() const;
 
        /** Transposed Green's function */
        Matrix<double> GFT() const;
 
        /**
        // Doubly stochastic approximation to the Green's function.
        // Arguments "tol" and "maxIter" will be passed to the
        // "makeCopulaSteps" method of the ArrayND class used to
        // build the approximation. If "maxIter" is 0, a less reliable
        // non-iterative method will be used that can produce some
        // negative entries in the matrix.
        */
        Matrix<double> doublyStochasticGF(double tol, unsigned maxIter) const;
 
        /**
        // Doubly stochastic approximation to the transposed Green's
        // function. Arguments "tol" and "maxIter" will be passed to
        // the "makeCopulaSteps" method of the ArrayND class used to
        // build the approximation. If "maxIter" is 0, a less reliable
        // non-iterative method will be used that can produce some
        // negative entries in the matrix.
        */
        Matrix<double> doublyStochasticGFT(double tol, unsigned maxIter) const;
 
        /**
        // Matrix representing the differencing scheme. Note that
        // this code is rather slow as it requires O(nx^3) operations.
        */
        Matrix<lapack_double> diffScheme() const;
 
        /**
        // Maximum absolute eigenvalue of the differencing scheme
        // matrix (for use in stability analysis). Note that this
        // code is rather slow as it requires O(nx^3) operations.
        */
        double maxAbsEigen() const;
 
    private:
        HeatEq1DNeumannBoundary();
 
        std::vector<lapack_double> gammas_;
        std::vector<lapack_double> D_;
        std::vector<lapack_double> DL_;
        std::vector<lapack_double> DU_;
        std::vector<lapack_double> DU2_;
        Matrix<lapack_double> m0_;
        Matrix<lapack_double> m1_;
        std::vector<int> IPIV_;
        double dt_;
        double intervalWidth_;
        unsigned long cnt_;
        unsigned nx_;
    };
}
 
#endif // NPSTAT_HEATEQ1DNEUMANNBOUNDARY_HH_