ufkc7e3d
/
mathnet-numerics
forked from tsai/mathnet-numerics
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

Remove substitution rule for [a, b] => [-1, 1] from DoubleExponential method.

uap-experimental
diluculo 7 years ago
parent
ab1f608345
commit
e0cd83b648
3 changed files with 42 additions and 63 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Unified View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 15
      src/Numerics.Tests/IntegrationTests/IntegrationTest.cs
  2. 66
      src/Numerics/Integrate.cs
  3. 24
      src/Numerics/Integration/GaussKronrodRule.cs

15
src/Numerics.Tests/IntegrationTests/IntegrationTest.cs
View File

@ -181,16 +181,13 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests
TargetAreaC, TargetAreaC,
Integrate.DoubleExponential(TargetFunctionC, StartC, StopC, 1e-10), Integrate.DoubleExponential(TargetFunctionC, StartC, StopC, 1e-10),
1e-10, 1e-10,
"DoubleExponential, Target 1e-10"); "DoubleExponential");
// integrate_(0)^(1) log(x) dx = -1
// Note that DoubleExponential returns -oo.
Assert.AreEqual( Assert.AreEqual(
TargetAreaD, TargetAreaD,
Integrate.OnClosedInterval(TargetFunctionD, StartD, StopD), Integrate.DoubleExponential(TargetFunctionD, StartD, StopD),
1e-10, 1e-10,
"Interval"); "DoubleExponential");
Assert.AreEqual( Assert.AreEqual(
TargetAreaD, TargetAreaD,
@ -474,19 +471,19 @@ namespace MathNet.Numerics.UnitTests.IntegrationTests
expected, expected,
factor * Integrate.DoubleExponential((x) => 1 / (1 + x * x), a, b), factor * Integrate.DoubleExponential((x) => 1 / (1 + x * x), a, b),
1e-10, 1e-10,
"DET Integral of sin(pi*x)/(pi*x) from -oo to oo"); "DET Integral of sin(pi*x)/(pi*x) from {0} to {1}", a, b);
Assert.AreEqual( Assert.AreEqual(
expected, expected,
factor * Integrate.GaussKronrod((x) => 1 / (1 + x * x), a, b), factor * Integrate.GaussKronrod((x) => 1 / (1 + x * x), a, b),
1e-10, 1e-10,
"GK Integral of sin(pi*x)/(pi*x) from -oo to oo"); "GK Integral of sin(pi*x)/(pi*x) from {0} to {1}", a, b);
Assert.AreEqual( Assert.AreEqual(
expected, expected,
factor * Integrate.GaussLegendre((x) => 1 / (1 + x * x), a, b, order: 128), factor * Integrate.GaussLegendre((x) => 1 / (1 + x * x), a, b, order: 128),
1e-10, 1e-10,
"GL Integral of sin(pi*x)/(pi*x) from -oo to oo"); "GL Integral of sin(pi*x)/(pi*x) from {0} to {1}", a, b);
} }
// integral_(-oo)^(oo) 1/(1 + j x^2) dx = -(-1)^(3/4) ¥ð // integral_(-oo)^(oo) 1/(1 + j x^2) dx = -(-1)^(3/4) ¥ð

66
src/Numerics/Integrate.cs
View File

@ -102,7 +102,7 @@ namespace MathNet.Numerics
// (-oo, oo) => [-1, 1] // (-oo, oo) => [-1, 1]
// //
// integral_{-oo}^{oo} f(x) dx = integral_{-1}^{1} f(g(t)) g'(t) dt // integral_(-oo)^(oo) f(x) dx = integral_(-1)^(1) f(g(t)) g'(t) dt
// g(t) = t / (1 - t^2) // g(t) = t / (1 - t^2)
// g'(t) = (1 + t^2) / (1 - t^2)^2 // g'(t) = (1 + t^2) / (1 - t^2)^2
if (double.IsInfinity(intervalBegin) && double.IsInfinity(intervalEnd)) if (double.IsInfinity(intervalBegin) && double.IsInfinity(intervalEnd))
@ -115,8 +115,8 @@ namespace MathNet.Numerics
} }
// [a, oo) => [0, 1] // [a, oo) => [0, 1]
// //
// integral_{a}^{oo} f(x) dx = integral_{0}^{oo} f(a + t^2) 2 t dt // 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 // = integral_(0)^(1) f(a + g(s)^2) 2 g(s) g'(s) ds
// g(s) = s / (1 - s) // g(s) = s / (1 - s)
// g'(s) = 1 / (1 - s)^2 // g'(s) = 1 / (1 - s)^2
else if (double.IsInfinity(intervalEnd)) else if (double.IsInfinity(intervalEnd))
@ -129,8 +129,8 @@ namespace MathNet.Numerics
} }
// (-oo, b] => [-1, 0] // (-oo, b] => [-1, 0]
// //
// integral_{-oo}^{b} f(x) dx = -integral_{-oo}^{0} f(b - t^2) 2 t dt // 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 // = -integral_(-1)^(0) f(b - g(s)^2) 2 g(s) g'(s) ds
// g(s) = s / (1 + s) // g(s) = s / (1 + s)
// g'(s) = 1 / (1 + s)^2 // g'(s) = 1 / (1 + s)^2
else if (double.IsInfinity(intervalBegin)) else if (double.IsInfinity(intervalBegin))
@ -141,18 +141,9 @@ namespace MathNet.Numerics
}; };
return DoubleExponentialTransformation.Integrate(u, -1, 0, targetAbsoluteError); return DoubleExponentialTransformation.Integrate(u, -1, 0, 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 else
{ {
Func<double, double> u = (t) => return DoubleExponentialTransformation.Integrate(f, intervalBegin, intervalEnd, targetAbsoluteError);
{
return f((intervalEnd - intervalBegin) / 4 * t * (3 - t * t) + (intervalEnd + intervalBegin) / 2) * 3 * (intervalEnd - intervalBegin) / 4 * (1 - t * t);
};
return DoubleExponentialTransformation.Integrate(u, -1, 1, targetAbsoluteError);
} }
} }
@ -178,7 +169,7 @@ namespace MathNet.Numerics
// (-oo, oo) => [-1, 1] // (-oo, oo) => [-1, 1]
// //
// integral_{-oo}^{oo} f(x) dx = integral_{-1}^{1} f(g(t)) g'(t) dt // integral_(-oo)^(oo) f(x) dx = integral_(-1)^(1) f(g(t)) g'(t) dt
// g(t) = t / (1 - t^2) // g(t) = t / (1 - t^2)
// g'(t) = (1 + t^2) / (1 - t^2)^2 // g'(t) = (1 + t^2) / (1 - t^2)^2
if (double.IsInfinity(intervalBegin) && double.IsInfinity(intervalEnd)) if (double.IsInfinity(intervalBegin) && double.IsInfinity(intervalEnd))
@ -191,8 +182,8 @@ namespace MathNet.Numerics
} }
// [a, oo) => [0, 1] // [a, oo) => [0, 1]
// //
// integral_{a}^{oo} f(x) dx = integral_{0}^{oo} f(a + t^2) 2 t dt // 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 // = integral_(0)^(1) f(a + g(s)^2) 2 g(s) g'(s) ds
// g(s) = s / (1 - s) // g(s) = s / (1 - s)
// g'(s) = 1 / (1 - s)^2 // g'(s) = 1 / (1 - s)^2
else if (double.IsInfinity(intervalEnd)) else if (double.IsInfinity(intervalEnd))
@ -205,8 +196,8 @@ namespace MathNet.Numerics
} }
// (-oo, b] => [-1, 0] // (-oo, b] => [-1, 0]
// //
// integral_{-oo}^{b} f(x) dx = -integral_{-oo}^{0} f(b - t^2) 2 t dt // 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 // = -integral_(-1)^(0) f(b - g(s)^2) 2 g(s) g'(s) ds
// g(s) = s / (1 + s) // g(s) = s / (1 + s)
// g'(s) = 1 / (1 + s)^2 // g'(s) = 1 / (1 + s)^2
else if (double.IsInfinity(intervalBegin)) else if (double.IsInfinity(intervalBegin))
@ -219,7 +210,7 @@ namespace MathNet.Numerics
} }
// [a, b] => [-1, 1] // [a, b] => [-1, 1]
// //
// integral_{a}^{b} f(x) dx = integral_{-1}^{1} f(g(t)) g'(t) dt // 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) = (b - a) * t * (3 - t^2) / 4 + (b + a) / 2
// g'(t) = 3 / 4 * (b - a) * (1 - t^2) // g'(t) = 3 / 4 * (b - a) * (1 - t^2)
else else
@ -292,7 +283,7 @@ namespace MathNet.Numerics
// (-oo, oo) => [-1, 1] // (-oo, oo) => [-1, 1]
// //
// integral_{-oo}^{oo} f(x) dx = integral_{-1}^{1} f(g(t)) g'(t) dt // integral_(-oo)^(oo) f(x) dx = integral_(-1)^(1) f(g(t)) g'(t) dt
// g(t) = t / (1 - t^2) // g(t) = t / (1 - t^2)
// g'(t) = (1 + t^2) / (1 - t^2)^2 // g'(t) = (1 + t^2) / (1 - t^2)^2
if (double.IsInfinity(intervalBegin) && double.IsInfinity(intervalEnd)) if (double.IsInfinity(intervalBegin) && double.IsInfinity(intervalEnd))
@ -305,8 +296,8 @@ namespace MathNet.Numerics
} }
// [a, oo) => [0, 1] // [a, oo) => [0, 1]
// //
// integral_{a}^{oo} f(x) dx = integral_{0}^{oo} f(a + t^2) 2 t dt // 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 // = integral_(0)^(1) f(a + g(s)^2) 2 g(s) g'(s) ds
// g(s) = s / (1 - s) // g(s) = s / (1 - s)
// g'(s) = 1 / (1 - s)^2 // g'(s) = 1 / (1 - s)^2
else if (double.IsInfinity(intervalEnd)) else if (double.IsInfinity(intervalEnd))
@ -319,8 +310,8 @@ namespace MathNet.Numerics
} }
// (-oo, b] => [-1, 0] // (-oo, b] => [-1, 0]
// //
// integral_{-oo}^{b} f(x) dx = -integral_{-oo}^{0} f(b - t^2) 2 t dt // 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 // = -integral_(-1)^(0) f(b - g(s)^2) 2 g(s) g'(s) ds
// g(s) = s / (1 + s) // g(s) = s / (1 + s)
// g'(s) = 1 / (1 + s)^2 // g'(s) = 1 / (1 + s)^2
else if (double.IsInfinity(intervalBegin)) else if (double.IsInfinity(intervalBegin))
@ -331,18 +322,9 @@ namespace MathNet.Numerics
}; };
return DoubleExponentialTransformation.ContourIntegrate(u, -1, 0, targetAbsoluteError); return DoubleExponentialTransformation.ContourIntegrate(u, -1, 0, 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 else
{ {
Func<double, Complex> u = (t) => return DoubleExponentialTransformation.ContourIntegrate(f, intervalBegin, intervalEnd, targetAbsoluteError);
{
return f((intervalEnd - intervalBegin) / 4 * t * (3 - t * t) + (intervalEnd + intervalBegin) / 2) * 3 * (intervalEnd - intervalBegin) / 4 * (1 - t * t);
};
return DoubleExponentialTransformation.ContourIntegrate(u, -1, 1, targetAbsoluteError);
} }
} }
@ -368,7 +350,7 @@ namespace MathNet.Numerics
// (-oo, oo) => [-1, 1] // (-oo, oo) => [-1, 1]
// //
// integral_{-oo}^{oo} f(x) dx = integral_{-1}^{1} f(g(t)) g'(t) dt // integral_(-oo)^(oo) f(x) dx = integral_(-1)^(1) f(g(t)) g'(t) dt
// g(t) = t / (1 - t^2) // g(t) = t / (1 - t^2)
// g'(t) = (1 + t^2) / (1 - t^2)^2 // g'(t) = (1 + t^2) / (1 - t^2)^2
if (double.IsInfinity(intervalBegin) && double.IsInfinity(intervalEnd)) if (double.IsInfinity(intervalBegin) && double.IsInfinity(intervalEnd))
@ -381,8 +363,8 @@ namespace MathNet.Numerics
} }
// [a, oo) => [0, 1] // [a, oo) => [0, 1]
// //
// integral_{a}^{oo} f(x) dx = integral_{0}^{oo} f(a + t^2) 2 t dt // 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 // = integral_(0)^(1) f(a + g(s)^2) 2 g(s) g'(s) ds
// g(s) = s / (1 - s) // g(s) = s / (1 - s)
// g'(s) = 1 / (1 - s)^2 // g'(s) = 1 / (1 - s)^2
else if (double.IsInfinity(intervalEnd)) else if (double.IsInfinity(intervalEnd))
@ -395,8 +377,8 @@ namespace MathNet.Numerics
} }
// (-oo, b] => [-1, 0] // (-oo, b] => [-1, 0]
// //
// integral_{-oo}^{b} f(x) dx = -integral_{-oo}^{0} f(b - t^2) 2 t dt // 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 // = -integral_(-1)^(0) f(b - g(s)^2) 2 g(s) g'(s) ds
// g(s) = s / (1 + s) // g(s) = s / (1 + s)
// g'(s) = 1 / (1 + s)^2 // g'(s) = 1 / (1 + s)^2
else if (double.IsInfinity(intervalBegin)) else if (double.IsInfinity(intervalBegin))
@ -409,7 +391,7 @@ namespace MathNet.Numerics
} }
// [a, b] => [-1, 1] // [a, b] => [-1, 1]
// //
// integral_{a}^{b} f(x) dx = integral_{-1}^{1} f(g(t)) g'(t) dt // 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) = (b - a) * t * (3 - t^2) / 4 + (b + a) / 2
// g'(t) = 3 / 4 * (b - a) * (1 - t^2) // g'(t) = 3 / 4 * (b - a) * (1 - t^2)
else else

24
src/Numerics/Integration/GaussKronrodRule.cs
View File

@ -240,7 +240,7 @@ namespace MathNet.Numerics.Integration
// (-oo, oo) => [-1, 1] // (-oo, oo) => [-1, 1]
// //
// integral_{-oo}^{oo} f(x) dx = integral_{-1}^{1} f(g(t)) g'(t) dt // integral_(-oo)^(oo) f(x) dx = integral_(-1)^(1) f(g(t)) g'(t) dt
// g(t) = t / (1 - t^2) // g(t) = t / (1 - t^2)
// g'(t) = (1 + t^2) / (1 - t^2)^2 // g'(t) = (1 + t^2) / (1 - t^2)^2
if ((intervalBegin < double.MinValue) && (intervalEnd > double.MaxValue)) if ((intervalBegin < double.MinValue) && (intervalEnd > double.MaxValue))
@ -253,8 +253,8 @@ namespace MathNet.Numerics.Integration
} }
// [a, oo) => [0, 1] // [a, oo) => [0, 1]
// //
// integral_{a}^{oo} f(x) dx = integral_{0}^{oo} f(a + t^2) 2 t dt // 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 // = integral_(0)^(1) f(a + g(s)^2) 2 g(s) g'(s) ds
// g(s) = s / (1 - s) // g(s) = s / (1 - s)
// g'(s) = 1 / (1 - s)^2 // g'(s) = 1 / (1 - s)^2
else if (intervalEnd > double.MaxValue) else if (intervalEnd > double.MaxValue)
@ -267,8 +267,8 @@ namespace MathNet.Numerics.Integration
} }
// (-oo, b] => [-1, 0] // (-oo, b] => [-1, 0]
// //
// integral_{-oo}^{b} f(x) dx = -integral_{-oo}^{0} f(b - t^2) 2 t dt // 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 // = -integral_(-1)^(0) f(b - g(s)^2) 2 g(s) g'(s) ds
// g(s) = s / (1 + s) // g(s) = s / (1 + s)
// g'(s) = 1 / (1 + s)^2 // g'(s) = 1 / (1 + s)^2
else if (intervalBegin < double.MinValue) else if (intervalBegin < double.MinValue)
@ -281,7 +281,7 @@ namespace MathNet.Numerics.Integration
} }
// [a, b] => [-1, 1] // [a, b] => [-1, 1]
// //
// integral_{a}^{b} f(x) dx = integral_{-1}^{1} f(g(t)) g'(t) dt // 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) = (b - a) * t * (3 - t^2) / 4 + (b + a) / 2
// g'(t) = 3 / 4 * (b - a) * (1 - t^2) // g'(t) = 3 / 4 * (b - a) * (1 - t^2)
else else
@ -326,7 +326,7 @@ namespace MathNet.Numerics.Integration
// (-oo, oo) => [-1, 1] // (-oo, oo) => [-1, 1]
// //
// integral_{-oo}^{oo} f(x) dx = integral_{-1}^{1} f(g(t)) g'(t) dt // integral_(-oo)^(oo) f(x) dx = integral_(-1)^(1) f(g(t)) g'(t) dt
// g(t) = t / (1 - t^2) // g(t) = t / (1 - t^2)
// g'(t) = (1 + t^2) / (1 - t^2)^2 // g'(t) = (1 + t^2) / (1 - t^2)^2
if ((intervalBegin < double.MinValue) && (intervalEnd > double.MaxValue)) if ((intervalBegin < double.MinValue) && (intervalEnd > double.MaxValue))
@ -339,8 +339,8 @@ namespace MathNet.Numerics.Integration
} }
// [a, oo) => [0, 1] // [a, oo) => [0, 1]
// //
// integral_{a}^{oo} f(x) dx = integral_{0}^{oo} f(a + t^2) 2 t dt // 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 // = integral_(0)^(1) f(a + g(s)^2) 2 g(s) g'(s) ds
// g(s) = s / (1 - s) // g(s) = s / (1 - s)
// g'(s) = 1 / (1 - s)^2 // g'(s) = 1 / (1 - s)^2
else if (intervalEnd > double.MaxValue) else if (intervalEnd > double.MaxValue)
@ -353,8 +353,8 @@ namespace MathNet.Numerics.Integration
} }
// (-oo, b] => [-1, 0] // (-oo, b] => [-1, 0]
// //
// integral_{-oo}^{b} f(x) dx = -integral_{-oo}^{0} f(b - t^2) 2 t dt // 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 // = -integral_(-1)^(0) f(b - g(s)^2) 2 g(s) g'(s) ds
// g(s) = s / (1 + s) // g(s) = s / (1 + s)
// g'(s) = 1 / (1 + s)^2 // g'(s) = 1 / (1 + s)^2
else if (intervalBegin < double.MinValue) else if (intervalBegin < double.MinValue)
@ -367,7 +367,7 @@ namespace MathNet.Numerics.Integration
} }
// [a, b] => [-1, 1] // [a, b] => [-1, 1]
// //
// integral_{a}^{b} f(x) dx = integral_{-1}^{1} f(g(t)) g'(t) dt // 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) = (b - a) * t * (3 - t^2) / 4 + (b + a) / 2
// g'(t) = 3 / 4 * (b - a) * (1 - t^2) // g'(t) = 3 / 4 * (b - a) * (1 - t^2)
else else

Write Preview
Loading…
Cancel
Save