doubleexpintegrator.hh 10.1 KB
 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- // vi: set et ts=4 sw=4 sts=4: /***************************************************************************** * See the file COPYING for full copying permissions. * * * * This program is free software: you can redistribute it and/or modify * * it under the terms of the GNU General Public License as published by * * the Free Software Foundation, either version 3 of the License, or * * (at your option) any later version. * * * * This program is distributed in the hope that it will be useful, * * but WITHOUT ANY WARRANTY; without even the implied warranty of * * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * * GNU General Public License for more details. * * * * You should have received a copy of the GNU General Public License * * along with this program. If not, see . * * * * This version is modified after the original version by John D. Cook * * see https://www.codeproject.com/ * * Articles/31550/Fast-Numerical-Integration * * which is licensed under BSD-2-clause, which reads as follows: * * Copyright John D. Cook * * * * Redistribution and use in source and binary forms, with or without * * modification, are permitted provided that the following * * conditions are met: * * 1. Redistributions of source code must retain the above * * copyright notice, this list of conditions * and the following disclaimer. * * 2. Redistributions in binary form must reproduce the above * * copyright notice, this list of conditions * * and the following disclaimer * * in the documentation and/or other materials * * provided with the distribution. * * * * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS * * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT * * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS * * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE * * COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, * * INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL * * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE * * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS * * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, * * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE * * OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, * * EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * * * *****************************************************************************/ #ifndef DUMUX_COMMON_DOUBLEEXP_INTEGRATOR_HH #define DUMUX_COMMON_DOUBLEEXP_INTEGRATOR_HH #include #include #include #include namespace Dumux { /*! * \brief Numerical integration in one dimension using the double exponential method of M. Mori. */ template class DoubleExponentialIntegrator { public: /*! * \brief Integrate an analytic function over a finite interval * \param f the integrand (invocable with a single scalar) * \param a lower limit of integration * \param b upper limit of integration * \param targetAbsoluteError desired bound on error * \param numFunctionEvaluations number of function evaluations used  Melanie Lipp committed Mar 19, 2020 77  * \param errorEstimate estimate for error in integration  78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206  * \return The value of the integral */ template::value>...> static Scalar integrate(const Function& f, const Scalar a, const Scalar b, const Scalar targetAbsoluteError, int& numFunctionEvaluations, Scalar& errorEstimate) { // Apply the linear change of variables x = ct + d // $$\int_a^b f(x) dx = c \int_{-1}^1 f( ct + d ) dt$$ // c = (b-a)/2, d = (a+b)/2 const Scalar c = 0.5*(b - a); const Scalar d = 0.5*(a + b); return integrateCore_(f, c, d, targetAbsoluteError, numFunctionEvaluations, errorEstimate, doubleExponentialIntegrationAbcissas, doubleExponentialIntegrationWeights); } /*! * \brief Integrate an analytic function over a finite interval. * \note This version overloaded to not require arguments passed in for * function evaluation counts or error estimates. * \param f the integrand (invocable with a single scalar) * \param a lower integral bound * \param b upper integral bound * \param targetAbsoluteError desired absolute error in the result * \return The value of the integral. */ template::value>...> static Scalar integrate(const Function& f, const Scalar a, const Scalar b, const Scalar targetAbsoluteError) { int numFunctionEvaluations; Scalar errorEstimate; return integrate(f, a, b, targetAbsoluteError, numFunctionEvaluations, errorEstimate); } private: // Integrate f(cx + d) with the given integration constants template static Scalar integrateCore_(const Function& f, const Scalar c, // slope of change of variables const Scalar d, // intercept of change of variables Scalar targetAbsoluteError, int& numFunctionEvaluations, Scalar& errorEstimate, const double* abcissas, const double* weights) { targetAbsoluteError /= c; // Offsets to where each level's integration constants start. // The last element is not a beginning but an end. static const int offsets[] = {1, 4, 7, 13, 25, 49, 97, 193}; static const int numLevels = sizeof(offsets)/sizeof(int) - 1; Scalar newContribution = 0.0; Scalar integral = 0.0; Scalar h = 1.0; errorEstimate = std::numeric_limits::max(); Scalar previousDelta, currentDelta = std::numeric_limits::max(); integral = f(c*abcissas[0] + d) * weights[0]; int i; for (i = offsets[0]; i != offsets[1]; ++i) integral += weights[i]*(f(c*abcissas[i] + d) + f(-c*abcissas[i] + d)); for (int level = 1; level != numLevels; ++level) { h *= 0.5; newContribution = 0.0; for (i = offsets[level]; i != offsets[level+1]; ++i) newContribution += weights[i]*(f(c*abcissas[i] + d) + f(-c*abcissas[i] + d)); newContribution *= h; // difference in consecutive integral estimates previousDelta = currentDelta; using std::abs; currentDelta = abs(0.5*integral - newContribution); integral = 0.5*integral + newContribution; // Once convergence kicks in, error is approximately squared at each step. // Determine whether we're in the convergent region by looking at the trend in the error. if (level == 1) continue; // previousDelta meaningless, so cannot check convergence. // Exact comparison with zero is harmless here. Could possibly be replaced with // a small positive upper limit on the size of currentDelta, but determining // that upper limit would be difficult. At worse, the loop is executed more // times than necessary. But no infinite loop can result since there is // an upper bound on the loop variable. if (currentDelta == 0.0) break; using std::log; const Scalar rate = log( currentDelta )/log( previousDelta ); // previousDelta != 0 or would have been kicked out previously if (rate > 1.9 && rate < 2.1) { // If convergence theory applied perfectly, r would be 2 in the convergence region. // r close to 2 is good enough. We expect the difference between this integral estimate // and the next one to be roughly delta^2. errorEstimate = currentDelta*currentDelta; } else { // Not in the convergence region. Assume only that error is decreasing. errorEstimate = currentDelta; } if (errorEstimate < 0.1*targetAbsoluteError) break; } numFunctionEvaluations = 2*i - 1; errorEstimate *= c; return c*integral; } }; } // end namespace Dumux #endif