diff --git a/src/Numerics/Differentiation/FiniteDifferenceCoefficients.cs b/src/Numerics/Differentiation/FiniteDifferenceCoefficients.cs
new file mode 100644
index 00000000..6c5de78f
--- /dev/null
+++ b/src/Numerics/Differentiation/FiniteDifferenceCoefficients.cs
@@ -0,0 +1,140 @@
+//
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2015 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+//
+
+using System;
+using MathNet.Numerics.LinearAlgebra.Double;
+
+namespace MathNet.Numerics.Differentiation
+{
+ ///
+ /// Class to calculate finite difference coefficients using Taylor series expansion method.
+ ///
+ ///
+ /// For n points, coefficients are calculated up to the maximum derivative order possible (n-1).
+ /// The current function value position specifies the "center" for surrounding coefficients.
+ /// Selecting the first, middle or last positions represent forward, backwards and central difference methods.
+ ///
+ ///
+ ///
+ public class FiniteDifferenceCoefficients
+ {
+ ///
+ /// Number of points for finite difference coefficients. Changing this value recalculates the coefficients table.
+ ///
+ public int Points
+ {
+ get { return _points; }
+ set
+ {
+ CalculateCoefficients(value);
+ _points = value;
+ }
+ }
+
+ private double[][,] _coefficients;
+ private int _points;
+
+ ///
+ /// Initializes a new instance of the class.
+ ///
+ /// Number of finite difference coefficients.
+ public FiniteDifferenceCoefficients(int points)
+ {
+ Points = points;
+ CalculateCoefficients(Points);
+ }
+
+ ///
+ /// Gets the finite difference coefficients for a specified center and order.
+ ///
+ /// Current function position with respect to coefficients. Must be within point range.
+ /// Order of finite difference coefficients.
+ /// Vector of finite difference coefficients.
+ public double[] GetCoefficients(int center, int order)
+ {
+ if (center >= _coefficients.Length)
+ throw new ArgumentOutOfRangeException("center", "Center position must be within the point range.");
+ if (order >= _coefficients.Length)
+ throw new ArgumentOutOfRangeException("order", "Maximum difference order is points-1.");
+
+ // Return proper row
+ var columns = _coefficients[center].GetLength(1);
+ var array = new double[columns];
+ for (int i = 0; i < columns; ++i)
+ array[i] = _coefficients[center][order, i];
+ return array;
+ }
+
+ ///
+ /// Gets the finite difference coefficients for all orders at a specified center.
+ ///
+ /// Current function position with respect to coefficients. Must be within point range.
+ /// Rectangular array of coefficients, with columns specifing order.
+ public double[,] GetCoefficientsForAllOrders(int center)
+ {
+ if (center >= _coefficients.Length)
+ throw new ArgumentOutOfRangeException("center", "Center position must be within the point range.");
+
+ return _coefficients[center];
+ }
+
+ private void CalculateCoefficients(int points)
+ {
+ var c = new double[points][,];
+
+ // For ever possible center given the number of points, compute ever possible coefficeint for all possible orders.
+ for (int center = 0; center < points; center++)
+ {
+ // Deltas matrix for center located at 'center'.
+ var A = new DenseMatrix(points);
+ var l = points - center - 1;
+ for (int row = points - 1; row >= 0; row--)
+ {
+ A[row, 0] = 1.0;
+ for (int col = 1; col < points; col++)
+ {
+ A[row, col] = A[row, col - 1] * l / col;
+ }
+ l -= 1;
+ }
+
+ c[center] = A.Inverse().ToArray();
+
+ // "Polish" results by rounding.
+ var fac = SpecialFunctions.Factorial(points);
+ for (int j = 0; j < points; j++)
+ for (int k = 0; k < points; k++)
+ c[center][j, k] = (Math.Round(c[center][j, k] * fac, MidpointRounding.AwayFromZero)) / fac;
+ }
+
+ _coefficients = c;
+ }
+ }
+}
diff --git a/src/Numerics/Differentiation/NumericalDerivative.cs b/src/Numerics/Differentiation/NumericalDerivative.cs
new file mode 100644
index 00000000..a0bfbf9a
--- /dev/null
+++ b/src/Numerics/Differentiation/NumericalDerivative.cs
@@ -0,0 +1,483 @@
+//
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2015 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+//
+
+using System;
+using System.Linq;
+
+
+namespace MathNet.Numerics.Differentiation
+{
+ ///
+ /// Type of finite different step size.
+ ///
+ public enum StepType
+ {
+ ///
+ /// The absolute step size value will be used in numerical derivatives, regardless of order or function parameters.
+ ///
+ Absolute,
+
+ ///
+ /// A base step size value, h, will be scaled according to the function input parameter. A common example is hx = h*(1+abs(x)), however
+ /// this may vary depending on implementation. This definition only guarantees that the only scaling will be relative to the
+ /// function input parameter and not the order of the finite difference derivative.
+ ///
+ RelativeX,
+
+ ///
+ /// A base step size value, eps (typically machine precision), is scaled according to the finite difference coefficient order
+ /// and function input parameter. The initial scaling according to finite different coefficient order can be thought of as producing a
+ /// base step size, h, that is equivalent to scaling. This stepsize is then scaled according to the function
+ /// input parameter. Although implementation may vary, an example of second order accurate scaling may be (eps)^(1/3)*(1+abs(x)).
+ ///
+ Relative
+ };
+
+ ///
+ /// Class to evaluate the numerical derivative of a function using finite difference approximations.
+ /// Variable point and center methods can be initialized .
+ /// This class can also be used to return function handles (delagates) for a fixed derivative order and variable.
+ /// It is possible to evaluate the derivative and partial derivative of univariate and multivariate functions respectively.
+ ///
+ public class NumericalDerivative
+ {
+ ///
+ /// Sets and gets the finite difference step size. This value is for each function evaluation if relative stepsize types are used.
+ /// If the base step size used in scaling is desired, see .
+ ///
+ ///
+ /// Setting then getting the StepSize may return a different value. This is not unusual since a user-defined step size is converted to a
+ /// base-2 representable number to improve finite difference accuracy.
+ ///
+ public double StepSize
+ {
+ get { return _stepSize; }
+ set
+ {
+ //Base 2 yields more accurate results...
+ var p = Math.Log(Math.Abs(value))/Math.Log(2);
+ _stepSize = Math.Pow(2, Math.Round(p));
+ }
+ }
+
+ ///
+ /// Sets and gets the base fininte difference step size. This assigned value to this parameter is only used if is set to RelativeX.
+ /// However, if the StepType is Relative, it will contain the base step size computed from based on the finite difference order.
+ ///
+ public double BaseStepSize
+ {
+ get { return _baseStepSize; }
+ set
+ {
+ //Base 2 yields more accurate results...
+ var p = Math.Log(Math.Abs(value)) / Math.Log(2);
+ _baseStepSize = Math.Pow(2, Math.Round(p));
+ }
+ }
+
+ ///
+ /// Sets and gets the base finite difference step size. This parameter is only used if is set to Relative.
+ /// By default this is set to machine epsilon, from which is computed.
+ ///
+ public double Epsilon
+ {
+ get { return _epsilon; }
+ set
+ {
+ //Base 2 yields more accurate results...
+ var p = Math.Log(Math.Abs(value)) / Math.Log(2);
+ _epsilon = Math.Pow(2, Math.Round(p));
+ }
+ }
+
+ ///
+ /// Sets and gets the location of the center point for the finite difference derivative.
+ ///
+ public int Center
+ {
+ get { return _center; }
+ set
+ {
+ if (value >= _points || value < 0)
+ throw new ArgumentOutOfRangeException("value", "Center must lie between 0 and points -1");
+ _center = value;
+ }
+ }
+
+ ///
+ /// Number of times a function is evaluated for numerical derivatives.
+ ///
+ public int Evaluations { get; private set; }
+
+ ///
+ /// Type of step size for computing finite differences. If set to absolute, dx = h.
+ /// If set to relative, dx = (1+abs(x))*h^(2/(order+1)). This provides accurate results when
+ /// h is approximately equal to the square-root of machine accuracy, epsilon.
+ ///
+ public StepType StepType
+ {
+ get { return _stepType; }
+ set { _stepType = value; }
+ }
+
+ private readonly int _points;
+ private int _center;
+ private double _stepSize = Math.Pow(2, -10);
+ private double _epsilon = Math.Pow(2, -52);
+ private double _baseStepSize = Math.Pow(2, -26);
+ private StepType _stepType = StepType.Relative;
+ private readonly FiniteDifferenceCoefficients _coefficients;
+
+ ///
+ /// Initializes a NumericalDerivative class with the default 3 point center difference method.
+ ///
+ public NumericalDerivative() : this(3, 1)
+ {
+ }
+
+ ///
+ /// Initialized a NumericalDerivative class.
+ ///
+ /// Number of points for finite difference derivatives.
+ /// Location of the center with respect to other points. Value ranges from zero to points-1.
+ public NumericalDerivative(int points, int center)
+ {
+ _center = center;
+ if (points < 2)
+ throw new ArgumentOutOfRangeException("points", "Points must be two or greater.");
+ _points = points;
+ Center = center;
+ _epsilon = CalculateMachineEpsilon();
+ _coefficients = new FiniteDifferenceCoefficients(points);
+ }
+
+ ///
+ /// Evaluates the derivative of equidistant points using the finite difference method.
+ ///
+ /// Vector of points StepSize apart.
+ /// Derivative order.
+ /// Finite difference step size.
+ /// Derivative of points of the specified order.
+ public double EvaluateDerivative(double[] points, int order, double stepSize)
+ {
+ if (points == null)
+ throw new ArgumentNullException("points");
+
+ if (order >= _points || order < 0)
+ throw new ArgumentOutOfRangeException("order", "Order must be between zero and points-1.");
+
+ var c = _coefficients.GetCoefficients(Center, order);
+ var result = c.Select((t, i) => t*points[i]).Sum();
+ result /= Math.Pow(stepSize, order);
+ return result;
+ }
+
+ ///
+ /// Evaluates the derivative of a scalar univariate function.
+ ///
+ ///
+ /// Supplying the optional argument currentValue will reduce the number of function evaluations
+ /// required to calculate the finite difference derivative.
+ ///
+ /// Function handle.
+ /// Point at which to compute the derivative.
+ /// Derivative order.
+ /// Current function value at center.
+ /// Function derivative at x of the specified order.
+ public double EvaluateDerivative(Func f, double x, int order, double? currentValue = null)
+ {
+ var c = _coefficients.GetCoefficients(Center, order);
+ var h = CalculateStepSize(_points, x, order);
+
+ var points = new double[_points];
+ for (int i = 0; i < _points; i++)
+ {
+ if (i == Center && currentValue.HasValue)
+ points[i] = currentValue.Value;
+ else if(c[i] != 0) // Only evaluate function if it will actually be used.
+ {
+ points[i] = f(x + (i - Center) * h);
+ Evaluations++;
+ }
+ }
+
+ return EvaluateDerivative(points, order, h);
+ }
+
+ ///
+ /// Creates a function handle for the derivative of a scalar univariate function.
+ ///
+ /// Input function handle.
+ /// Derivative order.
+ /// Function handle that evaluates the derivative of input function at a fixed order.
+ public Func CreateDerivativeFunctionHandle(Func f, int order)
+ {
+ return x => EvaluateDerivative(f, x, order);
+ }
+
+ ///
+ /// Evaluates the partial derivative of a multivariate function.
+ ///
+ /// Multivariate function handle.
+ /// Vector at which to evaluate the derivative.
+ /// Index of independent variable for partial derivative.
+ /// Derivative order.
+ /// Current function value at center.
+ /// Function partial derivative at x of the specified order.
+ public double EvaluatePartialDerivative(Func f, double[] x, int parameterIndex, int order, double? currentValue = null)
+ {
+ var xi = x[parameterIndex];
+ var c = _coefficients.GetCoefficients(Center, order);
+ var h = CalculateStepSize(_points, x[parameterIndex], order);
+
+ var points = new double[_points];
+ for (int i = 0; i < _points; i++)
+ {
+ if (i == Center && currentValue.HasValue)
+ points[i] = currentValue.Value;
+ else if(c[i] != 0) // Only evaluate function if it will actually be used.
+ {
+ x[parameterIndex] = xi + (i - Center) * h;
+ points[i] = f(x);
+ Evaluations++;
+ }
+ }
+
+ //restore original value
+ x[parameterIndex] = xi;
+ return EvaluateDerivative(points, order, h);
+ }
+
+ ///
+ /// Evaluates the partial derivatives of a multivariate function array.
+ ///
+ ///
+ /// This function assumes the input vector x is of the correct length for f.
+ ///
+ /// Multivariate vector function array handle.
+ /// Vector at which to evaluate the derivatives.
+ /// Index of independent variable for partial derivative.
+ /// Derivative order.
+ /// Current function value at center.
+ /// Vector of functions partial derivatives at x of the specified order.
+ public double[] EvaluatePartialDerivative(Func[] f, double[] x, int parameterIndex, int order, double?[] currentValue = null)
+ {
+ var df = new double[f.Length];
+ for (int i = 0; i < f.Length; i++)
+ {
+ if(currentValue != null && currentValue[i].HasValue)
+ df[i] = EvaluatePartialDerivative(f[i], x, parameterIndex, order, currentValue[i].Value);
+ else
+ df[i] = EvaluatePartialDerivative(f[i], x, parameterIndex, order);
+ }
+
+ return df;
+ }
+
+ ///
+ /// Creates a function handle for the partial derivative of a multivariate function.
+ ///
+ /// Input function handle.
+ /// Index of the independent variable for partial derivative.
+ /// Derivative order.
+ /// Function handle that evaluates partial derivative of input function at a fixed order.
+ public Func CreatePartialDerivativeFunctionHandle(Func f, int parameterIndex,
+ int order)
+ {
+ return x => EvaluatePartialDerivative(f, x, parameterIndex, order);
+ }
+
+ ///
+ /// Creates a function handle for the partial derivative of a vector multivariate function.
+ ///
+ /// Input function handle.
+ /// Index of the independent variable for partial derivative.
+ /// Derivative order.
+ /// Function handle that evaluates partial derivative of input function at fixed order.
+ public Func CreatePartialDerivativeFunctionHandle(Func[] f,
+ int parameterIndex,
+ int order)
+ {
+ return x => EvaluatePartialDerivative(f, x, parameterIndex, order);
+ }
+
+ ///
+ /// Evaluates the mixed partial derivative of variable order for multivariate functions.
+ ///
+ ///
+ /// This function recursively uses to evaluate mixed partial derivative.
+ /// Therefore, it is more efficient to call for higher order derivatives of
+ /// a single independent variable.
+ ///
+ /// Multivariate function handle.
+ /// Points at which to evaluate the derivative.
+ /// Vector of indices for the independent variables at descending derivative orders.
+ /// Highest order of differentiation.
+ /// Current function value at center.
+ /// Function mixed partial derivative at x of the specified order.
+ public double EvaluateMixedPartialDerivative(Func f, double[] x, int[] parameterIndex,
+ int order, double? currentValue = null)
+ {
+ if (parameterIndex.Length != order)
+ throw new ArgumentOutOfRangeException("parameterIndex",
+ "The number of parameters must match derivative order.");
+
+ if (order == 1)
+ return EvaluatePartialDerivative(f, x, parameterIndex[0], order, currentValue);
+
+ int reducedOrder = order - 1;
+ var reducedParameterIndex = new int[reducedOrder];
+ Array.Copy(parameterIndex, 0, reducedParameterIndex, 0, reducedOrder);
+
+ var points = new double[_points];
+ var currentParameterIndex = parameterIndex[order - 1];
+ var h = CalculateStepSize(_points, x[currentParameterIndex], order);
+
+ var xi = x[currentParameterIndex];
+ for (int i = 0; i < _points; i++)
+ {
+ x[currentParameterIndex] = xi + (i - Center)*h;
+ points[i] = EvaluateMixedPartialDerivative(f, x, reducedParameterIndex, reducedOrder);
+ }
+
+ // restore original value
+ x[currentParameterIndex] = xi;
+
+ // This will always be to the first order
+ return EvaluateDerivative(points, 1, h);
+ }
+
+ ///
+ /// Evaluates the mixed partial derivative of variable order for multivariate function arrays.
+ ///
+ ///
+ /// This function recursively uses to evaluate mixed partial derivative.
+ /// Therefore, it is more efficient to call for higher order derivatives of
+ /// a single independent variable.
+ ///
+ /// Multivariate function array handle.
+ /// Vector at which to evaluate the derivative.
+ /// Vector of indices for the independent variables at descending derivative orders.
+ /// Highest order of differentiation.
+ /// Current function value at center.
+ /// Function mixed partial derivatives at x of the specified order.
+ public double[] EvaluateMixedPartialDerivative(Func[] f, double[] x, int[] parameterIndex,
+ int order, double?[] currentValue = null)
+ {
+ var df = new double[f.Length];
+ for (int i = 0; i < f.Length; i++)
+ {
+ if(currentValue != null && currentValue[i].HasValue)
+ df[i] = EvaluateMixedPartialDerivative(f[i], x, parameterIndex, order, currentValue[i].Value);
+ else
+ df[i] = EvaluateMixedPartialDerivative(f[i], x, parameterIndex, order);
+ }
+
+ return df;
+ }
+
+ ///
+ /// Creates a function handle for the mixed partial derivative of a multivariate function.
+ ///
+ /// Input function handle.
+ /// Vector of indices for the independent variables at descending derivative orders.
+ /// Highest derivative order.
+ /// Function handle that evaluates the fixed mixed partial derivative of input function at fixed order.
+ public Func CreateMixedPartialDerivativeFunctionHandle(Func f,
+ int[] parameterIndex, int order)
+ {
+ return x => EvaluateMixedPartialDerivative(f, x, parameterIndex, order);
+ }
+
+ ///
+ /// Creates a function handle for the mixed partial derivative of a multivariate vector function.
+ ///
+ /// Input vector function handle.
+ /// Vector of indices for the independent variables at descending derivative orders.
+ /// Highest derivative order.
+ /// Function handle that evaluates the fixed mixed partial derivative of input function at fixed order.
+ public Func CreateMixedPartialDerivativeFunctionHandle(Func[] f,
+ int[] parameterIndex, int order)
+ {
+ return x => EvaluateMixedPartialDerivative(f, x, parameterIndex, order);
+ }
+
+ ///
+ /// Resets the evaluation counter.
+ ///
+ public void ResetEvaluations()
+ {
+ Evaluations = 0;
+ }
+
+ ///
+ /// Calculates machine epsilon - the smallest number that can be added to 1, yeilding a results different than 1.
+ /// This is also known as roundoff error.
+ ///
+ /// Machine epislon
+ public static double CalculateMachineEpsilon()
+ {
+ double eps = 1;
+
+ while ((1.0d + (eps / 2.0d)) > 1.0d)
+ eps /= 2.0d;
+
+ return eps;
+ }
+
+ private double[] CalculateStepSize(int points, double[] x, double order)
+ {
+ var h = new double[x.Length];
+ for (int i = 1; i < h.Length; i++)
+ h[i] = CalculateStepSize(points, x[i], order);
+
+ return h;
+ }
+
+ private double CalculateStepSize(int points, double x, double order)
+ {
+ // Step size relative to function input parameter
+ if (StepType == StepType.RelativeX)
+ {
+ StepSize = BaseStepSize*(1 + Math.Abs(x));
+ }
+ // Step size relative to function input parameter and order
+ else if (StepType == StepType.Relative)
+ {
+ var accuracy = points - order;
+ BaseStepSize = Math.Pow(Epsilon,(1/(accuracy + order)));
+ StepSize = BaseStepSize*(1 + Math.Abs(x));
+ }
+ // Do nothing for absolute step size.
+
+ return StepSize;
+ }
+ }
+}
diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj
index d88a75ae..166b4258 100644
--- a/src/Numerics/Numerics.csproj
+++ b/src/Numerics/Numerics.csproj
@@ -80,6 +80,8 @@
+
+
diff --git a/src/UnitTests/DifferentiationTests/FiniteDifferenceCoefficientsTests.cs b/src/UnitTests/DifferentiationTests/FiniteDifferenceCoefficientsTests.cs
new file mode 100644
index 00000000..dfeb0006
--- /dev/null
+++ b/src/UnitTests/DifferentiationTests/FiniteDifferenceCoefficientsTests.cs
@@ -0,0 +1,77 @@
+//
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+// Copyright (c) 2009-2015 Math.NET
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+//
+
+
+namespace MathNet.Numerics.UnitTests.DifferentiationTests
+{
+ using Differentiation;
+ using NUnit.Framework;
+
+ [TestFixture, Category("Differentiation")]
+ public class FiniteDifferenceCoefficientsTests
+ {
+ [Test]
+ public void CentralDifferenceFirstOrderThreePointTest()
+ {
+ double[] results = { -0.5, 0, 0.5 };
+ var finite = new FiniteDifferenceCoefficients(3);
+ var coeff = finite.GetCoefficients(1, 1);
+ Assert.AreEqual(results, coeff);
+ }
+
+ [Test]
+ public void CentralDifferenceSecondOrderFivePointsTest()
+ {
+ double[] results = { (double)-1 / 12, (double)4 / 3, (double)-5 / 2, (double)4 / 3, (double)-1 / 12 };
+ var finite = new FiniteDifferenceCoefficients(5);
+ var coeff = finite.GetCoefficients(2, 2);
+ for (int i = 0; i < coeff.Length; i++)
+ Assert.AreEqual(results[i], coeff[i]);
+ }
+
+ [Test]
+ public void ForwardDifferenceThirdOrderEightPointsTest()
+ {
+ double[] results = { (double)-967 / 120, (double)638 / 15, (double)-3929 / 40, (double)389 / 3,
+ (double)-2545 / 24, (double)268 / 5, (double)-1849 / 120, (double)29 / 15 };
+ var finite = new FiniteDifferenceCoefficients(8);
+ var coeff = finite.GetCoefficients(0, 3);
+ for (int i = 0; i < coeff.Length; i++)
+ Assert.AreEqual(results[i], coeff[i]);
+ }
+
+ [Test]
+ public void BackwardDifferenceThirdOrderFourPointsTest()
+ {
+ double[] results = { -1, 3, -3, 1 };
+ var finite = new FiniteDifferenceCoefficients(4);
+ var coeff = finite.GetCoefficients(3, 3);
+ for (int i = 0; i < coeff.Length; i++)
+ Assert.AreEqual(results[i], coeff[i]);
+
+ }
+ }
+}
diff --git a/src/UnitTests/DifferentiationTests/NumericalDerivativeTests.cs b/src/UnitTests/DifferentiationTests/NumericalDerivativeTests.cs
new file mode 100644
index 00000000..c8bf2c59
--- /dev/null
+++ b/src/UnitTests/DifferentiationTests/NumericalDerivativeTests.cs
@@ -0,0 +1,202 @@
+//
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+// Copyright (c) 2009-2015 Math.NET
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+//
+
+namespace MathNet.Numerics.UnitTests.DifferentiationTests
+{
+ using System;
+ using Differentiation;
+ using NUnit.Framework;
+
+ [TestFixture, Category("Differentiation")]
+ class NumericalDerivativeTests
+ {
+ [Test]
+ public void SinFirstDerivativeAtZeroTest()
+ {
+ Func f = Math.Sin;
+ var df = new NumericalDerivative();
+ Assert.AreEqual(1, df.EvaluateDerivative(f, 0, 1), 1e-10);
+ }
+
+ [Test]
+ public void CubicPolynomialThirdDerivativeAtAnyValTest()
+ {
+ Func f = x => 3 * Math.Pow(x, 3) + 2 * x - 6;
+ var df = new NumericalDerivative(5, 2);
+ Assert.AreEqual(18, df.EvaluateDerivative(f, 0, 3));
+ Assert.AreEqual(18, df.EvaluateDerivative(f, 10, 3));
+ df.Center = 0;
+ Assert.AreEqual(18, df.EvaluateDerivative(f, 0, 3));
+ Assert.AreEqual(18, df.EvaluateDerivative(f, 10, 3));
+ df.Center = 1;
+ Assert.AreEqual(18, df.EvaluateDerivative(f, 0, 3));
+ Assert.AreEqual(18, df.EvaluateDerivative(f, 10, 3));
+ df.Center = 2;
+ Assert.AreEqual(18, df.EvaluateDerivative(f, 0, 3));
+ Assert.AreEqual(18, df.EvaluateDerivative(f, 10, 3));
+ df.Center = 3;
+ Assert.AreEqual(18, df.EvaluateDerivative(f, 0, 3));
+ Assert.AreEqual(18, df.EvaluateDerivative(f, 10, 3));
+ df.Center = 4;
+ Assert.AreEqual(18, df.EvaluateDerivative(f, 0, 3));
+ Assert.AreEqual(18, df.EvaluateDerivative(f, 10, 3));
+ }
+
+ [Test]
+ public void CubicPolynomialFunctionValueTest()
+ {
+ Func f = x => 3 * Math.Pow(x, 3) + 2 * x - 6;
+ var current = f(2);
+ var df = new NumericalDerivative(3, 0);
+ Assert.AreEqual(38, df.EvaluateDerivative(f, 2, 1, current), 1e-8);
+ }
+
+ [Test]
+ public void CreateDerivativeFunctionHandleTest()
+ {
+ Func f = x => 3 * Math.Pow(x, 3) + 2 * x - 6;
+ var nd = new NumericalDerivative(5, 2);
+ var df = nd.CreateDerivativeFunctionHandle(f, 3);
+
+ Assert.AreEqual(18, df(0));
+
+ // Test new function with same nd class
+ Func f2 = x => 2 * Math.Pow(x, 3) + 2 * x - 6;
+ var df2 = nd.CreateDerivativeFunctionHandle(f2, 3);
+
+ Assert.AreEqual(12, df2(0));
+
+ // Original delegate not changed
+ Assert.AreEqual(18, df(0));
+
+ }
+
+ [Test]
+ public void ExponentialFunctionPartialDerivativeTest()
+ {
+ //Test Function
+ Func f = (x) => Math.Sin(x[0] * x[1]) + Math.Exp(-x[0] / 2) + x[1] / x[0];
+
+ //Analytical partial dfdx
+ Func dfdx =
+ (x) => Math.Cos(x[0] * x[1]) * x[1] - Math.Exp(-x[0] / 2) / 2 - x[1] / Math.Pow(x[0], 2);
+
+ //Analytical partial dfdy
+ Func dfdy = (x) => Math.Cos(x[0] * x[1]) * x[0] + 1 / x[0];
+
+ var df = new NumericalDerivative(3, 1);
+ var x1 = new double[] { 3, 3 };
+ Assert.AreEqual(dfdx(x1), df.EvaluatePartialDerivative(f, x1, 0, 1), 1e-8);
+
+ Assert.AreEqual(dfdy(x1), df.EvaluatePartialDerivative(f, x1, 1, 1), 1e-8);
+
+ var x2 = new double[] { 300, -50 };
+ df.StepType = StepType.Absolute;
+ Assert.AreEqual(dfdx(x2), df.EvaluatePartialDerivative(f, x2, 0, 1), 1e-5);
+ Assert.AreEqual(dfdy(x2), df.EvaluatePartialDerivative(f, x2, 1, 1), 1e-2);
+ }
+
+ [Test]
+ public void ExponentialFunctionPartialDerivativeCurrentValueTest()
+ {
+ //Test Function
+ Func f = (x) => Math.Sin(x[0] * x[1]) + Math.Exp(-x[0] / 2) + x[1] / x[0];
+
+ //Analytical partial dfdx
+ Func dfdx =
+ (x) => Math.Cos(x[0] * x[1]) * x[1] - Math.Exp(-x[0] / 2) / 2 - x[1] / Math.Pow(x[0], 2);
+
+ //Analytical partial dfdy
+ Func dfdy = (x) => Math.Cos(x[0] * x[1]) * x[0] + 1 / x[0];
+
+ // Current value
+ var x1 = new double[] { 3, 3 };
+ var current = f(x1);
+
+ var df = new NumericalDerivative(5, 2);
+ Assert.AreEqual(dfdx(x1), df.EvaluatePartialDerivative(f, x1, 0, 1, current), 1e-8);
+
+ Assert.AreEqual(dfdy(x1), df.EvaluatePartialDerivative(f, x1, 1, 1, current), 1e-8);
+ }
+
+ [Test]
+ public void RosenbrockFunctionMixedDerivativeOneVariableSecondOrderTest()
+ {
+ Func f = x => Math.Pow(1 - x[0], 2) + 100 * Math.Pow(x[1] - Math.Pow(x[0], 2), 2);
+ var df = new NumericalDerivative();
+ var x0 = new double[] { 2, 2 };
+ var parameterindex = new int[] { 0, 0 };
+ Assert.AreEqual(1602, df.EvaluatePartialDerivative(f, x0, 0, 1), 1e-6);
+ Assert.AreEqual(4002, df.EvaluateMixedPartialDerivative(f, x0, parameterindex, 2));
+ }
+
+ [Test]
+ public void RosenbrockFunctionMixedDerivativeTwoVariableSecondOrderTest()
+ {
+ Func f = x => Math.Pow(1 - x[0], 2) + 100 * Math.Pow(x[1] - Math.Pow(x[0], 2), 2);
+ var df = new NumericalDerivative();
+ var x0 = new double[] { 2, 2 };
+ var parameterIndex = new[] { 0, 1 };
+
+ Assert.AreEqual(-800, df.EvaluateMixedPartialDerivative(f, x0, parameterIndex, 2));
+ }
+
+ [Test]
+ public void VectorFunction1PartialDerivativeTest()
+ {
+ Func[] f =
+ {
+ (x) => Math.Pow(x[0],2) - 3*x[1],
+ (x) => x[1]*x[1] + 2*x[0]*x[1]
+ };
+
+ var x0 = new double[] { 2, 2 };
+ var g = new double[] { 4, 4 };
+
+ var df = new NumericalDerivative();
+ Assert.AreEqual(g, df.EvaluatePartialDerivative(f, x0, 0, 1));
+ Assert.AreEqual(new double[] { 2, 0 }, df.EvaluatePartialDerivative(f, x0, 0, 2));
+ Assert.AreEqual(new double[] { -3, 8 }, df.EvaluatePartialDerivative(f, x0, 1, 1));
+ }
+
+ [Test]
+ public void VectorFunctionMixedPartialDerivativeTest()
+ {
+ Func[] f =
+ {
+ (x) => Math.Pow(x[0],2) - 3*x[1],
+ (x) => x[1]*x[1] + 2*x[0]*x[1]
+ };
+
+ var x0 = new double[] { 2, 2 };
+
+ var df = new NumericalDerivative();
+ Assert.AreEqual(new double[] { 0, 2 }, df.EvaluateMixedPartialDerivative(f, x0, new int[] { 0, 1 }, 2));
+ Assert.AreEqual(new double[] { 0, 2 }, df.EvaluateMixedPartialDerivative(f, x0, new int[] { 1, 0 }, 2));
+
+ }
+ }
+}
diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj
index 1836e16d..1ac8eb41 100644
--- a/src/UnitTests/UnitTests.csproj
+++ b/src/UnitTests/UnitTests.csproj
@@ -77,6 +77,8 @@
+
+