Added numerical differentiation classes.

src/Numerics/Differentiation/FiniteDifferenceCoefficients.cs

@@ -0,0 +1,140 @@
+// <copyright file="FiniteDifferenceCoefficients.cs" company="Math.NET">
+// 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.
+// </copyright>
+
+using System;
+using MathNet.Numerics.LinearAlgebra.Double;
+
+namespace MathNet.Numerics.Differentiation
+{
+    /// <summary>
+    /// Class to calculate finite difference coefficients using Taylor series expansion method. 
+    /// <remarks>
+    /// <para>
+    /// 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.
+    /// </para> 
+    /// </remarks>
+    /// </summary>
+    public class FiniteDifferenceCoefficients
+    {
+        /// <summary>
+        /// Number of points for finite difference coefficients. Changing this value recalculates the coefficients table.
+        /// </summary>
+        public int Points
+        {
+            get { return _points; }
+            set
+            {
+                CalculateCoefficients(value);
+                _points = value;
+            }
+        }
+
+        private double[][,] _coefficients;
+        private int _points;
+
+        /// <summary>
+        /// Initializes a new instance of the <see cref="FiniteDifferenceCoefficients"/> class.
+        /// </summary>
+        /// <param name="points">Number of finite difference coefficients.</param>
+        public FiniteDifferenceCoefficients(int points)
+        {
+            Points = points;
+            CalculateCoefficients(Points);
+        }
+
+        /// <summary>
+        /// Gets the finite difference coefficients for a specified center and order.
+        /// </summary>
+        /// <param name="center">Current function position with respect to coefficients. Must be within point range.</param>
+        /// <param name="order">Order of finite difference coefficients.</param>
+        /// <returns>Vector of finite difference coefficients.</returns>
+        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;
+        }
+
+        /// <summary>
+        /// Gets the finite difference coefficients for all orders at a specified center.
+        /// </summary>
+        /// <param name="center">Current function position with respect to coefficients. Must be within point range.</param>
+        /// <returns>Rectangular array of coefficients, with columns specifing order.</returns>
+        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;
+        }
+    }
+}

src/Numerics/Differentiation/NumericalDerivative.cs

@@ -0,0 +1,483 @@
+// <copyright file="NumericalDerivative.cs" company="Math.NET">
+// 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.
+// </copyright>
+
+using System;
+using System.Linq;
+
+
+namespace MathNet.Numerics.Differentiation
+{
+    /// <summary>
+    /// Type of finite different step size.
+    /// </summary>
+    public enum StepType
+    {
+        /// <summary>
+        /// The absolute step size value will be used in numerical derivatives, regardless of order or function parameters.
+        /// </summary>
+        Absolute,
+
+        /// <summary>
+        /// 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.
+        /// </summary>
+        RelativeX,
+
+        /// <summary>
+        /// 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 <see cref="RelativeX"/> 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)).
+        /// </summary>
+        Relative
+    };
+
+    /// <summary>
+    /// Class to evaluate the numerical derivative of a function using finite difference approximations. 
+    /// Variable point and center methods can be initialized <seealso cref="FiniteDifferenceCoefficients"/>.
+    /// 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.
+    /// </summary>
+    public class NumericalDerivative
+    {
+        /// <summary>
+        /// 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 <see cref="Epsilon"/>. 
+        /// </summary>
+        /// <remarks>
+        /// 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. 
+        /// </remarks>
+        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));
+            }
+        }
+
+        /// <summary>
+        /// Sets and gets the base fininte difference step size. This assigned value to this parameter is only used if <see cref="StepType"/> is set to RelativeX.
+        /// However, if the StepType is Relative, it will contain the base step size computed from <see cref="Epsilon"/> based on the finite difference order.
+        /// </summary>
+        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));
+            }
+        }
+
+        /// <summary>
+        /// Sets and gets the base finite difference step size. This parameter is only used if <see cref="StepType"/> is set to Relative.
+        /// By default this is set to machine epsilon, from which <see cref="BaseStepSize"/> is computed.
+        /// </summary>
+        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));
+            }
+        }
+
+        /// <summary>
+        /// Sets and gets the location of the center point for the finite difference derivative.
+        /// </summary>
+        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;
+            }
+        }
+
+        /// <summary>
+        /// Number of times a function is evaluated for numerical derivatives.
+        /// </summary>
+        public int Evaluations { get; private set; }
+
+        /// <summary>
+        /// 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.
+        /// </summary>
+        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;
+
+        /// <summary>
+        /// Initializes a NumericalDerivative class with the default 3 point center difference method.
+        /// </summary>
+        public NumericalDerivative() : this(3, 1)
+        {
+        }
+
+        /// <summary>
+        /// Initialized a NumericalDerivative class.
+        /// </summary>
+        /// <param name="points">Number of points for finite difference derivatives.</param>
+        /// <param name="center">Location of the center with respect to other points. Value ranges from zero to points-1.</param>
+        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);
+        }
+
+        /// <summary>
+        /// Evaluates the derivative of equidistant points using the finite difference method.
+        /// </summary>
+        /// <param name="points">Vector of points StepSize apart.</param>
+        /// <param name="order">Derivative order.</param>
+        /// <param name="stepSize">Finite difference step size.</param>
+        /// <returns>Derivative of points of the specified order.</returns>
+        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;
+        }
+
+        /// <summary>
+        /// Evaluates the derivative of a scalar univariate function.
+        /// </summary>
+        /// <remarks>
+        /// Supplying the optional argument currentValue will reduce the number of function evaluations
+        /// required to calculate the finite difference derivative.
+        /// </remarks>
+        /// <param name="f">Function handle.</param>
+        /// <param name="x">Point at which to compute the derivative.</param>
+        /// <param name="order">Derivative order.</param>
+        /// <param name="currentValue">Current function value at center.</param>
+        /// <returns>Function derivative at x of the specified order.</returns>
+        public double EvaluateDerivative(Func<double, double> 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);
+        }
+
+        /// <summary>
+        /// Creates a function handle for the derivative of a scalar univariate function.
+        /// </summary>
+        /// <param name="f">Input function handle.</param>
+        /// <param name="order">Derivative order.</param>
+        /// <returns>Function handle that evaluates the derivative of input function at a fixed order.</returns>
+        public Func<double, double> CreateDerivativeFunctionHandle(Func<double, double> f, int order)
+        {
+            return x => EvaluateDerivative(f, x, order);
+        }
+
+        /// <summary>
+        /// Evaluates the partial derivative of a multivariate function.
+        /// </summary>
+        /// <param name="f">Multivariate function handle.</param>
+        /// <param name="x">Vector at which to evaluate the derivative.</param>
+        /// <param name="parameterIndex">Index of independent variable for partial derivative.</param>
+        /// <param name="order">Derivative order.</param>
+        /// <param name="currentValue">Current function value at center.</param>
+        /// <returns>Function partial derivative at x of the specified order.</returns>
+        public double EvaluatePartialDerivative(Func<double[], double> 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);
+        }
+
+        /// <summary>
+        /// Evaluates the partial derivatives of a multivariate function array.
+        /// </summary>
+        /// <remarks>
+        /// This function assumes the input vector x is of the correct length for f.
+        /// </remarks>
+        /// <param name="f">Multivariate vector function array handle.</param>
+        /// <param name="x">Vector at which to evaluate the derivatives.</param>
+        /// <param name="parameterIndex">Index of independent variable for partial derivative.</param>
+        /// <param name="order">Derivative order.</param>
+        /// <param name="currentValue">Current function value at center.</param>
+        /// <returns>Vector of functions partial derivatives at x of the specified order.</returns>
+        public double[] EvaluatePartialDerivative(Func<double[], double>[] 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;
+        }
+
+        /// <summary>
+        /// Creates a function handle for the partial derivative of a multivariate function.
+        /// </summary>
+        /// <param name="f">Input function handle.</param>
+        /// <param name="parameterIndex">Index of the independent variable for partial derivative.</param>
+        /// <param name="order">Derivative order.</param>
+        /// <returns>Function handle that evaluates partial derivative of input function at a fixed order.</returns>
+        public Func<double[], double> CreatePartialDerivativeFunctionHandle(Func<double[], double> f, int parameterIndex,
+                                                                            int order)
+        {
+            return x => EvaluatePartialDerivative(f, x, parameterIndex, order);
+        }
+
+        /// <summary>
+        /// Creates a function handle for the partial derivative of a vector multivariate function.
+        /// </summary>
+        /// <param name="f">Input function handle.</param>
+        /// <param name="parameterIndex">Index of the independent variable for partial derivative.</param>
+        /// <param name="order">Derivative order.</param>
+        /// <returns>Function handle that evaluates partial derivative of input function at fixed order.</returns>
+        public Func<double[], double[]> CreatePartialDerivativeFunctionHandle(Func<double[], double>[] f,
+                                                                            int parameterIndex,
+                                                                            int order)
+        {
+            return x => EvaluatePartialDerivative(f, x, parameterIndex, order);
+        }
+
+        /// <summary>
+        /// Evaluates the mixed partial derivative of variable order for multivariate functions.
+        /// </summary>
+        /// <remarks>
+        /// This function recursively uses <see cref="EvaluatePartialDerivative(Func&lt;double[], double&gt;, double[], int, int, double?)"/> to evaluate mixed partial derivative.
+        /// Therefore, it is more efficient to call <see cref="EvaluatePartialDerivative(Func&lt;double[], double&gt;, double[], int, int, double?)"/> for higher order derivatives of 
+        /// a single independent variable.
+        /// </remarks>
+        /// <param name="f">Multivariate function handle.</param>
+        /// <param name="x">Points at which to evaluate the derivative.</param>
+        /// <param name="parameterIndex">Vector of indices for the independent variables at descending derivative orders.</param>
+        /// <param name="order">Highest order of differentiation.</param>
+        /// <param name="currentValue">Current function value at center.</param>
+        /// <returns>Function mixed partial derivative at x of the specified order.</returns>
+        public double EvaluateMixedPartialDerivative(Func<double[], double> 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);
+        }
+
+        /// <summary>
+        /// Evaluates the mixed partial derivative of variable order for multivariate function arrays.
+        /// </summary>
+        /// <remarks>
+        /// This function recursively uses <see cref="EvaluatePartialDerivative(Func&lt;double[], double&gt;[], double[], int, int, double?[])"/> to evaluate mixed partial derivative.
+        /// Therefore, it is more efficient to call <see cref="EvaluatePartialDerivative(Func&lt;double[], double&gt;[], double[], int, int, double?[])"/> for higher order derivatives of 
+        /// a single independent variable.
+        /// </remarks>
+        /// <param name="f">Multivariate function array handle.</param>
+        /// <param name="x">Vector at which to evaluate the derivative.</param>
+        /// <param name="parameterIndex">Vector of indices for the independent variables at descending derivative orders.</param>
+        /// <param name="order">Highest order of differentiation.</param>
+        /// <param name="currentValue">Current function value at center.</param>
+        /// <returns>Function mixed partial derivatives at x of the specified order.</returns>
+        public double[] EvaluateMixedPartialDerivative(Func<double[], double>[] 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;
+        }
+
+        /// <summary>
+        /// Creates a function handle for the mixed partial derivative of a multivariate function.
+        /// </summary>
+        /// <param name="f">Input function handle.</param>
+        /// <param name="parameterIndex">Vector of indices for the independent variables at descending derivative orders.</param>
+        /// <param name="order">Highest derivative order.</param>
+        /// <returns>Function handle that evaluates the fixed mixed partial derivative of input function at fixed order.</returns>
+        public Func<double[], double> CreateMixedPartialDerivativeFunctionHandle(Func<double[], double> f,
+                                                                                 int[] parameterIndex, int order)
+        {
+            return x => EvaluateMixedPartialDerivative(f, x, parameterIndex, order);
+        }
+
+        /// <summary>
+        /// Creates a function handle for the mixed partial derivative of a multivariate vector function.
+        /// </summary>
+        /// <param name="f">Input vector function handle.</param>
+        /// <param name="parameterIndex">Vector of indices for the independent variables at descending derivative orders.</param>
+        /// <param name="order">Highest derivative order.</param>
+        /// <returns>Function handle that evaluates the fixed mixed partial derivative of input function at fixed order.</returns>
+        public Func<double[], double[]> CreateMixedPartialDerivativeFunctionHandle(Func<double[], double>[] f,
+                                                                                 int[] parameterIndex, int order)
+        {
+            return x => EvaluateMixedPartialDerivative(f, x, parameterIndex, order);
+        }
+
+        /// <summary>
+        /// Resets the evaluation counter.
+        /// </summary>
+        public void ResetEvaluations()
+        {
+            Evaluations = 0;
+        }
+
+        /// <summary>
+        /// 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.
+        /// </summary>
+        /// <returns>Machine epislon</returns>
+        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;
+        }
+    }
+}

src/Numerics/Numerics.csproj

@@ -80,6 +80,8 @@
     <Reference Include="System.Numerics" />
   </ItemGroup>
   <ItemGroup>
+    <Compile Include="Differentiation\FiniteDifferenceCoefficients.cs" />
+    <Compile Include="Differentiation\NumericalDerivative.cs" />
     <Compile Include="Distributions\Triangular.cs" />
     <Compile Include="Euclid.cs" />
     <Compile Include="Generate.cs" />

src/UnitTests/DifferentiationTests/FiniteDifferenceCoefficientsTests.cs

@@ -0,0 +1,77 @@
+// <copyright file="FiniteDifferenceCoefficientTests.cs" company="Math.NET">
+// 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.
+// </copyright>
+
+
+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]);
+
+        }
+    }
+}

src/UnitTests/DifferentiationTests/NumericalDerivativeTests.cs

@@ -0,0 +1,202 @@
+// <copyright file="NumericalDerivativeTests.cs" company="Math.NET">
+// 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.
+// </copyright>
+
+namespace MathNet.Numerics.UnitTests.DifferentiationTests
+{
+    using System;
+    using Differentiation;
+    using NUnit.Framework;
+
+    [TestFixture, Category("Differentiation")]
+    class NumericalDerivativeTests
+    {
+        [Test]
+        public void SinFirstDerivativeAtZeroTest()
+        {
+            Func<double, double> f = Math.Sin;
+            var df = new NumericalDerivative();
+            Assert.AreEqual(1, df.EvaluateDerivative(f, 0, 1), 1e-10);
+        }
+
+        [Test]
+        public void CubicPolynomialThirdDerivativeAtAnyValTest()
+        {
+            Func<double, double> 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<double, double> 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<double, double> 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<double, double> 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<double[], double> f = (x) => Math.Sin(x[0] * x[1]) + Math.Exp(-x[0] / 2) + x[1] / x[0];
+
+            //Analytical partial dfdx
+            Func<double[], double> 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<double[], double> 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<double[], double> f = (x) => Math.Sin(x[0] * x[1]) + Math.Exp(-x[0] / 2) + x[1] / x[0];
+
+            //Analytical partial dfdx
+            Func<double[], double> 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<double[], double> 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<double[], double> 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<double[], double> 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<double[], double>[] 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<double[], double>[] 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));
+
+        }
+    }
+}

src/UnitTests/UnitTests.csproj

@@ -77,6 +77,8 @@
     <Compile Include="ComplexTests\Complex32Test.TextHandling.cs" />
     <Compile Include="ComplexTests\ComplexTest.cs" />
     <Compile Include="ComplexTests\ComplexTest.TextHandling.cs" />
+    <Compile Include="DifferentiationTests\FiniteDifferenceCoefficientsTests.cs" />
+    <Compile Include="DifferentiationTests\NumericalDerivativeTests.cs" />
     <Compile Include="DistanceTests.cs" />
     <Compile Include="DistributionTests\CommonDistributionTests.cs" />
     <Compile Include="DistributionTests\Continuous\BetaTests.cs" />