Example 4: Central Differences

This example uses the same data as in the One-Sided Differences example. Agreement should be approximately the two-thirds power of machine precision. That agreement is achieved here. Generally this is the most accuracy one can expect using central divided differences. Note that using central differences requires essentially twice the number of evaluations of the function compared with obtaining one-sided differences. This can be a significant issue for functions that are expensive to evaluate. This example shows how to override EvaluateF . 
using System;
using Imsl.Math;


public class NumericalDerivativesEx4 : NumericalDerivatives
{
    static int m = 1, n = 2;
    static double a, b, c, v = 0.0;


    class NumericalDerivativesFcn : NumericalDerivatives.IFunction
    {
        public double[] F(int varIndex, double[] y) 
        {
            return new double[m];
        }
    }
    
    
    public NumericalDerivativesEx4(NumericalDerivatives.IFunction fcn) : 
        base(fcn) 
    {
    }
    
    
    // Override EvaluateF.
    public override double[] EvaluateF(int varIndex, double[] y)  
    {
        double[] valueF = new double[m];
     
        valueF[0] = a * Math.Exp(b * y[0]) + c * y[0] * y[1] * y[1];
        return valueF;
    }
    
    
    public static void Main(String[] args) 
    {
        NumericalDerivatives.DifferencingMethod[] options = 
            new NumericalDerivatives.DifferencingMethod[n];
        double u;
        double[] y = new double[n];
        double[] scale = new double[n];
        double[,] actual = new double[m, n];
        double[] re = new double[2];
    
        // Define data and point of evaluation:
        a = 2.5e6;
        b = 3.4e0;
        c = 4.5e0;
        y[0] = 2.1e0;
        y[1] = 3.2e0;
        
        // Machine precision, for measuring errors
        u = 2.220446049250313e-016;
        v = Math.Pow(3.0e0 * u, 2.0e0 / 3.0e0);
        
        // Set scaling:
        scale[0] = 1.0e0;
        // Increase scale to account for large value of a.
        scale[1] = 8.0e3;
        
        // Compute true values of partials.
        actual[0, 0] = a * b * Math.Exp(b * y[0]) + c * y[1] * y[1];
        actual[0, 1] = 2 * c * y[0] * y[1];
        
        options[0] = NumericalDerivatives.DifferencingMethod.Central;
        options[1] = NumericalDerivatives.DifferencingMethod.Central;
        
        // Set the increment used at the default value.
        scale[1] = 8.0e3;
        
        NumericalDerivativesEx4 derv = 
            new NumericalDerivativesEx4(new NumericalDerivativesFcn());
        derv.SetDifferencingMethods(options);
        derv.SetScalingFactors(scale);
        double[,] jacobian = derv.EvaluateJ(y);

        PrintMatrixFormat pmf = new PrintMatrixFormat();
        pmf.NumberFormat = "0.00";
        new PrintMatrix("Numerical gradient:").Print(pmf, jacobian);
        new PrintMatrix("Analytic gradient:").Print(pmf, actual);

        // Since the function is never evaluated at the
        // initial point, hold back until the request is made.
        
        // Check the relative accuracy of central differences.
        // They should be good to about two thirds-precision.
        jacobian[0, 0] = (jacobian[0, 0] - actual[0, 0]) / actual[0, 0];
        jacobian[0, 1] = (jacobian[0, 1] - actual[0, 1]) / actual[0, 1];
        re[0] = jacobian[0, 0];
        re[1] = jacobian[0, 1];
        
        Console.Out.WriteLine("Relative accuracy:");
        Console.Out.WriteLine("df/dy_1       df/dy_2");
        Console.Out.WriteLine(" {0:F2}v        {1:F2}v", 
            re[0]/v, re[1]/v);
        Console.Out.WriteLine("({0:#.###e+00})  ({1:#.###e+00})", 
            re[0], re[1]);
    }
}

Output

   Numerical gradient:
         0           1    
0  10722141354.39  60.48  

    Analytic gradient:
         0           1    
0  10722141353.42  60.48  

Relative accuracy:
df/dy_1       df/dy_2
 1.19v        0.27v
(9.1e-11)  (2.045e-11)
Link to C# source.