@ -43,31 +43,31 @@ namespace MathNet.Numerics.Integration
{
public class GaussKronrodRule
{
private readonly GaussPointPair gaussKronrodPoint ;
readonly GaussPointPair _ gaussKronrodPoint ;
/// <summary>
/// Getter for the order.
/// </summary>
public int Order = > gaussKronrodPoint . Order ;
public int Order = > _ gaussKronrodPoint . Order ;
/// <summary>
/// Getter that returns a clone of the array containing the Kronrod abscissas.
/// </summary>
public double [ ] KronrodAbscissas = > gaussKronrodPoint . Abscissas . Clone ( ) as double [ ] ;
public double [ ] KronrodAbscissas = > _ gaussKronrodPoint . Abscissas . Clone ( ) as double [ ] ;
/// <summary>
/// Getter that returns a clone of the array containing the Kronrod weights.
/// </summary>
public double [ ] KronrodWeights = > gaussKronrodPoint . Weights . Clone ( ) as double [ ] ;
public double [ ] KronrodWeights = > _ gaussKronrodPoint . Weights . Clone ( ) as double [ ] ;
/// <summary>
/// Getter that returns a clone of the array containing the Gauss weights.
/// </summary>
public double [ ] GaussWeights = > gaussKronrodPoint . SecondWeights . Clone ( ) as double [ ] ;
public double [ ] GaussWeights = > _ gaussKronrodPoint . SecondWeights . Clone ( ) as double [ ] ;
public GaussKronrodRule ( int order )
{
gaussKronrodPoint = GaussKronrodPointFactory . GetGaussPoint ( order ) ;
_ gaussKronrodPoint = GaussKronrodPointFactory . GetGaussPoint ( order ) ;
}
/// <summary>
@ -106,10 +106,7 @@ namespace MathNet.Numerics.Integration
// 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 ) ) ;
} ;
Func < double , double > u = ( t ) = > 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]
@ -120,10 +117,7 @@ namespace MathNet.Numerics.Integration
// 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 ) ) ;
} ;
Func < double , double > u = ( s ) = > 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]
@ -134,10 +128,7 @@ namespace MathNet.Numerics.Integration
// 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 ) ) ;
} ;
Func < double , double > u = ( s ) = > - 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]
@ -147,10 +138,7 @@ namespace MathNet.Numerics.Integration
// 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 ) ;
} ;
Func < double , double > u = ( t ) = > 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 ) ;
}
}
@ -192,185 +180,174 @@ namespace MathNet.Numerics.Integration
// 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 ) ) ;
} ;
Func < double , Complex > u = ( t ) = > 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 )
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 ) ) ;
} ;
Func < double , Complex > u = ( s ) = > 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 )
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 ) ) ;
} ;
Func < double , Complex > u = ( s ) = > - 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 ) ;
} ;
Func < double , Complex > u = ( t ) = > 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 )
static double integrate_non_adaptive_m1_1 ( Func < double , double > f , out double error , out double pL1 , GaussPointPair gaussKronrodPoint )
{
int gauss_s tart = 2 ;
int kronrod_s tart = 1 ;
int gauss_o rder = ( gaussKronrodPoint . Order - 1 ) / 2 ;
int gaussS tart = 2 ;
int kronrodS tart = 1 ;
int gaussO rder = ( gaussKronrodPoint . Order - 1 ) / 2 ;
double kronrod_r esult = 0d ;
double gauss_r esult = 0d ;
double kronrodR esult = 0d ;
double gaussR esult = 0d ;
double fp , fm ;
var KAbscissa = gaussKronrodPoint . Abscissas ;
var KWeights = gaussKronrodPoint . Weights ;
var GWeights = gaussKronrodPoint . SecondWeights ;
if ( ( gauss_o rder & 1 ) = = 1 )
if ( ( gaussO rder & 1 ) = = 1 )
{
fp = f ( 0 ) ;
kronrod_r esult = fp * KWeights [ 0 ] ;
gauss_r esult + = fp * GWeights [ 0 ] ;
kronrodR esult = fp * KWeights [ 0 ] ;
gaussR esult + = fp * GWeights [ 0 ] ;
}
else
{
fp = f ( 0 ) ;
kronrod_r esult = fp * KWeights [ 0 ] ;
gauss_s tart = 1 ;
kronrod_s tart = 2 ;
kronrodR esult = fp * KWeights [ 0 ] ;
gaussS tart = 1 ;
kronrodS tart = 2 ;
}
double L1 = Math . Abs ( kronrod_r esult ) ;
double L1 = Math . Abs ( kronrodR esult ) ;
for ( int i = gauss_s tart ; i < KAbscissa . Length ; i + = 2 )
for ( int i = gaussS tart ; i < KAbscissa . Length ; i + = 2 )
{
fp = f ( KAbscissa [ i ] ) ;
fm = f ( - KAbscissa [ i ] ) ;
kronrod_r esult + = ( fp + fm ) * KWeights [ i ] ;
kronrodR esult + = ( fp + fm ) * KWeights [ i ] ;
L1 + = ( Math . Abs ( fp ) + Math . Abs ( fm ) ) * KWeights [ i ] ;
gauss_r esult + = ( fp + fm ) * GWeights [ i / 2 ] ;
gaussR esult + = ( fp + fm ) * GWeights [ i / 2 ] ;
}
for ( int i = kronrod_s tart ; i < KAbscissa . Length ; i + = 2 )
for ( int i = kronrodS tart ; i < KAbscissa . Length ; i + = 2 )
{
fp = f ( KAbscissa [ i ] ) ;
fm = f ( - KAbscissa [ i ] ) ;
kronrod_r esult + = ( fp + fm ) * KWeights [ i ] ;
kronrodR esult + = ( fp + fm ) * KWeights [ i ] ;
L1 + = ( Math . Abs ( fp ) + Math . Abs ( fm ) ) * KWeights [ i ] ;
}
pL1 = L1 ;
error = Math . Max ( Math . Abs ( kronrod_result - gauss_r esult ) , Math . Abs ( kronrod_r esult * Precision . MachineEpsilon * 2d ) ) ;
return kronrod_r esult ;
error = Math . Max ( Math . Abs ( kronrodResult - gaussR esult ) , Math . Abs ( kronrodR esult * Precision . MachineEpsilon * 2d ) ) ;
return kronrodR esult ;
}
private static Complex contour_integrate_non_adaptive_m1_1 ( Func < double , Complex > f , out double error , out double pL1 , GaussPointPair gaussKronrodPoint )
static Complex contour_integrate_non_adaptive_m1_1 ( Func < double , Complex > f , out double error , out double pL1 , GaussPointPair gaussKronrodPoint )
{
int gauss_s tart = 2 ;
int kronrod_s tart = 1 ;
int gauss_o rder = ( gaussKronrodPoint . Order - 1 ) / 2 ;
int gaussS tart = 2 ;
int kronrodS tart = 1 ;
int gaussO rder = ( gaussKronrodPoint . Order - 1 ) / 2 ;
Complex kronrod_result = new Complex ( ) ;
Complex gauss_r esult = new Complex ( ) ;
Complex kronrodResult ;
Complex gaussR esult = new Complex ( ) ;
Complex fp , fm ;
var KAbscissa = gaussKronrodPoint . Abscissas ;
var KWeights = gaussKronrodPoint . Weights ;
var GWeights = gaussKronrodPoint . SecondWeights ;
if ( gauss_o rder . IsOdd ( ) )
if ( gaussO rder . IsOdd ( ) )
{
fp = f ( 0 ) ;
kronrod_r esult = fp * KWeights [ 0 ] ;
gauss_r esult + = fp * GWeights [ 0 ] ;
kronrodR esult = fp * KWeights [ 0 ] ;
gaussR esult + = fp * GWeights [ 0 ] ;
}
else
{
fp = f ( 0 ) ;
kronrod_r esult = fp * KWeights [ 0 ] ;
gauss_s tart = 1 ;
kronrod_s tart = 2 ;
kronrodR esult = fp * KWeights [ 0 ] ;
gaussS tart = 1 ;
kronrodS tart = 2 ;
}
double L1 = Complex . Abs ( kronrod_r esult ) ;
double L1 = Complex . Abs ( kronrodR esult ) ;
for ( int i = gauss_s tart ; i < KAbscissa . Length ; i + = 2 )
for ( int i = gaussS tart ; i < KAbscissa . Length ; i + = 2 )
{
fp = f ( KAbscissa [ i ] ) ;
fm = f ( - KAbscissa [ i ] ) ;
kronrod_r esult + = ( fp + fm ) * KWeights [ i ] ;
kronrodR esult + = ( fp + fm ) * KWeights [ i ] ;
L1 + = ( Complex . Abs ( fp ) + Complex . Abs ( fm ) ) * KWeights [ i ] ;
gauss_r esult + = ( fp + fm ) * GWeights [ i / 2 ] ;
gaussR esult + = ( fp + fm ) * GWeights [ i / 2 ] ;
}
for ( int i = kronrod_s tart ; i < KAbscissa . Length ; i + = 2 )
for ( int i = kronrodS tart ; i < KAbscissa . Length ; i + = 2 )
{
fp = f ( KAbscissa [ i ] ) ;
fm = f ( - KAbscissa [ i ] ) ;
kronrod_r esult + = ( fp + fm ) * KWeights [ i ] ;
kronrodR esult + = ( fp + fm ) * KWeights [ i ] ;
L1 + = ( Complex . Abs ( fp ) + Complex . Abs ( fm ) ) * KWeights [ i ] ;
}
pL1 = L1 ;
error = Math . Max ( Complex . Abs ( kronrod_result - gauss_r esult ) , Complex . Abs ( kronrod_r esult * Precision . MachineEpsilon * 2d ) ) ;
return kronrod_r esult ;
error = Math . Max ( Complex . Abs ( kronrodResult - gaussR esult ) , Complex . Abs ( kronrodR esult * Precision . MachineEpsilon * 2d ) ) ;
return kronrodR esult ;
}
private static double recursive_adaptive_integrate ( Func < double , double > f , double a , double b , int max_levels , double rel_tol , double abs_t ol , out double error , out double L1 , GaussPointPair gaussKronrodPoint )
static double recursive_adaptive_integrate ( Func < double , double > f , double a , double b , int maxLevels , double relTol , double absT ol , 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_l ocal, out L1 , gaussKronrodPoint ) ;
var r1 = integrate_non_adaptive_m1_1 ( ( x ) = > f ( scale * x + mean ) , out double errorL ocal, out L1 , gaussKronrodPoint ) ;
var estimate = scale * r1 ;
var tmp = estimate * rel_t ol ;
var abs_t ol1 = Math . Abs ( tmp ) ;
if ( abs_t ol = = 0 )
var tmp = estimate * relT ol ;
var absT ol1 = Math . Abs ( tmp ) ;
if ( absT ol = = 0 )
{
abs_tol = abs_t ol1 ;
absTol = absT ol1 ;
}
if ( max_l evels > 0 & & ( abs_tol1 < error_l ocal ) & & ( abs_tol < error_l ocal ) )
if ( maxL evels > 0 & & ( absTol1 < errorL ocal ) & & ( absTol < errorL ocal ) )
{
double mid = ( a + b ) / 2d ;
double L1_local ;
estimate = recursive_adaptive_integrate ( f , a , mid , max_l evels - 1 , rel_tol , abs_t ol / 2 , out error , out L1 , gaussKronrodPoint ) ;
estimate + = recursive_adaptive_integrate ( f , mid , b , max_l evels - 1 , rel_tol , abs_t ol / 2 , out error_l ocal , out L1_local , gaussKronrodPoint ) ;
error + = error_l ocal ;
estimate = recursive_adaptive_integrate ( f , a , mid , maxL evels - 1 , relTol , absT ol / 2 , out error , out L1 , gaussKronrodPoint ) ;
estimate + = recursive_adaptive_integrate ( f , mid , b , maxL evels - 1 , relTol , absT ol / 2 , out errorL ocal , out L1_local , gaussKronrodPoint ) ;
error + = errorL ocal ;
L1 + = L1_local ;
return estimate ;
}
L1 * = scale ;
error = error_l ocal ;
error = errorL ocal ;
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_t ol , out double error , out double L1 , GaussPointPair gaussKronrodPoint )
static Complex contour_recursive_adaptive_integrate ( Func < double , Complex > f , double a , double b , int maxLevels , double relTol , double absT ol , out double error , out double L1 , GaussPointPair gaussKronrodPoint )
{
double error_local ;
double mean = ( b + a ) / 2 ;
@ -379,19 +356,19 @@ namespace MathNet.Numerics.Integration
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_t ol ;
var abs_t ol1 = Complex . Abs ( tmp ) ;
if ( abs_t ol = = 0 )
var tmp = estimate * relT ol ;
var absT ol1 = Complex . Abs ( tmp ) ;
if ( absT ol = = 0 )
{
abs_tol = abs_t ol1 ;
absTol = absT ol1 ;
}
if ( max_l evels > 0 & & ( abs_t ol1 < error_local ) & & ( abs_t ol < error_local ) )
if ( maxL evels > 0 & & ( absT ol1 < error_local ) & & ( absT ol < error_local ) )
{
double mid = ( a + b ) / 2d ;
double L1_local ;
estimate = contour_recursive_adaptive_integrate ( f , a , mid , max_l evels - 1 , rel_tol , abs_t ol / 2 , out error , out L1 , gaussKronrodPoint ) ;
estimate + = contour_recursive_adaptive_integrate ( f , mid , b , max_l evels - 1 , rel_tol , abs_t ol / 2 , out error_local , out L1_local , gaussKronrodPoint ) ;
estimate = contour_recursive_adaptive_integrate ( f , a , mid , maxL evels - 1 , relTol , absT ol / 2 , out error , out L1 , gaussKronrodPoint ) ;
estimate + = contour_recursive_adaptive_integrate ( f , mid , b , maxL evels - 1 , relTol , absT ol / 2 , out error_local , out L1_local , gaussKronrodPoint ) ;
error + = error_local ;
L1 + = L1_local ;
return estimate ;
Christoph Ruegg
5 years ago