diff --git a/.paket/Paket.Restore.targets b/.paket/Paket.Restore.targets
index 818b4eca..952ad42f 100644
--- a/.paket/Paket.Restore.targets
+++ b/.paket/Paket.Restore.targets
@@ -5,6 +5,11 @@
$(MSBuildAllProjects);$(MSBuildThisFileFullPath)
+
+ $(MSBuildVersion)
+ 15.0.0
+ false
+ true
true
$(MSBuildThisFileDirectory)
@@ -59,7 +64,7 @@
-
+
true
true
@@ -73,37 +78,47 @@
+
+
+
-
+
true
$(NoWarn);NU1603;NU1604;NU1605;NU1608
+ false
+ true
-
-
- /usr/bin/shasum "$(PaketRestoreCacheFile)" | /usr/bin/awk '{ print $1 }'
- /usr/bin/shasum "$(PaketLockFilePath)" | /usr/bin/awk '{ print $1 }'
+
+
+
+
+
+
+ $([System.IO.File]::ReadAllText('$(PaketRestoreCacheFile)'))
+
+
+
+
+
+
+ $([System.Text.RegularExpressions.Regex]::Split(`%(Identity)`, `": "`)[0].Replace(`"`, ``).Replace(` `, ``))
+ $([System.Text.RegularExpressions.Regex]::Split(`%(Identity)`, `": "`)[1].Replace(`"`, ``).Replace(` `, ``))
+
+
+
+
+ %(PaketRestoreCachedKeyValue.Value)
+ %(PaketRestoreCachedKeyValue.Value)
-
-
-
-
-
-
-
-
-
-
-
-
-
- $([System.IO.File]::ReadAllText('$(PaketRestoreCacheFile)'))
- $([System.IO.File]::ReadAllText('$(PaketLockFilePath)'))
+
+
true
- false
+ false
true
@@ -116,7 +131,10 @@
+
+
+
@@ -126,7 +144,7 @@
-
+
$(PaketIntermediateOutputPath)\$(MSBuildProjectFile).paket.references.cached
@@ -163,6 +181,7 @@
+
@@ -224,8 +243,6 @@
false
- $(MSBuildVersion)
- 15.8.0
@@ -252,10 +269,12 @@
- <_NuspecFiles Include="$(AdjustedNuspecOutputPath)\*.$(PackageVersion).nuspec"/>
+ <_NuspecFiles Include="$(AdjustedNuspecOutputPath)\*.$(PackageVersion.Split(`+`)[0]).nuspec"/>
-
+
+
+
@@ -281,7 +300,7 @@
DevelopmentDependency="$(DevelopmentDependency)"
BuildOutputInPackage="@(_BuildOutputInPackage)"
TargetPathsToSymbols="@(_TargetPathsToSymbols)"
- SymbolPackageFormat="symbols.nupkg"
+ SymbolPackageFormat="$(SymbolPackageFormat)"
TargetFrameworks="@(_TargetFrameworks)"
AssemblyName="$(AssemblyName)"
PackageOutputPath="$(PackageOutputAbsolutePath)"
@@ -328,7 +347,7 @@
DevelopmentDependency="$(DevelopmentDependency)"
BuildOutputInPackage="@(_BuildOutputInPackage)"
TargetPathsToSymbols="@(_TargetPathsToSymbols)"
- SymbolPackageFormat="symbols.nupkg"
+ SymbolPackageFormat="$(SymbolPackageFormat)"
TargetFrameworks="@(_TargetFrameworks)"
AssemblyName="$(AssemblyName)"
PackageOutputPath="$(PackageOutputAbsolutePath)"
diff --git a/src/Numerics.Tests/IntegrationTests/IntegrationTest.cs b/src/Numerics.Tests/IntegrationTests/IntegrationTest.cs
index ad362cde..e5e55ca2 100644
--- a/src/Numerics.Tests/IntegrationTests/IntegrationTest.cs
+++ b/src/Numerics.Tests/IntegrationTests/IntegrationTest.cs
@@ -59,7 +59,7 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests
{
return Math.Exp(-x / 5) * (2 + Math.Sin(2 * y));
}
-
+
///
/// Test Function Start point.
///
@@ -311,5 +311,41 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests
GaussLegendreRule gaussLegendre = new GaussLegendreRule(StartA, StopA, order);
Assert.AreEqual(gaussLegendre.IntervalEnd, StopA);
}
+
+ // integral_(-oo)^(oo) exp(-x^2/2) dx = sqrt(2 ¥ð)
+ [TestCase(double.NegativeInfinity, double.PositiveInfinity, Constants.Sqrt2Pi)]
+ // integral_(-oo)^(0) exp(-x^2/2) dx = sqrt(¥ð/2)
+ [TestCase(double.NegativeInfinity, 0, Constants.SqrtPiOver2)]
+ // integral_(0)^(oo exp(-x^2/2) dx = sqrt(¥ð/2)
+ [TestCase(0, double.PositiveInfinity, Constants.SqrtPiOver2)]
+ // integral_(-1)^(1) exp(-x^2/2) dx = sqrt(2 ¥ð) erf(1/sqrt(2))
+ [TestCase(-1, 1, 1.7112487837842976063)]
+ // integral_(1)^(0) exp(-x^2/2) dx = -sqrt(¥ð/2) erf(1/sqrt(2))
+ [TestCase(1, 0, -0.85562439189214880317)]
+ public void TestGaussianIntegralBySubstitution(double a, double b, double expected)
+ {
+ Assert.AreEqual(
+ expected,
+ Integrate.OnOpenInterval((x) => Math.Exp(-x * x / 2), a, b),
+ 1e-10,
+ "Integral e^(-x^2 /2) from {0} to {1}", a, b);
+ }
+
+ // integral_(-oo)^(oo) sin(pi x) / (pi x) dx = 1 / pi integral_(-oo)^(oo) sin(x) / x dx
+ // = 1 / pi integral_(oo)^(oo) 1 / (1 + t^2) dt
+ // = 1
+ // or = 2 / pi integral_(0)^(oo) 1 / (1 + t^2) dt
+ // or = 2 / pi integral_(-oo)^(0) 1 / (1 + t^2) dt
+ [TestCase(double.NegativeInfinity, double.PositiveInfinity, 1, Constants.InvPi)]
+ [TestCase(0, double.PositiveInfinity, 1, Constants.TwoInvPi)]
+ [TestCase(double.NegativeInfinity, 0, 1, Constants.TwoInvPi)]
+ public void TestSincIntegralBySubstitution(double a, double b, double expected, double factor)
+ {
+ Assert.AreEqual(
+ expected,
+ factor * Integrate.OnOpenInterval((x) => 1 / (1 + x * x), a, b),
+ 1e-10,
+ "Integral sin(pi*x)/(pi*x) from -oo to oo");
+ }
}
}
diff --git a/src/Numerics/Integrate.cs b/src/Numerics/Integrate.cs
index 2d708fdf..461516e0 100644
--- a/src/Numerics/Integrate.cs
+++ b/src/Numerics/Integrate.cs
@@ -45,21 +45,85 @@ namespace MathNet.Numerics
/// Where the interval stops, inclusive and finite.
/// The expected relative accuracy of the approximation.
/// Approximation of the finite integral in the given interval.
- public static double OnClosedInterval(Func f, double intervalBegin, double intervalEnd, double targetAbsoluteError)
+ public static double OnClosedInterval(Func f, double intervalBegin, double intervalEnd, double targetAbsoluteError = 1E-8)
{
return DoubleExponentialTransformation.Integrate(f, intervalBegin, intervalEnd, targetAbsoluteError);
}
///
- /// Approximation of the definite integral of an analytic smooth function on a closed interval.
+ /// Approximation of the definite integral of an analytic smooth function by substitution. When either or both limits are infinite, the integrand is assumed rapidly decayed to zero as x -> infinity.
///
/// The analytic smooth function to integrate.
/// Where the interval starts, inclusive and finite.
/// Where the interval stops, inclusive and finite.
+ /// The expected relative accuracy of the approximation.
/// Approximation of the finite integral in the given interval.
- public static double OnClosedInterval(Func f, double intervalBegin, double intervalEnd)
+ public static double OnOpenInterval(Func f, double intervalBegin, double intervalEnd, double targetAbsoluteError = 1E-8)
{
- return DoubleExponentialTransformation.Integrate(f, intervalBegin, intervalEnd, 1e-8);
+ // Reference:
+ // 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 (intervalBegin > intervalEnd)
+ {
+ return -OnOpenInterval(f, intervalEnd, intervalBegin, targetAbsoluteError);
+ }
+
+ // (-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 (double.IsInfinity(intervalBegin) && double.IsInfinity(intervalEnd))
+ {
+ Func u = (t) =>
+ {
+ return f(t / (1 - t * t)) * (1 + t * t) / ((1 - t * t) * (1 - t * t));
+ };
+ return OnClosedInterval(u, -1d, 1d, targetAbsoluteError);
+ }
+ // [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 (double.IsInfinity(intervalEnd))
+ {
+ Func u = (s) =>
+ {
+ return 2 * s * f(intervalBegin + (s / (1 - s)) * (s / (1 - s))) / ((1 - s) * (1 - s) * (1 - s));
+ };
+ return OnClosedInterval(u, 0d, 1d, targetAbsoluteError);
+ }
+ // (-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 (double.IsInfinity(intervalBegin))
+ {
+ Func u = (s) =>
+ {
+ return -2 * s * f(intervalEnd - s / (1 + s) * (s / (1 + s))) / ((1 + s) * (1 + s) * (1 + s));
+ };
+ return OnClosedInterval(u, -1d, 0d, targetAbsoluteError);
+ }
+ // [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 u = (t) =>
+ {
+ return f((intervalEnd - intervalBegin) / 4 * t * (3 - t * t) + (intervalEnd + intervalBegin) / 2) * 3 * (intervalEnd - intervalBegin) / 4 * (1 - t * t);
+ };
+ return OnClosedInterval(u, -1d, 1d, targetAbsoluteError);
+ }
}
///