Numerical Differentiation

Binder Notebook

Numerical differentiation is used to estimate the derivative of a mathematical function using values of the function and perhaps other knowledge about the function.

Three-Point Differentiation

FSharp.Stats implements a three point differentiation method. This method takes a set of values and their function values. For a given value xT of the set, one defines three other points which should be considered to calculate the differential at xT. Here follows a small snippet.

First, we create our data. In this case our function is f(x) = x ^ 3.

open Plotly.NET
open FSharp.Stats

// x data
let xs = Array.init 100 (fun x -> float x / 8.)
// y data
let ys = Array.map (fun x -> x ** 3.) xs

Now we apply the threePointDifferentiation to every point starting with the second and ending with the second last. We can't do it for every point because for every point we need to have three other points in close proximity.

let y's = 
    [|
    for i = 1 to xs.Length - 2 do
        yield xs.[i],Integration.Differentiation.differentiateThreePoint xs ys i (i-1) (i) (i+1)
    |]

We compare the resulting values with the values of the known differential f'(x) = 3x^2, here called g(x)

open Plotly.NET

let comparisonChart = 
    [
    Chart.Point(xs,ys,Name = "f(x)")
    Chart.Point(y's,Name = "f'(x)")
    Chart.Point(Array.map (fun x -> x,(x ** 2.) * 3.) xs,Name = "g(x)")
    ]
    |> Chart.combine
    |> Chart.withTemplate ChartTemplates.lightMirrored

Two-Point Differentiation

To calculate the approximation for the derivative, a Two-Point Differentiation calculates the difference of f(x) at x and f(x) at x+h and correlates it to h. This will give better approximations the smaller h is. The function uses two different mathematical approaches to decrease the error, one for h > 1. and one for h < 1.


Data model


open FSharp.Stats.Integration.Differentiation.TwoPointDifferentiation
open FSharp.Stats.Integration.Differentiation

let testFunction x = x**2. 

let test1 = differentiate 0.5 testFunction 2.
//Result for test1 is: 4.00

let test2 = differentiate 3. testFunction 2.
//Result for test2 is: 7.00

let test3 = differentiate 0.1 testFunction 2.
//Result for test3 is: 4.00
  • The correct result for test1 (f(x) = x**2.) is assumed to be 4.
  • With a higher h (for test2 h = 3., compared to the 0.5 used in test1) the approximation error increases.

You can try and find an optimal h - value with the "differentiateOptimalHBy" - function. This function uses the first h-value it assumes to give good results for the numerical differentiation calculation. Due to this, possible error due to float precision is avoided. In the following example this is not really necessary, as values are quite big.

let hArray = [|0.1 .. 0.1 .. 2.|]

let test4 = differentiateOptimalHBy hArray testFunction 2.
//Result for test4 is: 4.00

If you want to use a presuggested hArray then you can use the "differentiateOptimalH" function. This function uses an array from 0.01 to 5e^-100 in [|0.01; 0.005; 0.001; 0.0005; 0.0001 ..|]-increments as hArray.

let test5 = differentiateOptimalH testFunction 2. 
//Result for test5 is: 4.00
namespace Plotly
namespace Plotly.NET
module Defaults from Plotly.NET
<summary> Contains mutable global default values. Changing these values will apply the default values to all consecutive Chart generations. </summary>
val mutable DefaultDisplayOptions: DisplayOptions
Multiple items
type DisplayOptions = inherit DynamicObj new: unit -> DisplayOptions static member addAdditionalHeadTags: additionalHeadTags: XmlNode list -> (DisplayOptions -> DisplayOptions) static member addDescription: description: XmlNode list -> (DisplayOptions -> DisplayOptions) static member combine: first: DisplayOptions -> second: DisplayOptions -> DisplayOptions static member getAdditionalHeadTags: displayOpts: DisplayOptions -> XmlNode list static member getDescription: displayOpts: DisplayOptions -> XmlNode list static member getPlotlyReference: displayOpts: DisplayOptions -> PlotlyJSReference static member init: ?AdditionalHeadTags: XmlNode list * ?Description: XmlNode list * ?PlotlyJSReference: PlotlyJSReference -> DisplayOptions static member initCDNOnly: unit -> DisplayOptions ...

--------------------
new: unit -> DisplayOptions
static member DisplayOptions.init: ?AdditionalHeadTags: Giraffe.ViewEngine.HtmlElements.XmlNode list * ?Description: Giraffe.ViewEngine.HtmlElements.XmlNode list * ?PlotlyJSReference: PlotlyJSReference -> DisplayOptions
type PlotlyJSReference = | CDN of string | Full | Require of string | NoReference
<summary> Sets how plotly is referenced in the head of html docs. </summary>
union case PlotlyJSReference.NoReference: PlotlyJSReference
Multiple items
namespace FSharp

--------------------
namespace Microsoft.FSharp
namespace FSharp.Stats
val xs: float[]
Multiple items
module Array from FSharp.Stats
<summary> Module to compute common statistical measure on array </summary>

--------------------
module Array from Microsoft.FSharp.Collections
<summary>Contains operations for working with arrays.</summary>
<remarks> See also <a href="https://docs.microsoft.com/dotnet/fsharp/language-reference/arrays">F# Language Guide - Arrays</a>. </remarks>


--------------------
type Array = new: unit -> Array static member geomspace: start: float * stop: float * num: int * ?IncludeEndpoint: bool -> float array static member linspace: start: float * stop: float * num: int * ?IncludeEndpoint: bool -> float[]

--------------------
new: unit -> Array
val init: count: int -> initializer: (int -> 'T) -> 'T[]
<summary>Creates an array given the dimension and a generator function to compute the elements.</summary>
<param name="count">The number of elements to initialize.</param>
<param name="initializer">The function to generate the initial values for each index.</param>
<returns>The created array.</returns>
<exception cref="T:System.ArgumentException">Thrown when count is negative.</exception>
<example id="init-1"><code lang="fsharp"> Array.init 4 (fun v -&gt; v + 5) </code> Evaluates to <c>[| 5; 6; 7; 8 |]</c></example>
<example id="init-2"><code lang="fsharp"> Array.init -5 (fun v -&gt; v + 5) </code> Throws <c>ArgumentException</c></example>
val x: int
Multiple items
val float: value: 'T -> float (requires member op_Explicit)
<summary>Converts the argument to 64-bit float. This is a direct conversion for all primitive numeric types. For strings, the input is converted using <c>Double.Parse()</c> with InvariantCulture settings. Otherwise the operation requires an appropriate static conversion method on the input type.</summary>
<param name="value">The input value.</param>
<returns>The converted float</returns>
<example id="float-example"><code lang="fsharp"></code></example>


--------------------
[<Struct>] type float = System.Double
<summary>An abbreviation for the CLI type <see cref="T:System.Double" />.</summary>
<category>Basic Types</category>


--------------------
type float<'Measure> = float
<summary>The type of double-precision floating point numbers, annotated with a unit of measure. The unit of measure is erased in compiled code and when values of this type are analyzed using reflection. The type is representationally equivalent to <see cref="T:System.Double" />.</summary>
<category index="6">Basic Types with Units of Measure</category>
val ys: float[]
val map: mapping: ('T -> 'U) -> array: 'T[] -> 'U[]
<summary>Builds a new array whose elements are the results of applying the given function to each of the elements of the array.</summary>
<param name="mapping">The function to transform elements of the array.</param>
<param name="array">The input array.</param>
<returns>The array of transformed elements.</returns>
<exception cref="T:System.ArgumentNullException">Thrown when the input array is null.</exception>
<example id="map-1"><code lang="fsharp"> let inputs = [| "a"; "bbb"; "cc" |] inputs |&gt; Array.map (fun x -&gt; x.Length) </code> Evaluates to <c>[| 1; 3; 2 |]</c></example>
val x: float
val y's: (float * float)[]
val i: int
property System.Array.Length: int with get
<summary>Gets the total number of elements in all the dimensions of the <see cref="T:System.Array" />.</summary>
<exception cref="T:System.OverflowException">The array is multidimensional and contains more than <see cref="F:System.Int32.MaxValue" /> elements.</exception>
<returns>The total number of elements in all the dimensions of the <see cref="T:System.Array" />; zero if there are no elements in the array.</returns>
namespace FSharp.Stats.Integration
module Differentiation from FSharp.Stats.Integration
<summary> In numerical analysis, numerical differentiation describes algorithms for estimating the derivative of a mathematical function using values of the function and perhaps other knowledge about the function. </summary>
val differentiateThreePoint: xValues: float[] -> yValues: float[] -> idxT: int -> idx0: int -> idx1: int -> idx2: int -> float
<summary>Three-Point Differentiation Helper</summary>
<remarks></remarks>
<param name="xValues">Sample Points t</param>
<param name="yValues">Sample Values x(t)</param>
<param name="idxT">Index of the point of the differentiation.</param>
<param name="idx0">idx0 Index of the first sample.</param>
<param name="idx1">Index of the second sample.</param>
<param name="idx2">Index of the third sample.</param>
<returns></returns>
<example><code></code></example>
val comparisonChart: GenericChart.GenericChart
type Chart = static member AnnotatedHeatmap: zData: seq<#seq<'a1>> * annotationText: seq<#seq<string>> * ?Name: string * ?ShowLegend: bool * ?Opacity: float * ?X: seq<'a3> * ?MultiX: seq<seq<'a3>> * ?XGap: int * ?Y: seq<'a4> * ?MultiY: seq<seq<'a4>> * ?YGap: int * ?Text: 'a5 * ?MultiText: seq<'a5> * ?ColorBar: ColorBar * ?ColorScale: Colorscale * ?ShowScale: bool * ?ReverseScale: bool * ?ZSmooth: SmoothAlg * ?Transpose: bool * ?UseWebGL: bool * ?ReverseYAxis: bool * ?UseDefaults: bool -> GenericChart (requires 'a1 :> IConvertible and 'a3 :> IConvertible and 'a4 :> IConvertible and 'a5 :> IConvertible) + 1 overload static member Area: x: seq<#IConvertible> * y: seq<#IConvertible> * ?ShowMarkers: bool * ?Name: string * ?ShowLegend: bool * ?Opacity: float * ?MultiOpacity: seq<float> * ?Text: 'a2 * ?MultiText: seq<'a2> * ?TextPosition: TextPosition * ?MultiTextPosition: seq<TextPosition> * ?MarkerColor: Color * ?MarkerColorScale: Colorscale * ?MarkerOutline: Line * ?MarkerSymbol: MarkerSymbol * ?MultiMarkerSymbol: seq<MarkerSymbol> * ?Marker: Marker * ?LineColor: Color * ?LineColorScale: Colorscale * ?LineWidth: float * ?LineDash: DrawingStyle * ?Line: Line * ?AlignmentGroup: string * ?OffsetGroup: string * ?StackGroup: string * ?Orientation: Orientation * ?GroupNorm: GroupNorm * ?FillColor: Color * ?FillPatternShape: PatternShape * ?FillPattern: Pattern * ?UseWebGL: bool * ?UseDefaults: bool -> GenericChart (requires 'a2 :> IConvertible) + 1 overload static member Bar: values: seq<#IConvertible> * ?Keys: seq<'a1> * ?MultiKeys: seq<seq<'a1>> * ?Name: string * ?ShowLegend: bool * ?Opacity: float * ?MultiOpacity: seq<float> * ?Text: 'a2 * ?MultiText: seq<'a2> * ?MarkerColor: Color * ?MarkerColorScale: Colorscale * ?MarkerOutline: Line * ?MarkerPatternShape: PatternShape * ?MultiMarkerPatternShape: seq<PatternShape> * ?MarkerPattern: Pattern * ?Marker: Marker * ?Base: #IConvertible * ?Width: 'a4 * ?MultiWidth: seq<'a4> * ?TextPosition: TextPosition * ?MultiTextPosition: seq<TextPosition> * ?UseDefaults: bool -> GenericChart (requires 'a1 :> IConvertible and 'a2 :> IConvertible and 'a4 :> IConvertible) + 1 overload static member BoxPlot: ?X: seq<'a0> * ?MultiX: seq<seq<'a0>> * ?Y: seq<'a1> * ?MultiY: seq<seq<'a1>> * ?Name: string * ?ShowLegend: bool * ?Text: 'a2 * ?MultiText: seq<'a2> * ?FillColor: Color * ?MarkerColor: Color * ?Marker: Marker * ?Opacity: float * ?WhiskerWidth: float * ?BoxPoints: BoxPoints * ?BoxMean: BoxMean * ?Jitter: float * ?PointPos: float * ?Orientation: Orientation * ?OutlineColor: Color * ?OutlineWidth: float * ?Outline: Line * ?AlignmentGroup: string * ?OffsetGroup: string * ?Notched: bool * ?NotchWidth: float * ?QuartileMethod: QuartileMethod * ?UseDefaults: bool -> GenericChart (requires 'a0 :> IConvertible and 'a1 :> IConvertible and 'a2 :> IConvertible) + 2 overloads static member Bubble: x: seq<#IConvertible> * y: seq<#IConvertible> * sizes: seq<int> * ?Name: string * ?ShowLegend: bool * ?Opacity: float * ?MultiOpacity: seq<float> * ?Text: 'a2 * ?MultiText: seq<'a2> * ?TextPosition: TextPosition * ?MultiTextPosition: seq<TextPosition> * ?MarkerColor: Color * ?MarkerColorScale: Colorscale * ?MarkerOutline: Line * ?MarkerSymbol: MarkerSymbol * ?MultiMarkerSymbol: seq<MarkerSymbol> * ?Marker: Marker * ?LineColor: Color * ?LineColorScale: Colorscale * ?LineWidth: float * ?LineDash: DrawingStyle * ?Line: Line * ?AlignmentGroup: string * ?OffsetGroup: string * ?StackGroup: string * ?Orientation: Orientation * ?GroupNorm: GroupNorm * ?UseWebGL: bool * ?UseDefaults: bool -> GenericChart (requires 'a2 :> IConvertible) + 1 overload static member Candlestick: ``open`` : seq<#IConvertible> * high: seq<#IConvertible> * low: seq<#IConvertible> * close: seq<#IConvertible> * ?X: seq<'a4> * ?MultiX: seq<seq<'a4>> * ?Name: string * ?ShowLegend: bool * ?Opacity: float * ?Text: 'a5 * ?MultiText: seq<'a5> * ?Line: Line * ?IncreasingColor: Color * ?Increasing: FinanceMarker * ?DecreasingColor: Color * ?Decreasing: FinanceMarker * ?WhiskerWidth: float * ?ShowXAxisRangeSlider: bool * ?UseDefaults: bool -> GenericChart (requires 'a4 :> IConvertible and 'a5 :> IConvertible) + 2 overloads static member Column: values: seq<#IConvertible> * ?Keys: seq<'a1> * ?MultiKeys: seq<seq<'a1>> * ?Name: string * ?ShowLegend: bool * ?Opacity: float * ?MultiOpacity: seq<float> * ?Text: 'a2 * ?MultiText: seq<'a2> * ?MarkerColor: Color * ?MarkerColorScale: Colorscale * ?MarkerOutline: Line * ?MarkerPatternShape: PatternShape * ?MultiMarkerPatternShape: seq<PatternShape> * ?MarkerPattern: Pattern * ?Marker: Marker * ?Base: #IConvertible * ?Width: 'a4 * ?MultiWidth: seq<'a4> * ?TextPosition: TextPosition * ?MultiTextPosition: seq<TextPosition> * ?UseDefaults: bool -> GenericChart (requires 'a1 :> IConvertible and 'a2 :> IConvertible and 'a4 :> IConvertible) + 1 overload static member Contour: zData: seq<#seq<'a1>> * ?Name: string * ?ShowLegend: bool * ?Opacity: float * ?X: seq<'a2> * ?MultiX: seq<seq<'a2>> * ?Y: seq<'a3> * ?MultiY: seq<seq<'a3>> * ?Text: 'a4 * ?MultiText: seq<'a4> * ?ColorBar: ColorBar * ?ColorScale: Colorscale * ?ShowScale: bool * ?ReverseScale: bool * ?Transpose: bool * ?ContourLineColor: Color * ?ContourLineDash: DrawingStyle * ?ContourLineSmoothing: float * ?ContourLine: Line * ?ContoursColoring: ContourColoring * ?ContoursOperation: ConstraintOperation * ?ContoursType: ContourType * ?ShowContourLabels: bool * ?ContourLabelFont: Font * ?Contours: Contours * ?FillColor: Color * ?NContours: int * ?UseDefaults: bool -> GenericChart (requires 'a1 :> IConvertible and 'a2 :> IConvertible and 'a3 :> IConvertible and 'a4 :> IConvertible) static member Funnel: x: seq<#IConvertible> * y: seq<#IConvertible> * ?Name: string * ?ShowLegend: bool * ?Opacity: float * ?Width: float * ?Offset: float * ?Text: 'a2 * ?MultiText: seq<'a2> * ?TextPosition: TextPosition * ?MultiTextPosition: seq<TextPosition> * ?Orientation: Orientation * ?AlignmentGroup: string * ?OffsetGroup: string * ?MarkerColor: Color * ?MarkerOutline: Line * ?Marker: Marker * ?TextInfo: TextInfo * ?ConnectorLineColor: Color * ?ConnectorLineStyle: DrawingStyle * ?ConnectorFillColor: Color * ?ConnectorLine: Line * ?Connector: FunnelConnector * ?InsideTextFont: Font * ?OutsideTextFont: Font * ?UseDefaults: bool -> GenericChart (requires 'a2 :> IConvertible) static member Heatmap: zData: seq<#seq<'a1>> * ?X: seq<'a2> * ?MultiX: seq<seq<'a2>> * ?Y: seq<'a3> * ?MultiY: seq<seq<'a3>> * ?Name: string * ?ShowLegend: bool * ?Opacity: float * ?XGap: int * ?YGap: int * ?Text: 'a4 * ?MultiText: seq<'a4> * ?ColorBar: ColorBar * ?ColorScale: Colorscale * ?ShowScale: bool * ?ReverseScale: bool * ?ZSmooth: SmoothAlg * ?Transpose: bool * ?UseWebGL: bool * ?ReverseYAxis: bool * ?UseDefaults: bool -> GenericChart (requires 'a1 :> IConvertible and 'a2 :> IConvertible and 'a3 :> IConvertible and 'a4 :> IConvertible) + 1 overload ...
static member Chart.Point: xy: seq<#System.IConvertible * #System.IConvertible> * ?Name: string * ?ShowLegend: bool * ?Opacity: float * ?MultiOpacity: seq<float> * ?Text: 'c * ?MultiText: seq<'c> * ?TextPosition: StyleParam.TextPosition * ?MultiTextPosition: seq<StyleParam.TextPosition> * ?MarkerColor: Color * ?MarkerColorScale: StyleParam.Colorscale * ?MarkerOutline: Line * ?MarkerSymbol: StyleParam.MarkerSymbol * ?MultiMarkerSymbol: seq<StyleParam.MarkerSymbol> * ?Marker: TraceObjects.Marker * ?AlignmentGroup: string * ?OffsetGroup: string * ?StackGroup: string * ?Orientation: StyleParam.Orientation * ?GroupNorm: StyleParam.GroupNorm * ?UseWebGL: bool * ?UseDefaults: bool -> GenericChart.GenericChart (requires 'c :> System.IConvertible)
static member Chart.Point: x: seq<#System.IConvertible> * y: seq<#System.IConvertible> * ?Name: string * ?ShowLegend: bool * ?Opacity: float * ?MultiOpacity: seq<float> * ?Text: 'c * ?MultiText: seq<'c> * ?TextPosition: StyleParam.TextPosition * ?MultiTextPosition: seq<StyleParam.TextPosition> * ?MarkerColor: Color * ?MarkerColorScale: StyleParam.Colorscale * ?MarkerOutline: Line * ?MarkerSymbol: StyleParam.MarkerSymbol * ?MultiMarkerSymbol: seq<StyleParam.MarkerSymbol> * ?Marker: TraceObjects.Marker * ?AlignmentGroup: string * ?OffsetGroup: string * ?StackGroup: string * ?Orientation: StyleParam.Orientation * ?GroupNorm: StyleParam.GroupNorm * ?UseWebGL: bool * ?UseDefaults: bool -> GenericChart.GenericChart (requires 'c :> System.IConvertible)
static member Chart.combine: gCharts: seq<GenericChart.GenericChart> -> GenericChart.GenericChart
static member Chart.withTemplate: template: Template -> (GenericChart.GenericChart -> GenericChart.GenericChart)
module ChartTemplates from Plotly.NET
val lightMirrored: Template
module GenericChart from Plotly.NET
<summary> Module to represent a GenericChart </summary>
val toChartHTML: gChart: GenericChart.GenericChart -> string
module TwoPointDifferentiation from FSharp.Stats.Integration.Differentiation
<summary> A two-point estimation is to compute the slope of a nearby secant line through two points. This gives an approximations of f'(x) at x respectively to two points "x and x+h"/"x-h and x+h"(depending on the used algorithm) of the function f. Choosing a small number h, h represents a small change in x, and it can be either positive or negative. </summary>
val testFunction: x: float -> float
val test1: float
val differentiate: h: float -> f: (float -> float) -> x: float -> float
<summary>Returns the approximation of f'(x) at x by calculating the two point differentiation.</summary>
<remarks></remarks>
<param name="h">window for the difference calculation</param>
<param name="f">f is the function for which to calculate numerical differentiation.</param>
<param name="x">x is the point at which the difference between "x and x+h"/"x-h and x+h" is calculated.</param>
<returns></returns>
<example><code></code></example>
val test2: float
val test3: float
val hArray: float[]
val test4: float
val differentiateOptimalHBy: hArr: float[] -> f: (float -> float) -> x: float -> float
<summary>Finds optimal h from all values given in hArr and calculates "differentiate" -function.</summary>
<remarks></remarks>
<param name="hArr"></param>
<param name="f">function for which numerical differentiation is calculated.</param>
<param name="x">x is the point at which numerical differentiation is calculated.</param>
<returns>Returns the approximation of f'(x) at x by calculating the two point differentiation.</returns>
<example><code></code></example>
val test5: float
val differentiateOptimalH: f: (float -> float) -> x: float -> float