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

Add Riemann–Liouville fractional derivative.

v4
diluculo 7 years ago
parent
510a82b8e0
commit
ae690ed557
2 changed files with 209 additions and 0 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. 75
      src/Numerics.Tests/FractionalCalculusTests/RiemannLiouvilleTests.cs
  2. 134
      src/Numerics/DifferIntegrate.cs

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;
}
}
}
}
Write Preview
Loading…
Cancel
Save