Index: trunk/npstat/nm/HeatEq1DNeumannBoundary.hh =================================================================== --- trunk/npstat/nm/HeatEq1DNeumannBoundary.hh (revision 692) +++ trunk/npstat/nm/HeatEq1DNeumannBoundary.hh (revision 693) @@ -1,94 +1,94 @@ #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 boundary). Useful mainly for +// (no heat transfer through the boundaries). Useful mainly for // generating doubly stochastic matrices. // // Author: I. Volobouev // // January 2020 */ #include #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_;} /** // 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 GF() 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 doublyStochasticGF(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 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 gammas_; std::vector D_; std::vector DL_; std::vector DU_; std::vector DU2_; Matrix m0_; Matrix m1_; std::vector IPIV_; double dt_; double intervalWidth_; unsigned long cnt_; unsigned nx_; }; } #endif // NPSTAT_HEATEQ1DNEUMANNBOUNDARY_HH_ Index: trunk/examples/Python/heat_equation_1d.py =================================================================== --- trunk/examples/Python/heat_equation_1d.py (revision 0) +++ trunk/examples/Python/heat_equation_1d.py (revision 693) @@ -0,0 +1,73 @@ +#!/usr/bin/env python3 +""" +This example illustrates the NPStat code that can be used to solve +one-dimensional heat equation with Neumann boundary conditions (no heat +transfer through the boundaries) and variable thermal diffusivity. The +Green's function of such an equation is a doubly stochastic matrix, as +it transforms a uniform initial temperature distribution into uniform +and also conserves energy. Therefore, this equation can be used to +generate variable bandwidth doubly stochastic filters useful for copula +smoothing. +""" + +__author__="Igor Volobouev (i.volobouev@ttu.edu)" +__version__="1.0" +__date__ ="Jan 31 2020" + +# Perform necessary imports +import npstat as ns +import numpy as np +import matplotlib.pyplot as plt +from scipy import interpolate +from npstat_utils import BoxND +from npstat_plotting import colormeshNumpyArray + +# Some parameters for the script +nDiscrete = 51 # Number of spatial discretization points +spatialStep = 1.0 # Spatial step size +boundaryDiffusivity = 5.0 # Thermal diffusivity at the boundaries of + # the interval +midpointDiffusivity = 1.0 # Thermal diffusivity at the center of the + # interval. The diffusivity will be modelled + # with a second order polynomial. +maxTimeSteps = 200 # Maximum number of time steps to make + +# Build the thermal diffusivity function +nIntervals = nDiscrete - 1 +xmax = nIntervals*spatialStep +thD = (boundaryDiffusivity, midpointDiffusivity, boundaryDiffusivity) +thDFunction = interpolate.interp1d((0.0, xmax/2.0, xmax), thD, 'quadratic') + +# Collection of thermal diffusivity values along the interval +thDValues = thDFunction(np.linspace(0.0, xmax, nDiscrete)) + +# Even though the heat equation is solved using the unconditionally +# stable Crank-Nicolson scheme, a high fidelity Green's function needs +# a small time step, like in the FTCS scheme. +maxDiffusivity = max(boundaryDiffusivity, midpointDiffusivity) +dt = spatialStep*spatialStep/2.0/maxDiffusivity + +# The heat equation solver +eq = ns.HeatEq1DNeumannBoundary(thDValues, spatialStep, dt) + +# Prepare to plot the Green's function +box = BoxND((0.0, xmax), (0.0, xmax)) +fig = plt.figure(figsize=(6, 6)) +ax = fig.add_subplot(111) +plt.draw() + +# Turn interactive mode on for dynamic updates +plt.ion() + +# Run the solver updating the plot +for counter in range(maxTimeSteps): + eq.step() + gf = ns.matrixToNumpy(eq.GF()) + colormeshNumpyArray(ax, gf, box) + ax.set_title("Heat Equation Green's Function," + " t = {:.1f}".format(eq.time())) + plt.pause(0.01) + +# Turn off the interactive mode +plt.ioff() +plt.show() Property changes on: trunk/examples/Python/heat_equation_1d.py ___________________________________________________________________ Added: svn:executable ## -0,0 +1 ## +* \ No newline at end of property Index: trunk/examples/Python/00README.txt =================================================================== --- trunk/examples/Python/00README.txt (revision 692) +++ trunk/examples/Python/00README.txt (revision 693) @@ -1,91 +1,95 @@ This directory contains the following examples: brute_force_kde_1d.py -- Brute-force kernel density estimation (KDE). The most accurate KDE method for small samples. brute_force_osde_1d.py -- Brute-force orthogonal series density estimation (OSDE). Works for small samples. comparison_density_estimation.py -- This script illustrates how to estimate 1-d densities via a comparison distribution. copula_2d_bernstein.py -- This script illustrates how to estimate a multivariate density by splitting it into a product of copula and marginals. The copula is estimated by smoothing it with Bernstein polynomials. density_interpolation_1d_1p.py -- Nonparametric interpolation of univariate densities as a function of a single predictor. density_interpolation_1d_Np.py -- Nonparametric interpolation of univariate densities as a function of multiple predictors. density_interpolation_Nd_Np.py -- Nonparametric interpolation of multivariate densities as a function of multiple predictors. discrete_osde_1d.py -- Orthogonal series density estimation (OSDE) for discretized samples. ems_unfolding_1d.py -- EMS (expectation-maximization with smoothing) unfolding of 1-d data. +heat_equation_1d.py -- Solve a 1-d heat equation with variable + thermal diffusivity and Neumann boundary + conditions. + johnson_curves.py -- Plot some curves from the Johnson system of distributions. kde_1d_by_fft.py -- Efficient kernel density estimation in 1-d with FFT. kde_2d_by_fft.py -- Efficient kernel density estimation in 2-d with FFT. Arbitrary-dimensional KDE would be essentially identical to this example. local_regression_1d.py -- Univariate local least squares polynomial regression. local_regression_2d.py -- Local least squares polynomial regression with two predictors. lorpe_1d_cv.py -- 1-d density estimation by local orthogonal polynomial expansion (LOrPE) with cross-validation. lorpe_1d.py -- 1-d density estimation by LOrPE with a simple convenience driver. lorpe_2d_cv.py -- Multivariate density estimation by local orthogonal polynomial expansion (LOrPE) with cross-validation. lorpe_2d.py -- Multivariate density estimation by LOrPE with a simple plug-in bandwidth guess. persistence.py -- Illustration of the NPStat object persistence via Geners archives. plot_density_1d_matplotlib.py -- Plot an AbsDistribution1D density using Matplotlib. plot_histo_1d_bars.py -- Plot a 1-d HistoND as barplot. plot_histo_1d_horizontal.py -- Plot a 1-d HistoND with bins oriented horizontally, with Matplotlib. plot_histo_1d_matplotlib.py -- Plot a 1-d HistoND with Matplotlib. plot_histo_2d_bars.py -- Plot a 2-d HistoND with Matplotlib, using bars. plot_histo_2d_colormap.py -- Plot a 2-d HistoND with Matplotlib, using color cells. tp_lorpe_2d_cv.py -- Multivariate LOrPE smoothing using a tensor product of 1-d LOrPE filters variable_bandwidth_kde_1d.py -- Kernel density estimation with variable bandwidth.