forked from tsai/mathnet-numerics
Christoph Ruegg
7 years ago
9 changed files with 2130 additions and 6 deletions
@ -0,0 +1,428 @@ |
|||||
|
// <copyright file="GaussKronrodRule.cs" company="Math.NET">
|
||||
|
// Math.NET Numerics, part of the Math.NET Project
|
||||
|
// http://numerics.mathdotnet.com
|
||||
|
// http://github.com/mathnet/mathnet-numerics
|
||||
|
//
|
||||
|
// Copyright (c) 2009-2019 Math.NET
|
||||
|
//
|
||||
|
// Permission is hereby granted, free of charge, to any person
|
||||
|
// obtaining a copy of this software and associated documentation
|
||||
|
// files (the "Software"), to deal in the Software without
|
||||
|
// restriction, including without limitation the rights to use,
|
||||
|
// copy, modify, merge, publish, distribute, sublicense, and/or sell
|
||||
|
// copies of the Software, and to permit persons to whom the
|
||||
|
// Software is furnished to do so, subject to the following
|
||||
|
// conditions:
|
||||
|
//
|
||||
|
// The above copyright notice and this permission notice shall be
|
||||
|
// included in all copies or substantial portions of the Software.
|
||||
|
//
|
||||
|
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
||||
|
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
|
||||
|
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
|
||||
|
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
|
||||
|
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
|
||||
|
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
|
||||
|
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
|
||||
|
// OTHER DEALINGS IN THE SOFTWARE.
|
||||
|
// </copyright>
|
||||
|
|
||||
|
// This file uses code from the Boost Project.
|
||||
|
// Copyright John Maddock 2017.
|
||||
|
// Copyright Nick Thompson 2017.
|
||||
|
// Use, modification and distribution are subject to the
|
||||
|
// Boost Software License, Version 1.0. (See accompanying file
|
||||
|
// 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 class GaussKronrodRule |
||||
|
{ |
||||
|
private readonly GaussPointPair gaussKronrodPoint; |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Getter for the order.
|
||||
|
/// </summary>
|
||||
|
public int Order |
||||
|
{ |
||||
|
get |
||||
|
{ |
||||
|
return gaussKronrodPoint.Order; |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Getter that returns a clone of the array containing the Kronrod abscissas.
|
||||
|
/// </summary>
|
||||
|
public double[] KronrodAbscissas |
||||
|
{ |
||||
|
get |
||||
|
{ |
||||
|
return gaussKronrodPoint.Abscissas.Clone() as double[]; |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Getter that returns a clone of the array containing the Kronrod weights.
|
||||
|
/// </summary>
|
||||
|
public double[] KronrodWeights |
||||
|
{ |
||||
|
get |
||||
|
{ |
||||
|
return gaussKronrodPoint.Weights.Clone() as double[]; |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Getter that returns a clone of the array containing the Gauss weights.
|
||||
|
/// </summary>
|
||||
|
public double[] GaussWeights |
||||
|
{ |
||||
|
get |
||||
|
{ |
||||
|
return gaussKronrodPoint.SecondWeights.Clone() as double[]; |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
public GaussKronrodRule(int order) |
||||
|
{ |
||||
|
gaussKronrodPoint = GaussKronrodPointFactory.GetGaussPoint(order); |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Performs adaptive Gauss-Kronrod quadrature on function f over the range (a,b)
|
||||
|
/// </summary>
|
||||
|
/// <param name="f">The analytic smooth function to integrate</param>
|
||||
|
/// <param name="intervalBegin">Where the interval starts</param>
|
||||
|
/// <param name="intervalEnd">Where the interval stops</param>
|
||||
|
/// <param name="error">The difference between the (N-1)/2 point Gauss approximation and the N-point Gauss-Kronrod approximation</param>
|
||||
|
/// <param name="L1Norm">The L1 norm of the result, if there is a significant difference between this and the returned value, then the result is likely to be ill-conditioned.</param>
|
||||
|
/// <param name="targetRelativeError">The maximum relative error in the result</param>
|
||||
|
/// <param name="maximumDepth">The maximum number of interval splittings permitted before stopping</param>
|
||||
|
/// <param name="order">The number of Gauss-Kronrod points. Pre-computed for 15, 21, 31, 41, 51 and 61 points</param>
|
||||
|
public static double Integrate(Func<double, double> f, double intervalBegin, double intervalEnd, out double error, out double L1Norm, double targetRelativeError = 1E-10, int maximumDepth = 15, int order = 15) |
||||
|
{ |
||||
|
// Formula used for variable subsitution from
|
||||
|
// 1. Shampine, L. F. (2008). Vectorized adaptive quadrature in MATLAB. Journal of Computational and Applied Mathematics, 211(2), 131-140.
|
||||
|
// 2. quadgk.m, GNU Octave
|
||||
|
|
||||
|
if (f == null) |
||||
|
{ |
||||
|
throw new ArgumentNullException(nameof(f)); |
||||
|
} |
||||
|
|
||||
|
if (intervalBegin > intervalEnd) |
||||
|
{ |
||||
|
return -Integrate(f, intervalEnd, intervalBegin, out error, out L1Norm, targetRelativeError, maximumDepth, order); |
||||
|
} |
||||
|
|
||||
|
GaussPointPair gaussKronrodPoint = GaussKronrodPointFactory.GetGaussPoint(order); |
||||
|
|
||||
|
// (-oo, oo) => [-1, 1]
|
||||
|
//
|
||||
|
// integral_(-oo)^(oo) f(x) dx = integral_(-1)^(1) f(g(t)) g'(t) dt
|
||||
|
// g(t) = t / (1 - t^2)
|
||||
|
// g'(t) = (1 + t^2) / (1 - t^2)^2
|
||||
|
if ((intervalBegin < double.MinValue) && (intervalEnd > double.MaxValue)) |
||||
|
{ |
||||
|
Func<double, double> u = (t) => |
||||
|
{ |
||||
|
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, gaussKronrodPoint); |
||||
|
} |
||||
|
// [a, oo) => [0, 1]
|
||||
|
//
|
||||
|
// integral_(a)^(oo) f(x) dx = integral_(0)^(oo) f(a + t^2) 2 t dt
|
||||
|
// = integral_(0)^(1) f(a + g(s)^2) 2 g(s) g'(s) ds
|
||||
|
// g(s) = s / (1 - s)
|
||||
|
// g'(s) = 1 / (1 - s)^2
|
||||
|
else if (intervalEnd > double.MaxValue) |
||||
|
{ |
||||
|
Func<double, double> u = (s) => |
||||
|
{ |
||||
|
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, gaussKronrodPoint); |
||||
|
} |
||||
|
// (-oo, b] => [-1, 0]
|
||||
|
//
|
||||
|
// integral_(-oo)^(b) f(x) dx = -integral_(-oo)^(0) f(b - t^2) 2 t dt
|
||||
|
// = -integral_(-1)^(0) f(b - g(s)^2) 2 g(s) g'(s) ds
|
||||
|
// g(s) = s / (1 + s)
|
||||
|
// g'(s) = 1 / (1 + s)^2
|
||||
|
else if (intervalBegin < double.MinValue) |
||||
|
{ |
||||
|
Func<double, double> u = (s) => |
||||
|
{ |
||||
|
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, gaussKronrodPoint); |
||||
|
} |
||||
|
// [a, b] => [-1, 1]
|
||||
|
//
|
||||
|
// integral_(a)^(b) f(x) dx = integral_(-1)^(1) f(g(t)) g'(t) dt
|
||||
|
// g(t) = (b - a) * t * (3 - t^2) / 4 + (b + a) / 2
|
||||
|
// g'(t) = 3 / 4 * (b - a) * (1 - t^2)
|
||||
|
else |
||||
|
{ |
||||
|
Func<double, double> u = (t) => |
||||
|
{ |
||||
|
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, gaussKronrodPoint); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Performs adaptive Gauss-Kronrod quadrature on function f over the range (a,b)
|
||||
|
/// </summary>
|
||||
|
/// <param name="f">The analytic smooth complex function to integrate, defined on the real axis.</param>
|
||||
|
/// <param name="intervalBegin">Where the interval starts</param>
|
||||
|
/// <param name="intervalEnd">Where the interval stops</param>
|
||||
|
/// <param name="error">The difference between the (N-1)/2 point Gauss approximation and the N-point Gauss-Kronrod approximation</param>
|
||||
|
/// <param name="L1Norm">The L1 norm of the result, if there is a significant difference between this and the returned value, then the result is likely to be ill-conditioned.</param>
|
||||
|
/// <param name="targetRelativeError">The maximum relative error in the result</param>
|
||||
|
/// <param name="maximumDepth">The maximum number of interval splittings permitted before stopping</param>
|
||||
|
/// <param name="order">The number of Gauss-Kronrod points. Pre-computed for 15, 21, 31, 41, 51 and 61 points</param>
|
||||
|
/// <returns></returns>
|
||||
|
public static Complex ContourIntegrate(Func<double, Complex> f, double intervalBegin, double intervalEnd, out double error, out double L1Norm, double targetRelativeError = 1E-10, int maximumDepth = 15, int order = 15) |
||||
|
{ |
||||
|
// Formula used for variable subsitution from
|
||||
|
// 1. Shampine, L. F. (2008). Vectorized adaptive quadrature in MATLAB. Journal of Computational and Applied Mathematics, 211(2), 131-140.
|
||||
|
// 2. quadgk.m, GNU Octave
|
||||
|
|
||||
|
if (f == null) |
||||
|
{ |
||||
|
throw new ArgumentNullException(nameof(f)); |
||||
|
} |
||||
|
|
||||
|
if (intervalBegin > intervalEnd) |
||||
|
{ |
||||
|
return -ContourIntegrate(f, intervalEnd, intervalBegin, out error, out L1Norm, targetRelativeError, maximumDepth, order); |
||||
|
} |
||||
|
|
||||
|
GaussPointPair gaussKronrodPoint = GaussKronrodPointFactory.GetGaussPoint(order); |
||||
|
|
||||
|
// (-oo, oo) => [-1, 1]
|
||||
|
//
|
||||
|
// integral_(-oo)^(oo) f(x) dx = integral_(-1)^(1) f(g(t)) g'(t) dt
|
||||
|
// g(t) = t / (1 - t^2)
|
||||
|
// g'(t) = (1 + t^2) / (1 - t^2)^2
|
||||
|
if ((intervalBegin < double.MinValue) && (intervalEnd > double.MaxValue)) |
||||
|
{ |
||||
|
Func<double, Complex> u = (t) => |
||||
|
{ |
||||
|
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, gaussKronrodPoint); |
||||
|
} |
||||
|
// [a, oo) => [0, 1]
|
||||
|
//
|
||||
|
// integral_(a)^(oo) f(x) dx = integral_(0)^(oo) f(a + t^2) 2 t dt
|
||||
|
// = integral_(0)^(1) f(a + g(s)^2) 2 g(s) g'(s) ds
|
||||
|
// g(s) = s / (1 - s)
|
||||
|
// g'(s) = 1 / (1 - s)^2
|
||||
|
else if (intervalEnd > double.MaxValue) |
||||
|
{ |
||||
|
Func<double, Complex> u = (s) => |
||||
|
{ |
||||
|
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, gaussKronrodPoint); |
||||
|
} |
||||
|
// (-oo, b] => [-1, 0]
|
||||
|
//
|
||||
|
// integral_(-oo)^(b) f(x) dx = -integral_(-oo)^(0) f(b - t^2) 2 t dt
|
||||
|
// = -integral_(-1)^(0) f(b - g(s)^2) 2 g(s) g'(s) ds
|
||||
|
// g(s) = s / (1 + s)
|
||||
|
// g'(s) = 1 / (1 + s)^2
|
||||
|
else if (intervalBegin < double.MinValue) |
||||
|
{ |
||||
|
Func<double, Complex> u = (s) => |
||||
|
{ |
||||
|
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, gaussKronrodPoint); |
||||
|
} |
||||
|
// [a, b] => [-1, 1]
|
||||
|
//
|
||||
|
// integral_(a)^(b) f(x) dx = integral_(-1)^(1) f(g(t)) g'(t) dt
|
||||
|
// g(t) = (b - a) * t * (3 - t^2) / 4 + (b + a) / 2
|
||||
|
// g'(t) = 3 / 4 * (b - a) * (1 - t^2)
|
||||
|
else |
||||
|
{ |
||||
|
Func<double, Complex> u = (t) => |
||||
|
{ |
||||
|
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, gaussKronrodPoint); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
private static double integrate_non_adaptive_m1_1(Func<double, double> f, out double error, out double pL1, GaussPointPair gaussKronrodPoint) |
||||
|
{ |
||||
|
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) |
||||
|
{ |
||||
|
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) |
||||
|
{ |
||||
|
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) |
||||
|
{ |
||||
|
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; |
||||
|
} |
||||
|
|
||||
|
private static Complex contour_integrate_non_adaptive_m1_1(Func<double, Complex> f, out double error, out double pL1, GaussPointPair gaussKronrodPoint) |
||||
|
{ |
||||
|
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()) |
||||
|
{ |
||||
|
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 = 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 * Precision.MachineEpsilon * 2d)); |
||||
|
return kronrod_result; |
||||
|
} |
||||
|
|
||||
|
private static double recursive_adaptive_integrate(Func<double, double> f, double a, double b, int max_levels, double rel_tol, double abs_tol, out double error, out double L1, GaussPointPair gaussKronrodPoint) |
||||
|
{ |
||||
|
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) |
||||
|
{ |
||||
|
abs_tol = abs_tol1; |
||||
|
} |
||||
|
|
||||
|
if (max_levels > 0 && (abs_tol1 < error_local) && (abs_tol < error_local)) |
||||
|
{ |
||||
|
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; |
||||
|
} |
||||
|
|
||||
|
private static Complex contour_recursive_adaptive_integrate(Func<double, Complex> f, double a, double b, int max_levels, double rel_tol, double abs_tol, out double error, out double L1, GaussPointPair gaussKronrodPoint) |
||||
|
{ |
||||
|
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) |
||||
|
{ |
||||
|
abs_tol = abs_tol1; |
||||
|
} |
||||
|
|
||||
|
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; |
||||
|
} |
||||
|
} |
||||
|
} |
||||
@ -0,0 +1,618 @@ |
|||||
|
using System; |
||||
|
using System.Collections.Generic; |
||||
|
using System.Linq; |
||||
|
using System.Numerics; |
||||
|
|
||||
|
namespace MathNet.Numerics.Integration.GaussRule |
||||
|
{ |
||||
|
/// <summary>
|
||||
|
/// Contains a method to compute the Gauss-Kronrod abscissas/weights and precomputed abscissas/weights for orders 15, 21, 31, 41, 51, 61.
|
||||
|
/// </summary>
|
||||
|
internal static partial class GaussKronrodPoint |
||||
|
{ |
||||
|
/// <summary>
|
||||
|
/// Precomputed abscissas/weights for orders 15, 21, 31, 41, 51, 61.
|
||||
|
/// </summary>
|
||||
|
internal static readonly Dictionary<int, GaussPointPair> PreComputed = new Dictionary<int, GaussPointPair> |
||||
|
{ |
||||
|
{ 15, new GaussPointPair(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 GaussPointPair(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 GaussPointPair(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 GaussPointPair(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 GaussPointPair(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 GaussPointPair(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, |
||||
|
}) |
||||
|
}, |
||||
|
}; |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Contains a method to compute the Gauss-Kronrod abscissas/weights.
|
||||
|
/// </summary>
|
||||
|
internal static partial class GaussKronrodPoint |
||||
|
{ |
||||
|
/// <summary>
|
||||
|
/// Computes the Gauss-Kronrod abscissas/weights and Gauss weights.
|
||||
|
/// </summary>
|
||||
|
/// <param name="order">Defines an Nth order Gauss-Kronrod rule. The order also defines the number of abscissas and weights for the rule.</param>
|
||||
|
/// <param name="eps">Required precision to compute the abscissas/weights.</param>
|
||||
|
/// <returns>Object containing the non-negative abscissas/weights, order.</returns>
|
||||
|
internal static GaussPointPair Generate(int order, double eps) |
||||
|
{ |
||||
|
int gaussOrder = (order - 1) / 2; |
||||
|
int gaussStart = gaussOrder.IsOdd() ? 0 : 1; |
||||
|
int kronrodStart = gaussOrder.IsOdd() ? 1 : 0; |
||||
|
|
||||
|
var gaussPoint = GaussLegendrePointFactory.GetGaussPoint(gaussOrder); |
||||
|
var gaussAbscissas = gaussPoint.Abscissas; |
||||
|
var gaussWeights = gaussPoint.Weights; |
||||
|
|
||||
|
// Calculate Kronrod polynomial in terms of Legendre polynomials
|
||||
|
// K(x) = c0*P(0, x) + c1*P(1, x) + ...
|
||||
|
|
||||
|
var c = StieltjesP(gaussOrder + 1); |
||||
|
|
||||
|
// Calculate Abscissas for Kronrod polynomial
|
||||
|
|
||||
|
int r = gaussOrder.IsOdd() ? (gaussOrder - 1) / 2 + 1 : gaussOrder / 2 + 1; |
||||
|
var kronrodAbscissas = new double[r]; |
||||
|
|
||||
|
for (int k = 1; k <= gaussOrder + 1; k = k + 2) |
||||
|
{ |
||||
|
var x0 = (1.0 - (1.0 - 1.0 / gaussOrder) / (8 * gaussOrder * gaussOrder)) * Math.Cos((k - 0.5) * Math.PI / (2.0 * gaussOrder + 1.0)); |
||||
|
var dx = 0d; |
||||
|
var j = 1; // iterations
|
||||
|
|
||||
|
// Newton iterations
|
||||
|
do |
||||
|
{ |
||||
|
var E = LegendreSeries(c, x0); |
||||
|
dx = E.Item1 / E.Item2; |
||||
|
x0 = x0 - dx; |
||||
|
j++; |
||||
|
} |
||||
|
while (Math.Abs(dx) > eps && j < 100); |
||||
|
|
||||
|
if (Math.Abs(x0) < Precision.MachineEpsilon) x0 = 0.0; |
||||
|
|
||||
|
kronrodAbscissas[(k - 1) / 2] = x0; |
||||
|
} |
||||
|
|
||||
|
// Concatenate two abscissas
|
||||
|
|
||||
|
var abscissas = new double[gaussAbscissas.Length + kronrodAbscissas.Length]; |
||||
|
gaussAbscissas.CopyTo(abscissas, 0); |
||||
|
kronrodAbscissas.CopyTo(abscissas, gaussAbscissas.Length); |
||||
|
abscissas = abscissas.OrderBy(v => v).ToArray(); |
||||
|
|
||||
|
// Calculate weights for abscissas
|
||||
|
|
||||
|
var weights = new double[gaussAbscissas.Length + kronrodAbscissas.Length]; |
||||
|
for (int i = gaussStart; i < abscissas.Length; i += 2) |
||||
|
{ |
||||
|
var x = abscissas[i]; |
||||
|
|
||||
|
var E = LegendreSeries(c, x); |
||||
|
var L = LegendreP(gaussOrder, x); |
||||
|
|
||||
|
var p = L.Item2; |
||||
|
var w2 = 2.0 / ((1.0 - x * x) * p * p); // Gauss weight
|
||||
|
weights[i] = w2 + 2.0 / ((gaussOrder + 1.0) * p * E.Item1); |
||||
|
} |
||||
|
for (int i = kronrodStart; i < abscissas.Length; i += 2) |
||||
|
{ |
||||
|
var x = abscissas[i]; |
||||
|
|
||||
|
var E = LegendreSeries(c, x); |
||||
|
var L = LegendreP(gaussOrder, x); |
||||
|
|
||||
|
weights[i] = 2.0 / ((gaussOrder + 1.0) * L.Item1 * E.Item2); |
||||
|
} |
||||
|
|
||||
|
return new GaussPointPair(order, abscissas, weights, gaussOrder, gaussWeights); |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Returns coefficients of a Stieltjes polynomial in terms of Legendre polynomials.
|
||||
|
/// </summary>
|
||||
|
internal static double[] StieltjesP(int order) |
||||
|
{ |
||||
|
// Reference:
|
||||
|
// 1. Patterson, Thomas NL. "The optimum addition of points to quadrature formulae." Mathematics of Computation 22.104 (1968): 847-856.
|
||||
|
// 2. Piessens, Robert, and Maria Branders. "A note on the optimal addition of abscissas to quadrature formulas of Gauss and Lobatto type." Mathematics of Computation (1974): 135-139.
|
||||
|
// 3. Legendre-Stieltjes Polynomials, Boost.Math
|
||||
|
//
|
||||
|
// Here, we are using Patterson algorithm, expanding the Stieltjes polynomial in terms of Legendre polynomials.
|
||||
|
//
|
||||
|
// Kronrod Polynomial K[n + 1, x] is expanded in terms of Legendre Polynomial P[n, x].
|
||||
|
//
|
||||
|
// K[n + 1, x] = sum_(n=1)^r a[i] P[2 * i - 1 - q, x]
|
||||
|
//
|
||||
|
// where P[n, x] is the Legendre polynomial of degree n,
|
||||
|
// [x] denotes the integer part of x,
|
||||
|
// q = n - 2[n/2]
|
||||
|
// r = [(n + 3)/2]
|
||||
|
//
|
||||
|
// The added n + 1 Kronrod abscissae is the roots of the Kronrod polynomial.
|
||||
|
|
||||
|
if (order == 1) // P(1, x)
|
||||
|
return new double[] { 0, 1 }; |
||||
|
else if (order == 2) // -2/5 * P(0, x) + P(2, x)
|
||||
|
return new double[] { -0.4, 0, 1 }; |
||||
|
else if (order == 3) // -9/14 * P(1, x) + P(3, x)
|
||||
|
return new double[] { 0, -0.642857142857142857142857142857, 0, 1 }; |
||||
|
else if (order == 4) // 14/891 * P(0, x) - 20/27 * P(2, x) + P(4, x)
|
||||
|
return new double[] { 0.0157126823793490460157126823793, 0, -0.740740740740740740740740740741, 0, 1 }; |
||||
|
else if (order == 5) // 135/12584 * P(1, x) - 35/44 * P(3, x) + P(5, x)
|
||||
|
return new double[] { 0, 0.0107279084551811824539097266370, 0, -0.795454545454545454545454545455, 0, 1 }; |
||||
|
|
||||
|
int n = order - 1; |
||||
|
int q = n.IsOdd() ? 1 : 0; |
||||
|
int r = n.IsOdd() ? (n - 1) / 2 + 2 : n / 2 + 1; |
||||
|
|
||||
|
double[] a = new double[r + 1]; |
||||
|
|
||||
|
// Calculate a[i] for i = 1, ..., r
|
||||
|
//
|
||||
|
// a[r] = 1;
|
||||
|
// a[r - 1] = -a[r] * S[r, 1] / S[r - 1, 1];
|
||||
|
// a[r - 2] = -a[r] * S[r, 2] / S[r - 2, 2] - a[r - 1] * S[r - 1, 2] / S[r - 2, 2];
|
||||
|
// ...
|
||||
|
// a[1] = -a[r] * S[r, r - 1] / S[1, r - 1] - a[r - 1] * S[r - 1, r - 1] / S[1, r - 1] - ... - a[2] * S[2, r - 1] / S[1, r - 1];
|
||||
|
//
|
||||
|
// S[i, k] / S[r - k, k] = S[i - 1, k] / S[r - k, k]
|
||||
|
// * ((n - q + 2 * (i + k - 1)) * (n + q + 2 * (k - i + 1)) * (n - 1 - q + 2 * (i - k)) * (2 * (k + i - 1) - 1 - q - n))
|
||||
|
// / ((n - q + 2 * (i - k)) * (2 * (k + i - 1) - q - n) * (n + 1 + q + 2 * (k - i)) * (n - 1 - q + 2 * (i + k)));
|
||||
|
|
||||
|
a[r] = 1.0; |
||||
|
for (int k = 1; k < r; k++) |
||||
|
{ |
||||
|
double ratio = 1.0; |
||||
|
a[r - k] = 0.0; |
||||
|
for (int i = r + 1 - k; i <= r; i++) |
||||
|
{ |
||||
|
double numerator = (n - q + 2 * (i + k - 1)) * (n + q + 2 * (k - i + 1)) * (n - 1 - q + 2 * (i - k)) * (2 * (k + i - 1) - 1 - q - n); |
||||
|
double denominator = (n - q + 2 * (i - k)) * (2 * (k + i - 1) - q - n) * (n + 1 + q + 2 * (k - i)) * (n - 1 - q + 2 * (i + k)); |
||||
|
ratio = ratio * numerator / denominator; |
||||
|
a[r - k] -= a[i] * ratio; |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
// K = sum c[k] P[k, x]
|
||||
|
|
||||
|
double[] c = new double[2 * r - q]; |
||||
|
for (int i = 1; i < a.Length; i++) |
||||
|
{ |
||||
|
c[2 * i - 1 - q] = a[i]; |
||||
|
} |
||||
|
|
||||
|
return c; |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Return value and derivative of a Legendre series at given points.
|
||||
|
/// </summary>
|
||||
|
internal static Tuple<double, double> LegendreSeries(double[] a, double x) |
||||
|
{ |
||||
|
// S = a[0]*P[0, x] + ... + a[k]*P[k, x] + ... + a[n]*P[n, x]
|
||||
|
// where P[k, x] is the Legendre polynomial of order k
|
||||
|
//
|
||||
|
// According to the Clenshaw algorithm, S can be written by
|
||||
|
// S = a[0] + x*b[1, x] - 1/2 * b[2,x]
|
||||
|
//
|
||||
|
// b[n + 1, x] = 0
|
||||
|
// b[n + 2, x] = 0
|
||||
|
// b[k, x] = a[k] + (2k + 1)/(k + 1)*x*b[k + 1, x] - (k + 1)/(k + 2)*b[k + 2, x]
|
||||
|
//
|
||||
|
// Derivative of S is given by
|
||||
|
// S' = b[1, x] + x*b'[1, x] - 1/2 * b'[2,x]
|
||||
|
//
|
||||
|
// b'[k, x] = (2k + 1)/(k + 1)*b[k + 1, x] + (2k + 1)/(k + 1)*x*b'[k + 1, x] - (k + 1)/(k + 2)*b'[k + 2, x]
|
||||
|
|
||||
|
if (a.Length == 1) |
||||
|
return new Tuple<double, double>(a[0], 0); |
||||
|
if (a.Length == 2) |
||||
|
return new Tuple<double, double>(a[0] + a[1] * x, a[1]); |
||||
|
|
||||
|
double b0 = 0.0, b1 = 0.0, b2 = 0.0; |
||||
|
double p0 = 0.0, p1 = 0.0, p2 = 0.0; |
||||
|
|
||||
|
for (int k = a.Length - 1; k >= 1; k--) |
||||
|
{ |
||||
|
b0 = a[k] + (2.0 * k + 1.0) / (k + 1.0) * x * b1 - (k + 1.0) / (k + 2.0) * b2; |
||||
|
p0 = (2.0 * k + 1.0) / (k + 1.0) * (b1 + x * p1) - (k + 1.0) / (k + 2.0) * p2; |
||||
|
|
||||
|
b2 = b1; |
||||
|
b1 = b0; |
||||
|
p2 = p1; |
||||
|
p1 = p0; |
||||
|
} |
||||
|
|
||||
|
var value = a[0] + b1 * x - 0.5 * b2; |
||||
|
var derivative = b1 + p1 * x - 0.5 * p2; |
||||
|
|
||||
|
return new Tuple<double, double>( value, derivative ); |
||||
|
} |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Return value and derivative of a Legendre polynomial of order at given points.
|
||||
|
/// </summary>
|
||||
|
internal static Tuple<double, double> LegendreP(int order, double x) |
||||
|
{ |
||||
|
// The Legendre polynomial, P[n, x], is defined by the recurrence relation:
|
||||
|
//
|
||||
|
// P[0, x] = 1
|
||||
|
// P[1, x] = x
|
||||
|
// (n + 1) * P[n + 1, x] = (2 * n + 1) * x * P[n, x] - n * P[n - 1, x]
|
||||
|
//
|
||||
|
// The derivative of the Legendre polynomial, P'[n, x] is given by
|
||||
|
// P'[0, x] = 0
|
||||
|
// P'[1, x] = 1
|
||||
|
// (n + 1) * P'[n + 1, x] = (2 * n + 1) * P[n, x] + (2 * n + 1) * x * P'[n, x] - n * P'[n - 1, x]
|
||||
|
// = (2 * n + 1) * (P[n, x] + x * P'[n, x]) - n * P'[n - 1, x]
|
||||
|
|
||||
|
if (order == 0) |
||||
|
return new Tuple<double, double>(1.0, 0.0); |
||||
|
if (order == 1) |
||||
|
return new Tuple<double, double>(x, 1.0); |
||||
|
|
||||
|
double b0 = 0.0, b1 = 1.0, b2 = 0.0; |
||||
|
double p0 = 0.0, p1 = 0.0, p2 = 0.0; |
||||
|
|
||||
|
for (int k = 1; k <= order; k++) |
||||
|
{ |
||||
|
b0 = (2.0 * k - 1.0) / k * x * b1 - (k - 1.0) / k * b2; // L(k, x)
|
||||
|
p0 = (2.0 * k - 1.0) / k * (b1 + x * p1) - (k - 1.0) / k * p2; // L'(k, x)
|
||||
|
|
||||
|
b2 = b1; |
||||
|
b1 = b0; |
||||
|
p2 = p1; |
||||
|
p1 = p0; |
||||
|
} |
||||
|
|
||||
|
var value = b0; |
||||
|
var derivative = p0; |
||||
|
|
||||
|
return new Tuple<double, double>(value, derivative); |
||||
|
} |
||||
|
} |
||||
|
} |
||||
@ -0,0 +1,34 @@ |
|||||
|
using System; |
||||
|
|
||||
|
namespace MathNet.Numerics.Integration.GaussRule |
||||
|
{ |
||||
|
/// <summary>
|
||||
|
/// Creates a Gauss-Kronrod point.
|
||||
|
/// </summary>
|
||||
|
internal static class GaussKronrodPointFactory |
||||
|
{ |
||||
|
[ThreadStatic] |
||||
|
private static GaussPointPair gaussKronrodPoint; |
||||
|
|
||||
|
/// <summary>
|
||||
|
/// Getter for the GaussKronrodPoint.
|
||||
|
/// </summary>
|
||||
|
/// <param name="order">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.</param>
|
||||
|
/// <returns>Object containing the non-negative abscissas/weights, and order.</returns>
|
||||
|
public static GaussPointPair GetGaussPoint(int order) |
||||
|
{ |
||||
|
// 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)) |
||||
|
{ |
||||
|
gaussKronrodPoint = GaussKronrodPoint.Generate(order, 1E-10); |
||||
|
} |
||||
|
} |
||||
|
|
||||
|
return gaussKronrodPoint; |
||||
|
} |
||||
|
} |
||||
|
} |
||||
@ -0,0 +1,40 @@ |
|||||
|
namespace MathNet.Numerics.Integration.GaussRule |
||||
|
{ |
||||
|
/// <summary>
|
||||
|
/// Contains two GaussPoint.
|
||||
|
/// </summary>
|
||||
|
internal class GaussPointPair |
||||
|
{ |
||||
|
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 double IntervalBegin { get; private set; } |
||||
|
|
||||
|
internal double IntervalEnd { get; private set; } |
||||
|
|
||||
|
internal GaussPointPair(double intervalBegin, double intervalEnd, int order, double[] abscissas, double[] weights, int secondOrder, double[] secondAbscissas, double[] secondWeights) |
||||
|
{ |
||||
|
IntervalBegin = intervalBegin; |
||||
|
IntervalEnd = intervalEnd; |
||||
|
Order = order; |
||||
|
Abscissas = abscissas; |
||||
|
Weights = weights; |
||||
|
SecondOrder = secondOrder; |
||||
|
SecondAbscissas = secondAbscissas; |
||||
|
SecondWeights = secondWeights; |
||||
|
} |
||||
|
|
||||
|
internal GaussPointPair(int order, double[] abscissas, double[] weights, int secondOrder, double[] secondWeights) |
||||
|
: this(-1, 1, order, abscissas, weights, secondOrder, null, secondWeights) |
||||
|
{ } |
||||
|
} |
||||
|
} |
||||
Loading…
Reference in new issue