001 /*
002 Copyright (C) 2009 Ueli Hofstetter
003
004 This source code is release under the BSD License.
005
006 This file is part of JQuantLib, a free-software/open-source library
007 for financial quantitative analysts and developers - http://jquantlib.org/
008
009 JQuantLib is free software: you can redistribute it and/or modify it
010 under the terms of the JQuantLib license. You should have received a
011 copy of the license along with this program; if not, please email
012 <jquant-devel@lists.sourceforge.net>. The license is also available online at
013 <http://www.jquantlib.org/index.php/LICENSE.TXT>.
014
015 This program is distributed in the hope that it will be useful, but WITHOUT
016 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
017 FOR A PARTICULAR PURPOSE. See the license for more details.
018
019 JQuantLib is based on QuantLib. http://quantlib.org/
020 When applicable, the original copyright notice follows this notice.
021 */
022
023 /*
024 Copyright (C) 2007 Francois du Vignaud
025 Copyright (C) 2003 Niels Elken Sonderby
026
027 This file is part of QuantLib, a free-software/open-source library
028 for financial quantitative analysts and developers - http://quantlib.org/
029
030 QuantLib is free software: you can redistribute it and/or modify it
031 under the terms of the QuantLib license. You should have received a
032 copy of the license along with this program; if not, please email
033 <quantlib-dev@lists.sf.net>. The license is also available online at
034 <http://quantlib.org/license.shtml>.
035
036 This program is distributed in the hope that it will be useful, but WITHOUT
037 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
038 FOR A PARTICULAR PURPOSE. See the license for more details.
039 */
040
041 package org.jquantlib.math.integrals;
042
043 import org.jquantlib.QL;
044 import org.jquantlib.math.Constants;
045 import org.jquantlib.math.Ops;
046
047 /**
048 * Integral of a 1-dimensional function using the Gauss-Kronrod methods
049 * <p>
050 * This class provide a non-adaptive integration procedure which
051 uses fixed Gauss-Kronrod abscissae to sample the integrand at
052 a maximum of 87 points. It is provided for fast integration
053 of smooth functions.
054 <p>
055 This function applies the Gauss-Kronrod 10-point, 21-point, 43-point
056 and 87-point integration rules in succession until an estimate of the
057 integral of f over (a, b) is achieved within the desired absolute and
058 relative error limits, epsabs and epsrel. The function returns the
059 final approximation, result, an estimate of the absolute error,
060 abserr and the number of function evaluations used, neval. The
061 Gauss-Kronrod rules are designed in such a way that each rule uses
062 all the results of its predecessors, in order to minimize the total
063 number of function evaluations.
064
065 * @author Ueli Hofstetter
066 */
067 public class GaussKronrodNonAdaptive extends KronrodIntegral {
068
069 private double relativeAccuracy_;
070
071 public double relativeAccuracy() {
072 return relativeAccuracy_;
073 }
074
075 public void setRelativeAccuracy(final double relativeAccuracy) {
076 this.relativeAccuracy_ = relativeAccuracy;
077 }
078
079 public GaussKronrodNonAdaptive(final double absoluteAccuracy, final int maxEvaluations, final double relativeAccuracy) {
080 super(absoluteAccuracy, maxEvaluations);
081 this.relativeAccuracy_ = relativeAccuracy;
082 }
083
084 @Override
085 public double integrate(final Ops.DoubleOp f, final double a, final double b) {
086 final double fv1[] = new double[5];
087 final double fv2[] = new double[5];
088 final double fv3[] = new double[5];
089 final double fv4[] = new double[5];
090 final double savfun[] = new double[21]; /* array of function values which have been computed */
091 double res10, res21, res43, res87; /* 10, 21, 43 and 87 point results */
092 double err;
093 double resAbs; /* approximation to the integral of abs(f) */
094 double resasc; /* approximation to the integral of abs(f-i/(b-a)) */
095 double result;
096 int k;
097
098 QL.require(a < b , "b must be greater than a"); // QA:[RG]::verified // TODO: message
099
100 final double halfLength = 0.5 * (b - a);
101 final double center = 0.5 * (b + a);
102 final double fCenter = f.op(center);
103
104 // Compute the integral using the 10- and 21-point formula.
105
106 res10 = 0;
107 res21 = w21b[5] * fCenter;
108 resAbs = w21b[5] * Math.abs(fCenter);
109
110 for (k = 0; k < 5; k++) {
111 final double abscissa = halfLength * x1[k];
112 final double fval1 = f.op(center + abscissa);
113 final double fval2 = f.op(center - abscissa);
114 final double fval = fval1 + fval2;
115 res10 += w10[k] * fval;
116 res21 += w21a[k] * fval;
117 resAbs += w21a[k] * (Math.abs(fval1) + Math.abs(fval2));
118 savfun[k] = fval;
119 fv1[k] = fval1;
120 fv2[k] = fval2;
121 }
122
123 for (k = 0; k < 5; k++) {
124 final double abscissa = halfLength * x2[k];
125 final double fval1 = f.op(center + abscissa);
126 final double fval2 = f.op(center - abscissa);
127 final double fval = fval1 + fval2;
128 res21 += w21b[k] * fval;
129 resAbs += w21b[k] * (Math.abs(fval1) + Math.abs(fval2));
130 savfun[k + 5] = fval;
131 fv3[k] = fval1;
132 fv4[k] = fval2;
133 }
134
135 result = res21 * halfLength;
136 resAbs *= halfLength;
137 final double mean = 0.5 * res21;
138 resasc = w21b[5] * Math.abs(fCenter - mean);
139
140 for (k = 0; k < 5; k++) {
141 resasc += w21a[k] * (Math.abs(fv1[k] - mean) + Math.abs(fv2[k] - mean))
142 + w21b[k] * (Math.abs(fv3[k] - mean) + Math.abs(fv4[k] - mean));
143 }
144
145 err = rescaleError((res21 - res10) * halfLength, resAbs, resasc);
146 resasc *= halfLength;
147
148 // test for convergence.
149 if (err < absoluteAccuracy() || err < relativeAccuracy() * Math.abs(result)) {
150 setAbsoluteError(err);
151 setNumberOfEvaluations(21);
152 return result;
153 }
154
155 /* compute the integral using the 43-point formula. */
156
157 res43 = w43b[11] * fCenter;
158
159 for (k = 0; k < 10; k++) {
160 res43 += savfun[k] * w43a[k];
161 }
162
163 for (k = 0; k < 11; k++) {
164 final double abscissa = halfLength * x3[k];
165 final double fval = (f.op(center + abscissa) + f.op(center - abscissa));
166 res43 += fval * w43b[k];
167 savfun[k + 10] = fval;
168 }
169
170 // test for convergence.
171
172 result = res43 * halfLength;
173 err = rescaleError((res43 - res21) * halfLength, resAbs, resasc);
174
175 if (err < absoluteAccuracy() || err < relativeAccuracy() * Math.abs(result)) {
176 setAbsoluteError(err);
177 setNumberOfEvaluations(43);
178 return result;
179 }
180
181 /* compute the integral using the 87-point formula. */
182
183 res87 = w87b[22] * fCenter;
184
185 for (k = 0; k < 21; k++) {
186 res87 += savfun[k] * w87a[k];
187 }
188
189 for (k = 0; k < 22; k++) {
190 final double abscissa = halfLength * x4[k];
191 res87 += w87b[k] * (f.op(center + abscissa) + f.op(center - abscissa));
192 }
193
194 // test for convergence.
195 result = res87 * halfLength;
196 err = rescaleError((res87 - res43) * halfLength, resAbs, resasc);
197
198 setAbsoluteError(err);
199 setNumberOfEvaluations(87);
200 return result;
201 }
202
203 static double rescaleError(double err, final double resultAbs, final double resultAsc) {
204 err = Math.abs(err);
205 if (resultAsc != 0 && err != 0) {
206 final double scale = Math.pow((200 * err / resultAsc), 1.5);
207 if (scale < 1) {
208 err = resultAsc * scale;
209 } else {
210 err = resultAsc;
211 }
212 }
213 if (resultAbs > Constants.QL_MIN_POSITIVE_REAL / (50 * Constants.QL_EPSILON)) {
214 final double min_err = 50 * Constants.QL_EPSILON * resultAbs;
215 if (min_err > err) {
216 err = min_err;
217 }
218 }
219 return err;
220 }
221
222 }