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

Merge pull request #662 from diluculo/fractionalCalculus 

Fractional calculus
v4
Christoph Ruegg 6 years ago
committed by GitHub
parent
b5fa485b0b ae690ed557
commit
0686faa293
No known key found for this signature in database GPG Key ID: 4AEE18F83AFDEB23
4 changed files with 211 additions and 2 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 2
      RELEASENOTES.md
  2. 75
      src/Numerics.Tests/FractionalCalculusTests/RiemannLiouvilleTests.cs
  3. 134
      src/Numerics/DifferIntegrate.cs
  4. 2
      src/Numerics/Integrate.cs

2
RELEASENOTES.md
View File

@ -7,7 +7,7 @@
### 4.9.0 - 2019-10-13
* Integration: Substitution for open intervals *~Jong Hyun Kim*
* Integration: Gaus-Kronrod Quadrature *~Jong Hyun Kim*
* Integration: Gauss-Kronrod Quadrature *~Jong Hyun Kim*
* Integration: Complex contour *~Jong Hyun Kim*
* Distributions: Three parameter Burr (Type XII) *~Philip Gruber*
* Distributions: Inverse Gaussian *~Philip Gruber*

75
src/Numerics.Tests/FractionalCalculusTests/RiemannLiouvilleTests.cs
View File

@ -0,0 +1,75 @@
using NUnit.Framework;
using System;
namespace MathNet.Numerics.Tests.FractionalCalculus
{
[TestFixture, Category("DifferIntegration")]
public class RiemannLiouvilleTests
{
[TestCase(3.0, 3.0, 0.0)] // D(3) x^2 = 0
[TestCase(2.5, 3.0, 0.651470015870559895448512830451)] // D(5/2) x^2 = Gamma(3)/Gamma(1/2)*x^(-1/2)
[TestCase(2.0, 3.0, 2.0)] // D(2) x^2 = 2
[TestCase(1.5, 3.0, 3.90882009522335937269107698271)] // D(3/2) x^2 = Gamma(3)/Gamma(3/2)*x^(1/2)
[TestCase(1.0, 3.0, 6.0)] // D(1) x^2 = 2 x
[TestCase(0.5, 3.0, 7.81764019044671874538215396541)] // D(1/2) x^2 = Gamma(3)/Gamma(5/2)*x^(3/2)
[TestCase(0.0, 3.0, 9.0)] // D(0) x^2 = x^2
[TestCase(-0.5, 3.0, 9.38116822853606249445858475850)] // D(-1/2) x^2 = Gamma(3)/Gamma(7/2)*x^(5/2)
[TestCase(-1.0, 3.0, 9.0)] // D(-1) x^2 = x^3/3
[TestCase(-1.5, 3.0, 8.04100133874519642382164407871)] // D(-3/2) x^2 = Gamma(3)/Gamma(9/2)*x^(7/2)
[TestCase(-2.0, 3.0, 6.75)] // D(-2) x^2 = x^4/12
[TestCase(-2.5, 3.0, 5.36066755916346428254776271914)] // D(-5/2) x^2 = Gamma(3)/Gamma(11/2)*x^(9/2)
[TestCase(-3.0, 3.0, 4.05)] // D(-3) x^2 = x^5/60
public void TestPolynomial(double n, double x, double expected)
{
Func<double, double> f = (t) => Math.Pow(t, 2);
Assert.LessOrEqual(
(expected - DifferIntegrate.DoubleExponential(f, x, n, x0: 0.0)) / expected,
1E-7,
"D({0}) x^2 (Double-Exponetial) where x = {1}", n, x);
Assert.LessOrEqual(
(expected - DifferIntegrate.GaussLegendre(f, x, n, x0: 0.0)) / expected,
1E-6,
"D({0}) x^2 (Gauss-Legendre) where x = {1}", n, x);
Assert.LessOrEqual(
(expected - DifferIntegrate.GaussKronrod(f, x, n, x0: 0.0)) / expected,
1E-4,
"D({0}) x^2 (Gauss-Kronrod) where x = {1}", n, x);
}
[TestCase(3.0, -1.5, 0.278538406900792139127369841989)] // D(3) exp(π x) = exp(π x) π^3,
[TestCase(2.5, -1.5, 0.157148467791413402880359987678)] // D(5/2) exp(π x) = exp(π x) π^(2.5),
[TestCase(2.0, -1.5, 0.0886615285984055199686857617477)] // D(2) exp(π x) = exp(π x) π^2,
[TestCase(1.5, -1.5, 0.0500219108966418944762345595265)] // D(3/2) exp(π x) = exp(π x) π^(1.5),
[TestCase(1.0, -1.5, 0.0282218410770393627236690193887)] // D(1) exp(π x) = exp(π x) π,
[TestCase(0.5, -1.5, 0.0159224687642057998088504370906)] // D(1/2) exp(π x) = exp(π x) π^(0.5),
[TestCase(0.0, -1.5, 0.00898329102112942788966495190793)] // D(0) exp(π x) = exp(π x)
[TestCase(-0.5, 3.5, 33631.1952282529497668583487591)] // D(-1/2) exp(π x) = exp(π x)/π^(0.5)
[TestCase(-1.0, 3.5, 18974.3700300413201713375621666)] // D(-1) exp(π x) = exp(π x)/π
[TestCase(-1.5, 3.5, 10705.1419253300403750707800464)] // D(-3/2) exp(π x) = exp(π x)/π^(1.5)
[TestCase(-2.0, 3.5, 6039.72956467158140885534431539)] // D(-2) exp(π x) = exp(π x)/π^2
[TestCase(-2.5, 3.5, 3407.55250783313088752769495213)] // D(-5/2) exp(π x) = exp(π x)/π^(2.5)
[TestCase(-3.0, 3.5, 1922.50563031148665828996231149)] // D(-3) exp(π x) = exp(π x)/π^3
public void TestExponential(double n, double x, double expected)
{
Func<double, double> f = (t) => Math.Exp(Math.PI * t);
Assert.LessOrEqual(
(expected - DifferIntegrate.DoubleExponential(f, x, n, x0: double.NegativeInfinity)) / expected,
1E-9,
"D({0}) exp(π x) (Double-Exponential) where x = {1}", n, x);
Assert.LessOrEqual(
(expected - DifferIntegrate.GaussLegendre(f, x, n, x0: double.NegativeInfinity)) / expected,
1E-6,
"D({0}) exp(π x) (Gauss-Legendre) where x = {1}", n, x);
Assert.LessOrEqual(
(expected - DifferIntegrate.GaussKronrod(f, x, n, x0: double.NegativeInfinity)) / expected,
1E-7,
"D({0}) exp(π x) (Gauss-Kronrod) where x = {1}", n, x);
}
}
}

134
src/Numerics/DifferIntegrate.cs
View File

@ -0,0 +1,134 @@
using System;
namespace MathNet.Numerics
{
public static class DifferIntegrate
{
/// <summary>
/// Evaluates the Riemann-Liouville fractional derivative that uses the double exponential integration.
/// </summary>
/// <remarks>
/// <para>order = 1.0 : normal derivative</para>
/// <para>order = 0.5 : semi-derivative</para>
/// <para>order = -0.5 : semi-integral</para>
/// <para>order = -1.0 : normal integral</para>
/// </remarks>
/// <param name="f">The analytic smooth function to differintegrate.</param>
/// <param name="x">The evaluation point.</param>
/// <param name="order">The order of fractional derivative.</param>
/// <param name="x0">The reference point of integration.</param>
/// <param name="targetAbsoluteError">The expected relative accuracy of the Double-Exponential integration.</param>
/// <returns>Approximation of the differintegral of order n at x.</returns>
public static double DoubleExponential(Func<double, double> f, double x, double order, double x0 = 0, double targetAbsoluteError = 1E-10)
{
// The Riemann–Liouville fractional derivative of f(x) of order n is defined as
// \,_{x_0}{\mathbb{D}}^n_xf(x) = \frac{1}{\Gamma(m-n)} \frac{d^m}{dx^m} \int_{x_0}^{x} (x-t)^{m-n-1} f(t) dt
// where m is the smallest interger greater than n.
// see https://en.wikipedia.org/wiki/Differintegral
if (Math.Abs(order) < double.Epsilon)
{
return f(x);
}
else if (order > 0 && Math.Abs(order - (int)order) < double.Epsilon)
{
return Differentiate.Derivative(f, x, (int)order);
}
else
{
int m = (int)Math.Ceiling(order) + 1;
if (m < 1) m = 1;
double r = m - order - 1;
Func<double, double> g = (v) => Integrate.DoubleExponential((t) => Math.Pow(v - t, r) * f(t), x0, v, targetAbsoluteError: targetAbsoluteError);
double numerator = Differentiate.Derivative(g, x, m);
double denominator = SpecialFunctions.Gamma(m - order);
return numerator / denominator;
}
}
/// <summary>
/// Evaluates the Riemann-Liouville fractional derivative that uses the Gauss-Legendre integration.
/// </summary>
/// <remarks>
/// <para>order = 1.0 : normal derivative</para>
/// <para>order = 0.5 : semi-derivative</para>
/// <para>order = -0.5 : semi-integral</para>
/// <para>order = -1.0 : normal integral</para>
/// </remarks>
/// <param name="f">The analytic smooth function to differintegrate.</param>
/// <param name="x">The evaluation point.</param>
/// <param name="order">The order of fractional derivative.</param>
/// <param name="x0">The reference point of integration.</param>
/// <param name="gaussLegendrePoints">The number of Gauss-Legendre points.</param>
/// <returns>Approximation of the differintegral of order n at x.</returns>
public static double GaussLegendre(Func<double, double> f, double x, double order, double x0 = 0, int gaussLegendrePoints = 128)
{
// The Riemann–Liouville fractional derivative of f(x) of order n is defined as
// \,_{x_0}{\mathbb{D}}^n_xf(x) = \frac{1}{\Gamma(m-n)} \frac{d^m}{dx^m} \int_{x_0}^{x} (x-t)^{m-n-1} f(t) dt
// where m is the smallest interger greater than n.
// see https://en.wikipedia.org/wiki/Differintegral
if (Math.Abs(order) < double.Epsilon)
{
return f(x);
}
else if (order > 0 && Math.Abs(order - (int)order) < double.Epsilon)
{
return Differentiate.Derivative(f, x, (int)order);
}
else
{
int m = (int)Math.Ceiling(order) + 1;
if (m < 1) m = 1;
double r = m - order - 1;
Func<double, double> g = (v) => Integrate.GaussLegendre((t) => Math.Pow(v - t, r) * f(t), x0, v, order: gaussLegendrePoints);
double numerator = Differentiate.Derivative(g, x, m);
double denominator = SpecialFunctions.Gamma(m - order);
return numerator / denominator;
}
}
/// <summary>
/// Evaluates the Riemann-Liouville fractional derivative that uses the Gauss-Kronrod integration.
/// </summary>
/// <remarks>
/// <para>order = 1.0 : normal derivative</para>
/// <para>order = 0.5 : semi-derivative</para>
/// <para>order = -0.5 : semi-integral</para>
/// <para>order = -1.0 : normal integral</para>
/// </remarks>
/// <param name="f">The analytic smooth function to differintegrate.</param>
/// <param name="x">The evaluation point.</param>
/// <param name="order">The order of fractional derivative.</param>
/// <param name="x0">The reference point of integration.</param>
/// <param name="targetRelativeError">The expected relative accuracy of the Gauss-Kronrod integration.</param>
/// <param name="gaussKronrodPoints">The number of Gauss-Kronrod points. Pre-computed for 15, 21, 31, 41, 51 and 61 points.</param>
/// <returns>Approximation of the differintegral of order n at x.</returns>
public static double GaussKronrod(Func<double, double> f, double x, double order, double x0 = 0, double targetRelativeError = 1E-10, int gaussKronrodPoints = 15)
{
// The Riemann–Liouville fractional derivative of f(x) of order n is defined as
// \,_{x_0}{\mathbb{D}}^n_xf(x) = \frac{1}{\Gamma(m-n)} \frac{d^m}{dx^m} \int_{x_0}^{x} (x-t)^{m-n-1} f(t) dt
// where m is the smallest interger greater than n.
// see https://en.wikipedia.org/wiki/Differintegral
if (Math.Abs(order) < double.Epsilon)
{
return f(x);
}
else if (order > 0 && Math.Abs(order - (int)order) < double.Epsilon)
{
return Differentiate.Derivative(f, x, (int)order);
}
else
{
int m = (int)Math.Ceiling(order) + 1;
if (m < 1) m = 1;
double r = m - order - 1;
Func<double, double> g = (v) => Integrate.GaussKronrod((t) => Math.Pow(v - t, r) * f(t), x0, v, targetRelativeError: targetRelativeError, order: gaussKronrodPoints);
double numerator = Differentiate.Derivative(g, x, m);
double denominator = SpecialFunctions.Gamma(m - order);
return numerator / denominator;
}
}
}
}

2
src/Numerics/Integrate.cs
View File

@ -243,7 +243,7 @@ namespace MathNet.Numerics
/// <param name="intervalEnd">Where the interval stops.</param>
/// <param name="targetRelativeError">The expected relative accuracy of the approximation.</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, 31, 41, 51 and 61 points.</param>
/// <param name="order">The number of Gauss-Kronrod points. Pre-computed for 15, 21, 31, 41, 51 and 61 points.</param>
/// <returns>Approximation of the finite integral in the given interval.</returns>
public static double GaussKronrod(Func<double, double> f, double intervalBegin, double intervalEnd, double targetRelativeError = 1E-8, int maximumDepth = 15, int order = 15)
{

Write Preview
Loading…
Cancel
Save