diff --git a/src/Numerics/Integration/GaussKronrodRule.cs b/src/Numerics/Integration/GaussKronrodRule.cs
index 9d798275..f3e95d6c 100644
--- a/src/Numerics/Integration/GaussKronrodRule.cs
+++ b/src/Numerics/Integration/GaussKronrodRule.cs
@@ -35,178 +35,63 @@
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
// https://github.com/boostorg/math/blob/develop/include/boost/math/quadrature/gauss_kronrod.hpp
+using MathNet.Numerics.Integration.GaussRule;
using System;
using System.Numerics;
namespace MathNet.Numerics.Integration
{
- public static class GaussKronrodRule
+ public class GaussKronrodRule
{
- const double epsilon = 2.2204460492503131e-016;
+ private readonly GaussPoints gaussKronrodPoint;
///
- /// The number of Gauss-Kronrod points. Pre-computed for 15, 31, 41, 51 and 61 points.
+ /// Getter for the order.
///
- static int Order = 15;
-
- static double integrate_non_adaptive_m1_1(Func f, out double error, out double pL1)
+ public int Order
{
- int gauss_start = 2;
- int kronrod_start = 1;
- int gauss_order = ((int)Order - 1) / 2;
-
- double kronrod_result = 0d;
- double gauss_result = 0d;
- double fp, fm;
-
- var KAbscissa = KronrodAbscissa();
- var KWeights = KronrodWeights();
- var GWeights = GaussWeights();
-
- if ((gauss_order & 1) == 1)
- {
- fp = f(0);
- kronrod_result = fp * KWeights[0];
- gauss_result += fp * GWeights[0];
- }
- else
- {
- fp = f(0);
- kronrod_result = fp * KWeights[0];
- gauss_start = 1;
- kronrod_start = 2;
- }
- double L1 = Math.Abs(kronrod_result);
-
- for (int i = gauss_start; i < KAbscissa.Length; i += 2)
+ get
{
- fp = f(KAbscissa[i]);
- fm = f(-KAbscissa[i]);
- kronrod_result += (fp + fm) * KWeights[i];
- L1 += (Math.Abs(fp) + Math.Abs(fm)) * KWeights[i];
- gauss_result += (fp + fm) * GWeights[i / 2];
+ return gaussKronrodPoint.Order;
}
- for (int i = kronrod_start; i < KAbscissa.Length; i += 2)
- {
- fp = f(KAbscissa[i]);
- fm = f(-KAbscissa[i]);
- kronrod_result += (fp + fm) * KWeights[i];
- L1 += (Math.Abs(fp) + Math.Abs(fm)) * KWeights[i];
- }
- pL1 = L1;
- error = Math.Max(Math.Abs(kronrod_result - gauss_result), Math.Abs(kronrod_result * epsilon * 2d));
- return kronrod_result;
}
- static Complex contour_integrate_non_adaptive_m1_1(Func f, out double error, out double pL1)
+ ///
+ /// Getter that returns a clone of the array containing the Kronrod abscissas.
+ ///
+ public double[] KronrodAbscissas
{
- int gauss_start = 2;
- int kronrod_start = 1;
- int gauss_order = ((int)Order - 1) / 2;
-
- Complex kronrod_result = new Complex();
- Complex gauss_result = new Complex();
- Complex fp, fm;
-
- var KAbscissa = KronrodAbscissa();
- var KWeights = KronrodWeights();
- var GWeights = GaussWeights();
-
- if (gauss_order.IsOdd())
+ get
{
- fp = f(0);
- kronrod_result = fp * KWeights[0];
- gauss_result += fp * GWeights[0];
- }
- else
- {
- fp = f(0);
- kronrod_result = fp * KWeights[0];
- gauss_start = 1;
- kronrod_start = 2;
+ return gaussKronrodPoint.Abscissas.Clone() as double[];
}
- double L1 = Complex.Abs(kronrod_result);
-
- for (int i = gauss_start; i < KAbscissa.Length; i += 2)
- {
- fp = f(KAbscissa[i]);
- fm = f(-KAbscissa[i]);
- kronrod_result += (fp + fm) * KWeights[i];
- L1 += (Complex.Abs(fp) + Complex.Abs(fm)) * KWeights[i];
- gauss_result += (fp + fm) * GWeights[i / 2];
- }
- for (int i = kronrod_start; i < KAbscissa.Length; i += 2)
- {
- fp = f(KAbscissa[i]);
- fm = f(-KAbscissa[i]);
- kronrod_result += (fp + fm) * KWeights[i];
- L1 += (Complex.Abs(fp) + Complex.Abs(fm)) * KWeights[i];
- }
- pL1 = L1;
- error = Math.Max(Complex.Abs(kronrod_result - gauss_result), Complex.Abs(kronrod_result * epsilon * 2d));
- return kronrod_result;
}
- static double recursive_adaptive_integrate(Func f, double a, double b, int max_levels, double rel_tol, double abs_tol, out double error, out double L1)
+ ///
+ /// Getter that returns a clone of the array containing the Kronrod weights.
+ ///
+ public double[] KronrodWeights
{
- double error_local;
- double mean = (b + a) / 2;
- double scale = (b - a) / 2;
-
- var r1 = integrate_non_adaptive_m1_1((x) => f(scale * x + mean), out error_local, out L1);
- var estimate = scale * r1;
-
- var tmp = estimate * rel_tol;
- var abs_tol1 = Math.Abs(tmp);
- if (abs_tol == 0)
- {
- abs_tol = abs_tol1;
- }
-
- if (max_levels > 0 && (abs_tol1 < error_local) && (abs_tol < error_local))
+ get
{
- double mid = (a + b) / 2d;
- double L1_local;
- estimate = recursive_adaptive_integrate(f, a, mid, max_levels - 1, rel_tol, abs_tol / 2, out error, out L1);
- estimate += recursive_adaptive_integrate(f, mid, b, max_levels - 1, rel_tol, abs_tol / 2, out error_local, out L1_local);
- error += error_local;
- L1 += L1_local;
- return estimate;
+ return gaussKronrodPoint.Weights.Clone() as double[];
}
- L1 *= scale;
- error = error_local;
- return estimate;
}
- static Complex contour_recursive_adaptive_integrate(Func f, double a, double b, int max_levels, double rel_tol, double abs_tol, out double error, out double L1)
+ ///
+ /// Getter that returns a clone of the array containing the Gauss weights.
+ ///
+ public double[] GaussWeights
{
- double error_local;
- double mean = (b + a) / 2;
- double scale = (b - a) / 2;
-
- var r1 = contour_integrate_non_adaptive_m1_1((x) => f(scale * x + mean), out error_local, out L1);
- var estimate = scale * r1;
-
- var tmp = estimate * rel_tol;
- var abs_tol1 = Complex.Abs(tmp);
- if (abs_tol == 0)
+ get
{
- abs_tol = abs_tol1;
+ return gaussKronrodPoint.SecondWeights.Clone() as double[];
}
+ }
- if (max_levels > 0 && (abs_tol1 < error_local) && (abs_tol < error_local))
- {
- double mid = (a + b) / 2d;
- double L1_local;
- estimate = contour_recursive_adaptive_integrate(f, a, mid, max_levels - 1, rel_tol, abs_tol / 2, out error, out L1);
- estimate += contour_recursive_adaptive_integrate(f, mid, b, max_levels - 1, rel_tol, abs_tol / 2, out error_local, out L1_local);
- error += error_local;
- L1 += L1_local;
- return estimate;
- }
- L1 *= scale;
- error = error_local;
- return estimate;
+ public GaussKronrodRule(int order)
+ {
+ gaussKronrodPoint = GaussKronrodPointFactory.GetGaussPoint(order);
}
///
@@ -236,7 +121,7 @@ namespace MathNet.Numerics.Integration
return -Integrate(f, intervalEnd, intervalBegin, out error, out L1Norm, targetRelativeError, maximumDepth, order);
}
- Order = order;
+ GaussPoints gaussKronrodPoint = GaussKronrodPointFactory.GetGaussPoint(order);
// (-oo, oo) => [-1, 1]
//
@@ -249,7 +134,7 @@ namespace MathNet.Numerics.Integration
{
return f(t / (1 - t * t)) * (1 + t * t) / ((1 - t * t) * (1 - t * t));
};
- return recursive_adaptive_integrate(u, -1, 1, maximumDepth, targetRelativeError, 0, out error, out L1Norm);
+ return recursive_adaptive_integrate(u, -1, 1, maximumDepth, targetRelativeError, 0, out error, out L1Norm, gaussKronrodPoint);
}
// [a, oo) => [0, 1]
//
@@ -263,7 +148,7 @@ namespace MathNet.Numerics.Integration
{
return 2 * s * f(intervalBegin + (s / (1 - s)) * (s / (1 - s))) / ((1 - s) * (1 - s) * (1 - s));
};
- return recursive_adaptive_integrate(u, 0, 1, maximumDepth, targetRelativeError, 0, out error, out L1Norm);
+ return recursive_adaptive_integrate(u, 0, 1, maximumDepth, targetRelativeError, 0, out error, out L1Norm, gaussKronrodPoint);
}
// (-oo, b] => [-1, 0]
//
@@ -277,7 +162,7 @@ namespace MathNet.Numerics.Integration
{
return -2 * s * f(intervalEnd - s / (1 + s) * (s / (1 + s))) / ((1 + s) * (1 + s) * (1 + s));
};
- return recursive_adaptive_integrate(u, -1, 0, maximumDepth, targetRelativeError, 0, out error, out L1Norm);
+ return recursive_adaptive_integrate(u, -1, 0, maximumDepth, targetRelativeError, 0, out error, out L1Norm, gaussKronrodPoint);
}
// [a, b] => [-1, 1]
//
@@ -290,7 +175,7 @@ namespace MathNet.Numerics.Integration
{
return f((intervalEnd - intervalBegin) / 4 * t * (3 - t * t) + (intervalEnd + intervalBegin) / 2) * 3 * (intervalEnd - intervalBegin) / 4 * (1 - t * t);
};
- return recursive_adaptive_integrate(u, -1, 1, maximumDepth, targetRelativeError, 0d, out error, out L1Norm);
+ return recursive_adaptive_integrate(u, -1, 1, maximumDepth, targetRelativeError, 0d, out error, out L1Norm, gaussKronrodPoint);
}
}
@@ -322,7 +207,7 @@ namespace MathNet.Numerics.Integration
return -ContourIntegrate(f, intervalEnd, intervalBegin, out error, out L1Norm, targetRelativeError, maximumDepth, order);
}
- Order = order;
+ GaussPoints gaussKronrodPoint = GaussKronrodPointFactory.GetGaussPoint(order);
// (-oo, oo) => [-1, 1]
//
@@ -335,7 +220,7 @@ namespace MathNet.Numerics.Integration
{
return f(t / (1 - t * t)) * (1 + t * t) / ((1 - t * t) * (1 - t * t));
};
- return contour_recursive_adaptive_integrate(u, -1, 1, maximumDepth, targetRelativeError, 0, out error, out L1Norm);
+ return contour_recursive_adaptive_integrate(u, -1, 1, maximumDepth, targetRelativeError, 0, out error, out L1Norm, gaussKronrodPoint);
}
// [a, oo) => [0, 1]
//
@@ -349,7 +234,7 @@ namespace MathNet.Numerics.Integration
{
return 2 * s * f(intervalBegin + (s / (1 - s)) * (s / (1 - s))) / ((1 - s) * (1 - s) * (1 - s));
};
- return contour_recursive_adaptive_integrate(u, 0, 1, maximumDepth, targetRelativeError, 0, out error, out L1Norm);
+ return contour_recursive_adaptive_integrate(u, 0, 1, maximumDepth, targetRelativeError, 0, out error, out L1Norm, gaussKronrodPoint);
}
// (-oo, b] => [-1, 0]
//
@@ -363,7 +248,7 @@ namespace MathNet.Numerics.Integration
{
return -2 * s * f(intervalEnd - s / (1 + s) * (s / (1 + s))) / ((1 + s) * (1 + s) * (1 + s));
};
- return contour_recursive_adaptive_integrate(u, -1, 0, maximumDepth, targetRelativeError, 0, out error, out L1Norm);
+ return contour_recursive_adaptive_integrate(u, -1, 0, maximumDepth, targetRelativeError, 0, out error, out L1Norm, gaussKronrodPoint);
}
// [a, b] => [-1, 1]
//
@@ -376,428 +261,168 @@ namespace MathNet.Numerics.Integration
{
return f((intervalEnd - intervalBegin) / 4 * t * (3 - t * t) + (intervalEnd + intervalBegin) / 2) * 3 * (intervalEnd - intervalBegin) / 4 * (1 - t * t);
};
- return contour_recursive_adaptive_integrate(u, -1, 1, maximumDepth, targetRelativeError, 0d, out error, out L1Norm);
+ return contour_recursive_adaptive_integrate(u, -1, 1, maximumDepth, targetRelativeError, 0d, out error, out L1Norm, gaussKronrodPoint);
}
}
- #region Pre-computed Abscissa and weights
-
- static double[] KronrodAbscissa()
+ private static double integrate_non_adaptive_m1_1(Func f, out double error, out double pL1, GaussPoints gaussKronrodPoint)
{
- switch (Order)
+ int gauss_start = 2;
+ int kronrod_start = 1;
+ int gauss_order = (gaussKronrodPoint.Order - 1) / 2;
+
+ double kronrod_result = 0d;
+ double gauss_result = 0d;
+ double fp, fm;
+
+ var KAbscissa = gaussKronrodPoint.Abscissas;
+ var KWeights = gaussKronrodPoint.Weights;
+ var GWeights = gaussKronrodPoint.SecondWeights;
+
+ if ((gauss_order & 1) == 1)
{
- default:
- case 15:
- return PrecomputedKronrodAbscissas[0];
- case 21:
- return PrecomputedKronrodAbscissas[1];
- case 31:
- return PrecomputedKronrodAbscissas[2];
- case 41:
- return PrecomputedKronrodAbscissas[3];
- case 51:
- return PrecomputedKronrodAbscissas[4];
- case 61:
- return PrecomputedKronrodAbscissas[5];
+ fp = f(0);
+ kronrod_result = fp * KWeights[0];
+ gauss_result += fp * GWeights[0];
}
- }
-
- static double[] KronrodWeights()
- {
- switch (Order)
+ else
{
- default:
- case 15:
- return PrecomputedKronrodWeights[0];
- case 21:
- return PrecomputedKronrodWeights[1];
- case 31:
- return PrecomputedKronrodWeights[2];
- case 41:
- return PrecomputedKronrodWeights[3];
- case 51:
- return PrecomputedKronrodWeights[4];
- case 61:
- return PrecomputedKronrodWeights[5];
+ fp = f(0);
+ kronrod_result = fp * KWeights[0];
+ gauss_start = 1;
+ kronrod_start = 2;
}
- }
+ double L1 = Math.Abs(kronrod_result);
- static double[] GaussWeights()
- {
- switch (Order)
+ for (int i = gauss_start; i < KAbscissa.Length; i += 2)
+ {
+ fp = f(KAbscissa[i]);
+ fm = f(-KAbscissa[i]);
+ kronrod_result += (fp + fm) * KWeights[i];
+ L1 += (Math.Abs(fp) + Math.Abs(fm)) * KWeights[i];
+ gauss_result += (fp + fm) * GWeights[i / 2];
+ }
+ for (int i = kronrod_start; i < KAbscissa.Length; i += 2)
{
- default:
- case 15:
- return PrecomputedGaussWeights[0];
- case 21:
- return PrecomputedGaussWeights[1];
- case 31:
- return PrecomputedGaussWeights[2];
- case 41:
- return PrecomputedGaussWeights[3];
- case 51:
- return PrecomputedGaussWeights[4];
- case 61:
- return PrecomputedGaussWeights[5];
+ fp = f(KAbscissa[i]);
+ fm = f(-KAbscissa[i]);
+ kronrod_result += (fp + fm) * KWeights[i];
+ L1 += (Math.Abs(fp) + Math.Abs(fm)) * KWeights[i];
}
+ pL1 = L1;
+ error = Math.Max(Math.Abs(kronrod_result - gauss_result), Math.Abs(kronrod_result * Precision.MachineEpsilon * 2d));
+ return kronrod_result;
}
- ///
- /// precomputed abscissa vector per order 15, 21, 31, 41, 51 and 61
- ///
- static readonly double[][] PrecomputedKronrodAbscissas =
+ private static Complex contour_integrate_non_adaptive_m1_1(Func f, out double error, out double pL1, GaussPoints gaussKronrodPoint)
{
- new[] // 15-point Gauss-Kronrod
- {
- 0.00000000000000000e+00,
- 2.07784955007898468e-01,
- 4.05845151377397167e-01,
- 5.86087235467691130e-01,
- 7.41531185599394440e-01,
- 8.64864423359769073e-01,
- 9.49107912342758525e-01,
- 9.91455371120812639e-01,
- },
- new[] // 21-point Gauss-Kronrod
- {
- 0.00000000000000000e+00,
- 1.48874338981631211e-01,
- 2.94392862701460198e-01,
- 4.33395394129247191e-01,
- 5.62757134668604683e-01,
- 6.79409568299024406e-01,
- 7.80817726586416897e-01,
- 8.65063366688984511e-01,
- 9.30157491355708226e-01,
- 9.73906528517171720e-01,
- 9.95657163025808081e-01,
- },
- new[] // 31-point Gauss-Kronrod
+ int gauss_start = 2;
+ int kronrod_start = 1;
+ int gauss_order = (gaussKronrodPoint.Order - 1) / 2;
+
+ Complex kronrod_result = new Complex();
+ Complex gauss_result = new Complex();
+ Complex fp, fm;
+
+ var KAbscissa = gaussKronrodPoint.Abscissas;
+ var KWeights = gaussKronrodPoint.Weights;
+ var GWeights = gaussKronrodPoint.SecondWeights;
+
+ if (gauss_order.IsOdd())
{
- 0.00000000000000000e+00,
- 1.01142066918717499e-01,
- 2.01194093997434522e-01,
- 2.99180007153168812e-01,
- 3.94151347077563370e-01,
- 4.85081863640239681e-01,
- 5.70972172608538848e-01,
- 6.50996741297416971e-01,
- 7.24417731360170047e-01,
- 7.90418501442465933e-01,
- 8.48206583410427216e-01,
- 8.97264532344081901e-01,
- 9.37273392400705904e-01,
- 9.67739075679139134e-01,
- 9.87992518020485428e-01,
- 9.98002298693397060e-01,
- },
- new[] // 41-point Gauss-Kronrod
+ fp = f(0);
+ kronrod_result = fp * KWeights[0];
+ gauss_result += fp * GWeights[0];
+ }
+ else
{
- 0.00000000000000000e+00,
- 7.65265211334973338e-02,
- 1.52605465240922676e-01,
- 2.27785851141645078e-01,
- 3.01627868114913004e-01,
- 3.73706088715419561e-01,
- 4.43593175238725103e-01,
- 5.10867001950827098e-01,
- 5.75140446819710315e-01,
- 6.36053680726515025e-01,
- 6.93237656334751385e-01,
- 7.46331906460150793e-01,
- 7.95041428837551198e-01,
- 8.39116971822218823e-01,
- 8.78276811252281976e-01,
- 9.12234428251325906e-01,
- 9.40822633831754754e-01,
- 9.63971927277913791e-01,
- 9.81507877450250259e-01,
- 9.93128599185094925e-01,
- 9.98859031588277664e-01,
- },
- new[] // 51-point Gauss-Kronrod
+ fp = f(0);
+ kronrod_result = fp * KWeights[0];
+ gauss_start = 1;
+ kronrod_start = 2;
+ }
+ double L1 = Complex.Abs(kronrod_result);
+
+ for (int i = gauss_start; i < KAbscissa.Length; i += 2)
{
- 0.00000000000000000e+00,
- 6.15444830056850789e-02,
- 1.22864692610710396e-01,
- 1.83718939421048892e-01,
- 2.43866883720988432e-01,
- 3.03089538931107830e-01,
- 3.61172305809387838e-01,
- 4.17885382193037749e-01,
- 4.73002731445714961e-01,
- 5.26325284334719183e-01,
- 5.77662930241222968e-01,
- 6.26810099010317413e-01,
- 6.73566368473468364e-01,
- 7.17766406813084388e-01,
- 7.59259263037357631e-01,
- 7.97873797998500059e-01,
- 8.33442628760834001e-01,
- 8.65847065293275595e-01,
- 8.94991997878275369e-01,
- 9.20747115281701562e-01,
- 9.42974571228974339e-01,
- 9.61614986425842512e-01,
- 9.76663921459517511e-01,
- 9.88035794534077248e-01,
- 9.95556969790498098e-01,
- 9.99262104992609834e-01,
- },
- new[] // 61-point Gauss-Kronrod
+ fp = f(KAbscissa[i]);
+ fm = f(-KAbscissa[i]);
+ kronrod_result += (fp + fm) * KWeights[i];
+ L1 += (Complex.Abs(fp) + Complex.Abs(fm)) * KWeights[i];
+ gauss_result += (fp + fm) * GWeights[i / 2];
+ }
+ for (int i = kronrod_start; i < KAbscissa.Length; i += 2)
{
- 0.00000000000000000e+00,
- 5.14718425553176958e-02,
- 1.02806937966737030e-01,
- 1.53869913608583547e-01,
- 2.04525116682309891e-01,
- 2.54636926167889846e-01,
- 3.04073202273625077e-01,
- 3.52704725530878113e-01,
- 4.00401254830394393e-01,
- 4.47033769538089177e-01,
- 4.92480467861778575e-01,
- 5.36624148142019899e-01,
- 5.79345235826361692e-01,
- 6.20526182989242861e-01,
- 6.60061064126626961e-01,
- 6.97850494793315797e-01,
- 7.33790062453226805e-01,
- 7.67777432104826195e-01,
- 7.99727835821839083e-01,
- 8.29565762382768397e-01,
- 8.57205233546061099e-01,
- 8.82560535792052682e-01,
- 9.05573307699907799e-01,
- 9.26200047429274326e-01,
- 9.44374444748559979e-01,
- 9.60021864968307512e-01,
- 9.73116322501126268e-01,
- 9.83668123279747210e-01,
- 9.91630996870404595e-01,
- 9.96893484074649540e-01,
- 9.99484410050490638e-01,
+ fp = f(KAbscissa[i]);
+ fm = f(-KAbscissa[i]);
+ kronrod_result += (fp + fm) * KWeights[i];
+ L1 += (Complex.Abs(fp) + Complex.Abs(fm)) * KWeights[i];
}
- };
+ pL1 = L1;
+ error = Math.Max(Complex.Abs(kronrod_result - gauss_result), Complex.Abs(kronrod_result * Precision.MachineEpsilon * 2d));
+ return kronrod_result;
+ }
- ///
- /// precomputed weight vector per order 15, 21, 31, 41, 51 and 61
- ///
- static readonly double[][] PrecomputedKronrodWeights =
+ private static double recursive_adaptive_integrate(Func f, double a, double b, int max_levels, double rel_tol, double abs_tol, out double error, out double L1, GaussPoints gaussKronrodPoint)
{
- new[] // 15-point Gauss-Kronrod integration
- {
- 2.09482141084727828e-01,
- 2.04432940075298892e-01,
- 1.90350578064785410e-01,
- 1.69004726639267903e-01,
- 1.40653259715525919e-01,
- 1.04790010322250184e-01,
- 6.30920926299785533e-02,
- 2.29353220105292250e-02,
- },
- new[] // 21-point Gauss-Kronrod integration
- {
- 1.49445554002916906e-01,
- 1.47739104901338491e-01,
- 1.42775938577060081e-01,
- 1.34709217311473326e-01,
- 1.23491976262065851e-01,
- 1.09387158802297642e-01,
- 9.31254545836976055e-02,
- 7.50396748109199528e-02,
- 5.47558965743519960e-02,
- 3.25581623079647275e-02,
- 1.16946388673718743e-02,
- },
- new[] // 31-point Gauss-Kronrod integration
- {
- 1.01330007014791549e-01,
- 1.00769845523875595e-01,
- 9.91735987217919593e-02,
- 9.66427269836236785e-02,
- 9.31265981708253212e-02,
- 8.85644430562117706e-02,
- 8.30805028231330210e-02,
- 7.68496807577203789e-02,
- 6.98541213187282587e-02,
- 6.20095678006706403e-02,
- 5.34815246909280873e-02,
- 4.45897513247648766e-02,
- 3.53463607913758462e-02,
- 2.54608473267153202e-02,
- 1.50079473293161225e-02,
- 5.37747987292334899e-03,
- },
- new[] // 41-point Gauss-Kronrod integration
- {
- 7.66007119179996564e-02,
- 7.63778676720807367e-02,
- 7.57044976845566747e-02,
- 7.45828754004991890e-02,
- 7.30306903327866675e-02,
- 7.10544235534440683e-02,
- 6.86486729285216193e-02,
- 6.58345971336184221e-02,
- 6.26532375547811680e-02,
- 5.91114008806395724e-02,
- 5.51951053482859947e-02,
- 5.09445739237286919e-02,
- 4.64348218674976747e-02,
- 4.16688733279736863e-02,
- 3.66001697582007980e-02,
- 3.12873067770327990e-02,
- 2.58821336049511588e-02,
- 2.03883734612665236e-02,
- 1.46261692569712530e-02,
- 8.60026985564294220e-03,
- 3.07358371852053150e-03,
- },
- new[] // 51-point Gauss-Kronrod integration
+ double error_local;
+ double mean = (b + a) / 2;
+ double scale = (b - a) / 2;
+
+ var r1 = integrate_non_adaptive_m1_1((x) => f(scale * x + mean), out error_local, out L1, gaussKronrodPoint);
+ var estimate = scale * r1;
+
+ var tmp = estimate * rel_tol;
+ var abs_tol1 = Math.Abs(tmp);
+ if (abs_tol == 0)
{
- 6.15808180678329351e-02,
- 6.14711898714253167e-02,
- 6.11285097170530483e-02,
- 6.05394553760458629e-02,
- 5.97203403241740600e-02,
- 5.86896800223942080e-02,
- 5.74371163615678329e-02,
- 5.59508112204123173e-02,
- 5.42511298885454901e-02,
- 5.23628858064074759e-02,
- 5.02776790807156720e-02,
- 4.79825371388367139e-02,
- 4.55029130499217889e-02,
- 4.28728450201700495e-02,
- 4.00838255040323821e-02,
- 3.71162714834155436e-02,
- 3.40021302743293378e-02,
- 3.07923001673874889e-02,
- 2.74753175878517378e-02,
- 2.40099456069532162e-02,
- 2.04353711458828355e-02,
- 1.68478177091282982e-02,
- 1.32362291955716748e-02,
- 9.47397338617415161e-03,
- 5.56193213535671376e-03,
- 1.98738389233031593e-03,
- },
- new[] // 61-point Gauss-Kronrod integration
+ abs_tol = abs_tol1;
+ }
+
+ if (max_levels > 0 && (abs_tol1 < error_local) && (abs_tol < error_local))
{
- 5.14947294294515676e-02,
- 5.14261285374590259e-02,
- 5.12215478492587722e-02,
- 5.08817958987496065e-02,
- 5.04059214027823468e-02,
- 4.97956834270742064e-02,
- 4.90554345550297789e-02,
- 4.81858617570871291e-02,
- 4.71855465692991539e-02,
- 4.60592382710069881e-02,
- 4.48148001331626632e-02,
- 4.34525397013560693e-02,
- 4.19698102151642461e-02,
- 4.03745389515359591e-02,
- 3.86789456247275930e-02,
- 3.68823646518212292e-02,
- 3.49793380280600241e-02,
- 3.29814470574837260e-02,
- 3.09072575623877625e-02,
- 2.87540487650412928e-02,
- 2.65099548823331016e-02,
- 2.41911620780806014e-02,
- 2.18280358216091923e-02,
- 1.94141411939423812e-02,
- 1.69208891890532726e-02,
- 1.43697295070458048e-02,
- 1.18230152534963417e-02,
- 9.27327965951776343e-03,
- 6.63070391593129217e-03,
- 3.89046112709988405e-03,
- 1.38901369867700762e-03,
- },
- };
+ double mid = (a + b) / 2d;
+ double L1_local;
+ estimate = recursive_adaptive_integrate(f, a, mid, max_levels - 1, rel_tol, abs_tol / 2, out error, out L1, gaussKronrodPoint);
+ estimate += recursive_adaptive_integrate(f, mid, b, max_levels - 1, rel_tol, abs_tol / 2, out error_local, out L1_local, gaussKronrodPoint);
+ error += error_local;
+ L1 += L1_local;
+ return estimate;
+ }
+ L1 *= scale;
+ error = error_local;
+ return estimate;
+ }
- ///
- /// precomputed Gauss weight vector per order 7, 10, 15, 20, 25 and 30
- ///
- static readonly double[][] PrecomputedGaussWeights =
+ private static Complex contour_recursive_adaptive_integrate(Func f, double a, double b, int max_levels, double rel_tol, double abs_tol, out double error, out double L1, GaussPoints gaussKronrodPoint)
{
- new [] // 7-point Gauss
- {
- 4.17959183673469388e-01,
- 3.81830050505118945e-01,
- 2.79705391489276668e-01,
- 1.29484966168869693e-01,
- },
- new[] // 10-point Gauss
- {
- 2.95524224714752870e-01,
- 2.69266719309996355e-01,
- 2.19086362515982044e-01,
- 1.49451349150580593e-01,
- 6.66713443086881376e-02,
- },
- new[] // 15-point Gauss
- {
- 2.02578241925561273e-01,
- 1.98431485327111576e-01,
- 1.86161000015562211e-01,
- 1.66269205816993934e-01,
- 1.39570677926154314e-01,
- 1.07159220467171935e-01,
- 7.03660474881081247e-02,
- 3.07532419961172684e-02,
- },
- new[] // 20-point Gauss
- {
- 1.52753387130725851e-01,
- 1.49172986472603747e-01,
- 1.42096109318382051e-01,
- 1.31688638449176627e-01,
- 1.18194531961518417e-01,
- 1.01930119817240435e-01,
- 8.32767415767047487e-02,
- 6.26720483341090636e-02,
- 4.06014298003869413e-02,
- 1.76140071391521183e-02,
- },
- new[] // 25-point Gauss
- {
- 1.23176053726715451e-01,
- 1.22242442990310042e-01,
- 1.19455763535784772e-01,
- 1.14858259145711648e-01,
- 1.08519624474263653e-01,
- 1.00535949067050644e-01,
- 9.10282619829636498e-02,
- 8.01407003350010180e-02,
- 6.80383338123569172e-02,
- 5.49046959758351919e-02,
- 4.09391567013063127e-02,
- 2.63549866150321373e-02,
- 1.13937985010262879e-02,
- },
- new[] // 30-point Gauss
+ double error_local;
+ double mean = (b + a) / 2;
+ double scale = (b - a) / 2;
+
+ var r1 = contour_integrate_non_adaptive_m1_1((x) => f(scale * x + mean), out error_local, out L1, gaussKronrodPoint);
+ var estimate = scale * r1;
+
+ var tmp = estimate * rel_tol;
+ var abs_tol1 = Complex.Abs(tmp);
+ if (abs_tol == 0)
{
- 1.02852652893558840e-01,
- 1.01762389748405505e-01,
- 9.95934205867952671e-02,
- 9.63687371746442596e-02,
- 9.21225222377861287e-02,
- 8.68997872010829798e-02,
- 8.07558952294202154e-02,
- 7.37559747377052063e-02,
- 6.59742298821804951e-02,
- 5.74931562176190665e-02,
- 4.84026728305940529e-02,
- 3.87991925696270496e-02,
- 2.87847078833233693e-02,
- 1.84664683110909591e-02,
- 7.96819249616660562e-03,
+ abs_tol = abs_tol1;
}
- };
- #endregion Pre-computed Abscissa and weights
+ if (max_levels > 0 && (abs_tol1 < error_local) && (abs_tol < error_local))
+ {
+ double mid = (a + b) / 2d;
+ double L1_local;
+ estimate = contour_recursive_adaptive_integrate(f, a, mid, max_levels - 1, rel_tol, abs_tol / 2, out error, out L1, gaussKronrodPoint);
+ estimate += contour_recursive_adaptive_integrate(f, mid, b, max_levels - 1, rel_tol, abs_tol / 2, out error_local, out L1_local, gaussKronrodPoint);
+ error += error_local;
+ L1 += L1_local;
+ return estimate;
+ }
+ L1 *= scale;
+ error = error_local;
+ return estimate;
+ }
}
}
diff --git a/src/Numerics/Integration/GaussRule/GaussKronrodPoint.cs b/src/Numerics/Integration/GaussRule/GaussKronrodPoint.cs
new file mode 100644
index 00000000..8659663a
--- /dev/null
+++ b/src/Numerics/Integration/GaussRule/GaussKronrodPoint.cs
@@ -0,0 +1,382 @@
+using System;
+using System.Collections.Generic;
+
+namespace MathNet.Numerics.Integration.GaussRule
+{
+ ///
+ /// Contains a method to compute the Gauss-Kronrod abscissas/weights and precomputed abscissas/weights for orders 15, 21, 31, 41, 51, 61.
+ ///
+ internal static partial class GaussKronrodPoint
+ {
+ ///
+ /// Precomputed abscissas/weights for orders 15, 21, 31, 41, 51, 61.
+ ///
+ internal static readonly Dictionary PreComputed = new Dictionary
+ {
+ { 15, new GaussPoints(15,
+ new[] // 15-point Gauss-Kronrod Abscissa
+ {
+ 0.00000000000000000e+00,
+ 2.07784955007898468e-01,
+ 4.05845151377397167e-01,
+ 5.86087235467691130e-01,
+ 7.41531185599394440e-01,
+ 8.64864423359769073e-01,
+ 9.49107912342758525e-01,
+ 9.91455371120812639e-01,
+ },
+ new[] // 15-point Gauss-Kronrod Weights
+ {
+ 2.09482141084727828e-01,
+ 2.04432940075298892e-01,
+ 1.90350578064785410e-01,
+ 1.69004726639267903e-01,
+ 1.40653259715525919e-01,
+ 1.04790010322250184e-01,
+ 6.30920926299785533e-02,
+ 2.29353220105292250e-02,
+ }, 7,
+ new[] // 7-point Gauss Weights
+ {
+ 4.17959183673469388e-01,
+ 3.81830050505118945e-01,
+ 2.79705391489276668e-01,
+ 1.29484966168869693e-01,
+ })
+ },
+ { 21, new GaussPoints(21,
+ new[] // 21-point Gauss-Kronrod Abscissa
+ {
+ 0.00000000000000000e+00,
+ 1.48874338981631211e-01,
+ 2.94392862701460198e-01,
+ 4.33395394129247191e-01,
+ 5.62757134668604683e-01,
+ 6.79409568299024406e-01,
+ 7.80817726586416897e-01,
+ 8.65063366688984511e-01,
+ 9.30157491355708226e-01,
+ 9.73906528517171720e-01,
+ 9.95657163025808081e-01,
+ },
+ new[] // 21-point Gauss-Kronrod Weights
+ {
+ 1.49445554002916906e-01,
+ 1.47739104901338491e-01,
+ 1.42775938577060081e-01,
+ 1.34709217311473326e-01,
+ 1.23491976262065851e-01,
+ 1.09387158802297642e-01,
+ 9.31254545836976055e-02,
+ 7.50396748109199528e-02,
+ 5.47558965743519960e-02,
+ 3.25581623079647275e-02,
+ 1.16946388673718743e-02,
+ }, 10,
+ new[] // 10-point Gauss Weights
+ {
+ 2.95524224714752870e-01,
+ 2.69266719309996355e-01,
+ 2.19086362515982044e-01,
+ 1.49451349150580593e-01,
+ 6.66713443086881376e-02,
+ })
+ },
+ { 31, new GaussPoints(31,
+ new[] // 31-point Gauss-Kronrod Abscissa
+ {
+ 0.00000000000000000e+00,
+ 1.01142066918717499e-01,
+ 2.01194093997434522e-01,
+ 2.99180007153168812e-01,
+ 3.94151347077563370e-01,
+ 4.85081863640239681e-01,
+ 5.70972172608538848e-01,
+ 6.50996741297416971e-01,
+ 7.24417731360170047e-01,
+ 7.90418501442465933e-01,
+ 8.48206583410427216e-01,
+ 8.97264532344081901e-01,
+ 9.37273392400705904e-01,
+ 9.67739075679139134e-01,
+ 9.87992518020485428e-01,
+ 9.98002298693397060e-01,
+ },
+ new[] // 31-point Gauss-Kronrod Weights
+ {
+ 1.01330007014791549e-01,
+ 1.00769845523875595e-01,
+ 9.91735987217919593e-02,
+ 9.66427269836236785e-02,
+ 9.31265981708253212e-02,
+ 8.85644430562117706e-02,
+ 8.30805028231330210e-02,
+ 7.68496807577203789e-02,
+ 6.98541213187282587e-02,
+ 6.20095678006706403e-02,
+ 5.34815246909280873e-02,
+ 4.45897513247648766e-02,
+ 3.53463607913758462e-02,
+ 2.54608473267153202e-02,
+ 1.50079473293161225e-02,
+ 5.37747987292334899e-03,
+ }, 15,
+ new[] // 15-point Gauss Weights
+ {
+ 2.02578241925561273e-01,
+ 1.98431485327111576e-01,
+ 1.86161000015562211e-01,
+ 1.66269205816993934e-01,
+ 1.39570677926154314e-01,
+ 1.07159220467171935e-01,
+ 7.03660474881081247e-02,
+ 3.07532419961172684e-02,
+ })
+ },
+ { 41, new GaussPoints(41,
+ new[] // 41-point Gauss-Kronrod Abscissa
+ {
+ 0.00000000000000000e+00,
+ 7.65265211334973338e-02,
+ 1.52605465240922676e-01,
+ 2.27785851141645078e-01,
+ 3.01627868114913004e-01,
+ 3.73706088715419561e-01,
+ 4.43593175238725103e-01,
+ 5.10867001950827098e-01,
+ 5.75140446819710315e-01,
+ 6.36053680726515025e-01,
+ 6.93237656334751385e-01,
+ 7.46331906460150793e-01,
+ 7.95041428837551198e-01,
+ 8.39116971822218823e-01,
+ 8.78276811252281976e-01,
+ 9.12234428251325906e-01,
+ 9.40822633831754754e-01,
+ 9.63971927277913791e-01,
+ 9.81507877450250259e-01,
+ 9.93128599185094925e-01,
+ 9.98859031588277664e-01,
+ },
+ new[] // 41-point Gauss-Kronrod Weights
+ {
+ 7.66007119179996564e-02,
+ 7.63778676720807367e-02,
+ 7.57044976845566747e-02,
+ 7.45828754004991890e-02,
+ 7.30306903327866675e-02,
+ 7.10544235534440683e-02,
+ 6.86486729285216193e-02,
+ 6.58345971336184221e-02,
+ 6.26532375547811680e-02,
+ 5.91114008806395724e-02,
+ 5.51951053482859947e-02,
+ 5.09445739237286919e-02,
+ 4.64348218674976747e-02,
+ 4.16688733279736863e-02,
+ 3.66001697582007980e-02,
+ 3.12873067770327990e-02,
+ 2.58821336049511588e-02,
+ 2.03883734612665236e-02,
+ 1.46261692569712530e-02,
+ 8.60026985564294220e-03,
+ 3.07358371852053150e-03,
+ }, 20,
+ new[] // 20-point Gauss Weights
+ {
+ 1.52753387130725851e-01,
+ 1.49172986472603747e-01,
+ 1.42096109318382051e-01,
+ 1.31688638449176627e-01,
+ 1.18194531961518417e-01,
+ 1.01930119817240435e-01,
+ 8.32767415767047487e-02,
+ 6.26720483341090636e-02,
+ 4.06014298003869413e-02,
+ 1.76140071391521183e-02,
+ })
+ },
+ { 51, new GaussPoints(51,
+ new[] // 51-point Gauss-Kronrod Abscissa
+ {
+ 0.00000000000000000e+00,
+ 6.15444830056850789e-02,
+ 1.22864692610710396e-01,
+ 1.83718939421048892e-01,
+ 2.43866883720988432e-01,
+ 3.03089538931107830e-01,
+ 3.61172305809387838e-01,
+ 4.17885382193037749e-01,
+ 4.73002731445714961e-01,
+ 5.26325284334719183e-01,
+ 5.77662930241222968e-01,
+ 6.26810099010317413e-01,
+ 6.73566368473468364e-01,
+ 7.17766406813084388e-01,
+ 7.59259263037357631e-01,
+ 7.97873797998500059e-01,
+ 8.33442628760834001e-01,
+ 8.65847065293275595e-01,
+ 8.94991997878275369e-01,
+ 9.20747115281701562e-01,
+ 9.42974571228974339e-01,
+ 9.61614986425842512e-01,
+ 9.76663921459517511e-01,
+ 9.88035794534077248e-01,
+ 9.95556969790498098e-01,
+ 9.99262104992609834e-01,
+ },
+ new[] // 51-point Gauss-Kronrod Weights
+ {
+ 6.15808180678329351e-02,
+ 6.14711898714253167e-02,
+ 6.11285097170530483e-02,
+ 6.05394553760458629e-02,
+ 5.97203403241740600e-02,
+ 5.86896800223942080e-02,
+ 5.74371163615678329e-02,
+ 5.59508112204123173e-02,
+ 5.42511298885454901e-02,
+ 5.23628858064074759e-02,
+ 5.02776790807156720e-02,
+ 4.79825371388367139e-02,
+ 4.55029130499217889e-02,
+ 4.28728450201700495e-02,
+ 4.00838255040323821e-02,
+ 3.71162714834155436e-02,
+ 3.40021302743293378e-02,
+ 3.07923001673874889e-02,
+ 2.74753175878517378e-02,
+ 2.40099456069532162e-02,
+ 2.04353711458828355e-02,
+ 1.68478177091282982e-02,
+ 1.32362291955716748e-02,
+ 9.47397338617415161e-03,
+ 5.56193213535671376e-03,
+ 1.98738389233031593e-03,
+ }, 25,
+ new[] // 25-point Gauss Weights
+ {
+ 1.23176053726715451e-01,
+ 1.22242442990310042e-01,
+ 1.19455763535784772e-01,
+ 1.14858259145711648e-01,
+ 1.08519624474263653e-01,
+ 1.00535949067050644e-01,
+ 9.10282619829636498e-02,
+ 8.01407003350010180e-02,
+ 6.80383338123569172e-02,
+ 5.49046959758351919e-02,
+ 4.09391567013063127e-02,
+ 2.63549866150321373e-02,
+ 1.13937985010262879e-02,
+ })
+ },
+ { 61, new GaussPoints(61,
+ new[] // 61-point Gauss-Kronrod Abscissa
+ {
+ 0.00000000000000000e+00,
+ 5.14718425553176958e-02,
+ 1.02806937966737030e-01,
+ 1.53869913608583547e-01,
+ 2.04525116682309891e-01,
+ 2.54636926167889846e-01,
+ 3.04073202273625077e-01,
+ 3.52704725530878113e-01,
+ 4.00401254830394393e-01,
+ 4.47033769538089177e-01,
+ 4.92480467861778575e-01,
+ 5.36624148142019899e-01,
+ 5.79345235826361692e-01,
+ 6.20526182989242861e-01,
+ 6.60061064126626961e-01,
+ 6.97850494793315797e-01,
+ 7.33790062453226805e-01,
+ 7.67777432104826195e-01,
+ 7.99727835821839083e-01,
+ 8.29565762382768397e-01,
+ 8.57205233546061099e-01,
+ 8.82560535792052682e-01,
+ 9.05573307699907799e-01,
+ 9.26200047429274326e-01,
+ 9.44374444748559979e-01,
+ 9.60021864968307512e-01,
+ 9.73116322501126268e-01,
+ 9.83668123279747210e-01,
+ 9.91630996870404595e-01,
+ 9.96893484074649540e-01,
+ 9.99484410050490638e-01,
+ },
+ new[] // 61-point Gauss-Kronrod Weights
+ {
+ 5.14947294294515676e-02,
+ 5.14261285374590259e-02,
+ 5.12215478492587722e-02,
+ 5.08817958987496065e-02,
+ 5.04059214027823468e-02,
+ 4.97956834270742064e-02,
+ 4.90554345550297789e-02,
+ 4.81858617570871291e-02,
+ 4.71855465692991539e-02,
+ 4.60592382710069881e-02,
+ 4.48148001331626632e-02,
+ 4.34525397013560693e-02,
+ 4.19698102151642461e-02,
+ 4.03745389515359591e-02,
+ 3.86789456247275930e-02,
+ 3.68823646518212292e-02,
+ 3.49793380280600241e-02,
+ 3.29814470574837260e-02,
+ 3.09072575623877625e-02,
+ 2.87540487650412928e-02,
+ 2.65099548823331016e-02,
+ 2.41911620780806014e-02,
+ 2.18280358216091923e-02,
+ 1.94141411939423812e-02,
+ 1.69208891890532726e-02,
+ 1.43697295070458048e-02,
+ 1.18230152534963417e-02,
+ 9.27327965951776343e-03,
+ 6.63070391593129217e-03,
+ 3.89046112709988405e-03,
+ 1.38901369867700762e-03,
+ }, 30,
+ new[] // 30-point Gauss Weights
+ {
+ 1.02852652893558840e-01,
+ 1.01762389748405505e-01,
+ 9.95934205867952671e-02,
+ 9.63687371746442596e-02,
+ 9.21225222377861287e-02,
+ 8.68997872010829798e-02,
+ 8.07558952294202154e-02,
+ 7.37559747377052063e-02,
+ 6.59742298821804951e-02,
+ 5.74931562176190665e-02,
+ 4.84026728305940529e-02,
+ 3.87991925696270496e-02,
+ 2.87847078833233693e-02,
+ 1.84664683110909591e-02,
+ 7.96819249616660562e-03,
+ })
+ },
+ };
+ }
+
+ ///
+ /// Contains a method to compute the Gauss-Kronrod abscissas/weights.
+ ///
+ internal static partial class GaussKronrodPoint
+ {
+ ///
+ /// Computes the Gauss-Kronrod abscissas/weights and Gauss weights.
+ ///
+ /// Defines an Nth order Gauss-Kronrod rule. The order also defines the number of abscissas and weights for the rule.
+ /// Required precision to compute the abscissas/weights.
+ /// Object containing the non-negative abscissas/weights, order.
+ internal static GaussPoints Generate(int order, double targetAbsoluteTolerance = 1E-10)
+ {
+ throw new NotSupportedException(string.Format("The order of Gauss Kronrod, {0} is not yet supported", order));
+ }
+ }
+}
diff --git a/src/Numerics/Integration/GaussRule/GaussKronrodPointFactory.cs b/src/Numerics/Integration/GaussRule/GaussKronrodPointFactory.cs
new file mode 100644
index 00000000..ee294d37
--- /dev/null
+++ b/src/Numerics/Integration/GaussRule/GaussKronrodPointFactory.cs
@@ -0,0 +1,36 @@
+using System;
+
+namespace MathNet.Numerics.Integration.GaussRule
+{
+ ///
+ /// Creates a Gauss-Kronrod point.
+ ///
+ internal static class GaussKronrodPointFactory
+ {
+ [ThreadStatic]
+ private static GaussPoints gaussKronrodPoint;
+
+ ///
+ /// Getter for the GaussKronrodPoint.
+ ///
+ /// Defines an Nth order Gauss-Kronrod rule. Precomputed Gauss-Kronrod abscissas/weights for orders 15, 21, 31, 41, 51, 61 are used, otherwise they're calculated on the fly.
+ /// Required precision to compute the abscissas/weights. 1e-10 is usually fine.
+ /// Object containing the non-negative abscissas/weights, and order.
+ public static GaussPoints GetGaussPoint(int order, double targetAbsoluteTolerance = 1E-10)
+ {
+ // Try to get the GaussKronrodPoint from the cached static field.
+ bool gaussKronrodPointIsCached = gaussKronrodPoint != null && gaussKronrodPoint.Order == order;
+ if (!gaussKronrodPointIsCached)
+ {
+ // Try to find the GaussKronrodPoint in the precomputed dictionary.
+ if (!GaussKronrodPoint.PreComputed.TryGetValue(order, out gaussKronrodPoint))
+ {
+ //Not yet supported!
+ gaussKronrodPoint = GaussKronrodPoint.Generate(order, targetAbsoluteTolerance);
+ }
+ }
+
+ return gaussKronrodPoint;
+ }
+ }
+}
diff --git a/src/Numerics/Integration/GaussRule/GaussPoints.cs b/src/Numerics/Integration/GaussRule/GaussPoints.cs
new file mode 100644
index 00000000..309a7247
--- /dev/null
+++ b/src/Numerics/Integration/GaussRule/GaussPoints.cs
@@ -0,0 +1,35 @@
+namespace MathNet.Numerics.Integration.GaussRule
+{
+ ///
+ /// Contains two set of abscissas, weights, and order.
+ ///
+ internal class GaussPoints
+ {
+ internal int Order { get; private set; }
+
+ internal double[] Abscissas { get; private set; }
+
+ internal double[] Weights { get; private set; }
+
+ internal int SecondOrder { get; private set; }
+
+ internal double[] SecondAbscissas { get; private set; }
+
+ internal double[] SecondWeights { get; private set; }
+
+ internal GaussPoints(int order, double[] abscissas, double[] weights, int secondOrder, double[] secondAbscissas, double[] secondWeights)
+ {
+ Order = order;
+ Abscissas = abscissas;
+ Weights = weights;
+ SecondOrder = secondOrder;
+ SecondAbscissas = secondAbscissas;
+ SecondWeights = secondWeights;
+
+ }
+
+ internal GaussPoints(int order, double[] abscissas, double[] weights, int secondOrder, double[] secondWeights)
+ : this(order, abscissas, weights, secondOrder, null, secondWeights)
+ { }
+ }
+}