Example 3: Accumulation Of A Component

This example uses the same data as in the One-Sided Differences example. An alternate examination of the function f(y_1 ,y_2 ) = a\exp (by_1 ) + cy_1 y_2^2 shows that the first term on the right-hand side need be evaluated just when computing the first partial. The additive term cy_2^2 occurs when computing the partial with respect to y_1. Also the first term does not depend on the second variable. Thus the first term can be left out of the function evaluation when computing the partial with respect to y_2, potentially avoiding cancellation errors. The input values of array options allow NumericalDerivatives to use these facts and obtain greater accuracy using a minimum number of computations of the exponential function. 
using System;
using Imsl.Math;

public class NumericalDerivativesEx3 : NumericalDerivatives.IJacobian  
{    
    static int m = 1, n = 2;
    static double a, b, c;
    

    public double[] F(int varIndex, double[] y) 
    {
        double[] tmp = new double[m];
                
        if (varIndex != 2) 
        {
            tmp[0] = a * Math.Exp(b * y[0]);
        } 
        else 
        {
            // This is the function value for the partial wrt y_2.
            tmp[0] = c * y[0] * y[1] * y[1];
        }
                
        return tmp;
    }
         
   
    public double[,] Jacobian(double[] y) 
    {
        double[,] tmp = new double[m, n];
                
        // Start with part of the derivative that is known.
        tmp[0, 0] = c * y[1] * y[1];
                
        return tmp;
    }
    

    public static void Main(String[] args) 
    {
        NumericalDerivatives.DifferencingMethod[] options = 
            new NumericalDerivatives.DifferencingMethod[n];
        double u;
        double[] y = new double[n];
        double[] valueF = new double[m];
        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;
        
        // Precision, for measuring errors
        u = Math.Sqrt(2.220446049250313e-016);
        
        // 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.Accumulate;
        options[1] = NumericalDerivatives.DifferencingMethod.OneSided;
        
        valueF[0] = a * Math.Exp(b * y[0]);
        scale[1] = 1.0e0;
        
        NumericalDerivatives derv = 
            new NumericalDerivatives(new NumericalDerivativesEx3());
        derv.SetDifferencingMethods(options);
        derv.SetScalingFactors(scale);
        derv.SetInitialF(valueF);
        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);
   
        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}u        {1:F2}u", 
            re[0]/u, re[1]/u);
        Console.Out.WriteLine("({0:#.###e+00})  ({1:#.###e+00})", 
            re[0], re[1]);
    }
}

Output

   Numerical gradient:
         0           1    
0  10722141710.08  60.48  

    Analytic gradient:
         0           1    
0  10722141353.42  60.48  

Relative accuracy:
df/dy_1       df/dy_2
 2.23u        -0.51u
(3.326e-08)  (-7.569e-09)
Link to C# source.