doubleexpintegrator.hh 10.1 KB
Newer Older
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 <http://www.gnu.org/licenses/>.   *
 *                                                                           *
 *   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 <cmath>
#include <limits>

#include <dune/common/std/type_traits.hh>
#include <dumux/common/doubleexpintegrationconstants.hh>

namespace Dumux {

/*!
 * \brief Numerical integration in one dimension using the double exponential method of M. Mori.
 */
template<class Scalar>
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's avatar
Melanie Lipp committed
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<class Function,
             typename std::enable_if_t<Dune::Std::is_invocable<Function, Scalar>::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<class Function,
             typename std::enable_if_t<Dune::Std::is_invocable<Function, Scalar>::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<class Function>
    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<Scalar>::max();
        Scalar previousDelta, currentDelta = std::numeric_limits<Scalar>::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