In this post I develop a class for numerically differentiating a function using equally spaced values of that function. An overview of this can be found on Wikipedia at Numerical Differentiation.
The Taylor series of a function expanded about x0 can be written as shown below.
f(x) = (f0(x0)/0!)(x-x0)0 + (f1(x0)/1!)(x-x0)1 + (f2(x0)/2!)(x-x0)2 + …. + (fn(x0)/n!) (x-x0)n
For the differentiation, I will have nNumPoints function values each equally spaced along the x axis with a spacing distance of h.
Consider the case of n points. In this case the following Taylor series can be written.
In this case, the terms on the left, which are the equally spaced function values, are known. What I want to solve for is the derivatives that appear on the right side of these equations.
Converting this to matrix form gives the following matrix equation.
F = T D
Where F is the vector shown below.
f(x0+ 0h)
f(x0 + 1h)
f(x0 + 2h)
+
…
+
f(x0 + nh)
D is the vector that I want to solve for shown below.
h0f(x0)
h1f1(x0)
h2f2(x0)
+
…
+
hnfn(x0)
And T is the matrix below.
00/0! 01/1! 02/2! … 0n/n!
10/0! 11/1! 12/2! + … 1n/n!
20/0! 21/1! 22/2! + … 2n/n!
+
…
+
n0/0! n1/1! n2/2! … nn/n!
Inversion of the above matrix and multiplication of it by the functional values gives the derivatives.
The code to do this is shown below. The class MatrixNumeric used in the code is from a previous post creating a Matrix class.
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Diagnostics;
using NumericalMethods.FourthBlog;
namespace NumericalMethods
{
public class Derivatives
{
// _coefficients is the array of differential coefficients matrices.
// The index corresponds to the position from the left edge
// of the points.
// I.e _coefficients[0] is for a matrix with three points in it corresponds to
// the left most point.
// The coefficients of the derivatives go down by row. I.e. the first row
// is the functional value, the second row is for the first derivative of the functional
// value, the third row is the second derivative of the functional value.
// The columns correspond to the points themselves.
private Matrix[] _coefficients;
private Derivatives()
{
}
public Derivatives(int numberOfPoints)
: this()
{
SolveCoefs(numberOfPoints);
}
public void SolveCoefs(int numberOfPoints)
{
_coefficients = new Matrix[numberOfPoints];
for (int i = 0; i < numberOfPoints; i++)
{
Matrix deltsMatrix = new Matrix(numberOfPoints, numberOfPoints);
for (int j = 0; j < numberOfPoints; j++)
{
double delt = (double)(j - i);
double HTerm = 1.0;
for (int k = 0; k < numberOfPoints; k++)
{
deltsMatrix[j, k] = HTerm / Factorial(k);
HTerm *= delt;
}
}
_coefficients[i] = deltsMatrix.Invert();
double numPointsFactorial = Factorial(numberOfPoints);
for (int j = 0; j < numberOfPoints; j++)
{
for (int k = 0; k < numberOfPoints; k++)
{
_coefficients[i][j, k] = (Math.Round(_coefficients[i][j, k] * numPointsFactorial)) / numPointsFactorial;
}
}
}
}
private static double Factorial(int value)
{
double result = 1.0;
for (int i = 1; i <= value; i++)
{
result *= (double)i;
}
return result;
}
/// <summary>
/// Computes the derivative of a function.
/// </summary>
/// <param name="points">Equally spaced function value points</param>
/// <param name="order">The order of the derivative to take</param>
/// <param name="variablePosition">The position in the array of function values to take the derivative at.</param>
/// <param name="step">The x axis step size.</param>
/// <returns></returns>
public double ComputeDerivative(double[] points, int order, int variablePosition, double step)
{
Debug.Assert(points.Length == _coefficients.Length);
Debug.Assert(order < _coefficients.Length);
double result = 0.0;
for (int i = 0; i < _coefficients.Length; i++)
{
result += _coefficients[variablePosition][order, i] * points[i];
}
result /= Math.Pow(step, order);
return result;
}
public double ComputePartialDerivative(Func<double> function, Parameter parameter, int order)
{
int numberOfPoints = _coefficients.Length;
double result = 0.0;
double originalValue = parameter;
double[] points = new double[numberOfPoints];
double derivativeStepSize = parameter.DerivativeStepSize;
int centerPoint = (numberOfPoints - 1) / 2;
for (int i = 0; i < numberOfPoints; i++)
{
parameter.Value = originalValue + ((double)(i - centerPoint)) * derivativeStepSize;
points[i] = function();
}
result = ComputeDerivative(points, order, centerPoint, derivativeStepSize);
parameter.Value = originalValue;
return result;
}
public double ComputePartialDerivative(Func<double> function, Parameter parameter, int order, double currentFunctionValue)
{
int numberOfPoints = _coefficients.Length;
double result = 0.0;
double originalValue = parameter;
double[] points = new double[numberOfPoints];
double derivativeStepSize = parameter.DerivativeStepSize;
int centerPoint = (numberOfPoints - 1) / 2;
for (int i = 0; i < numberOfPoints; i++)
{
if (i != centerPoint)
{
parameter.Value = originalValue + ((double)(i - centerPoint)) * derivativeStepSize;
points[i] = function();
}
else
{
points[i] = currentFunctionValue;
}
}
result = ComputeDerivative(points, order, centerPoint, derivativeStepSize);
parameter.Value = originalValue;
return result;
}
public double[] ComputePartialDerivatives(Func<double> function, Parameter parameter, int[] derivativeOrders)
{
int numberOfPoints = _coefficients.Length;
double[] result = new double[derivativeOrders.Length];
double originalValue = parameter;
double[] points = new double[numberOfPoints];
double derivativeStepSize = parameter.DerivativeStepSize;
int centerPoint = (numberOfPoints - 1) / 2;
for (int i = 0; i < numberOfPoints; i++)
{
parameter.Value = originalValue + ((double)(i - centerPoint)) * derivativeStepSize;
points[i] = function();
}
for (int i = 0; i < derivativeOrders.Length; i++)
{
result[i] = ComputeDerivative(points, derivativeOrders[i], centerPoint, derivativeStepSize);
}
parameter.Value = originalValue;
return result;
}
public double[] ComputePartialDerivatives(Func<double> function, Parameter parameter, int[] derivativeOrders, double currentFunctionValue)
{
int numberOfPoints = _coefficients.Length;
double[] result = new double[derivativeOrders.Length];
double originalValue = parameter;
double[] points = new double[numberOfPoints];
double derivativeStepSize = parameter.DerivativeStepSize;
int centerPoint = (numberOfPoints - 1) / 2;
for (int i = 0; i < numberOfPoints; i++)
{
if (i != centerPoint)
{
parameter.Value = originalValue + ((double)(i - centerPoint)) * derivativeStepSize;
points[i] = function();
}
else
{
points[i] = currentFunctionValue;
}
}
for (int i = 0; i < derivativeOrders.Length; i++)
{
result[i] = ComputeDerivative(points, derivativeOrders[i], centerPoint, derivativeStepSize);
}
parameter.Value = originalValue;
return result;
}
}
}
Subscribe