Example 1: One-Sided Differences

A simple use of NumericalDerivatives is shown. The gradient of the function f(y_1 ,y_2 ) = a\exp (by_1 ) + cy_1 y_2^2 is required for values a = 2.5e6, b = 3.4, c = 4.5, y_1 = 2.1, y_2 = 3.2.

The numerical gradient is compared to the analytic gradient, cast as a 1 by 2 Jacobian:

grad(f) = \left[ {\begin{array}{*{20}c} {ab\exp (by_1 ) + cy_2^2 ,} & {2cy_1 y_2 } \\ \end{array} } \right]
This analytic gradient is expected to approximately agree with the numerical differentiation gradient. Relative agreement should be approximately the square root of machine precision. That is achieved here. Generally this is the most accuracy one can expect using one-sided divided differences.

import com.imsl.math.*;
import java.text.*;

public class NumericalDerivativesEx1 {
    static int m = 1, n = 2;
    static double a, b, c;
    

    public static void main(String args[]) {
        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;
        
        // Precision for measuring errors
        u = Math.sqrt(2.220446049250313e-016);
        
        // Set scaling:
        scale[0] = 1.e0;
        // Increase scale to account for large value of a.
        scale[1] = 8.e3;
        
        // 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];
               
        // This sets the function value used in forming one-sided
        // differences.
        NumericalDerivatives.Function fcn = new NumericalDerivatives.Function() {
            public double[] f(int varIndex, double[] y) {
                double[] tmp = new double[m];    
                tmp[0] = a * Math.exp(b * y[0]) + c * y[0] * y[1] * y[1];
                return tmp;
            }
        };
        
        NumericalDerivatives derv = new NumericalDerivatives(fcn);
        derv.setScalingFactors(scale);
        double[][] jacobian = derv.evaluateJ(y);
        
        NumberFormat nf = NumberFormat.getInstance();
        nf.setMaximumFractionDigits(2);
        nf.setMinimumFractionDigits(2);
        
        PrintMatrixFormat pmf = new PrintMatrixFormat();
        pmf.setNumberFormat(nf);
        new PrintMatrix("Numerical gradient:").print(pmf, jacobian);
        new PrintMatrix("Analytic gradient:").print(pmf, actual);
        
        // Check the relative accuracy of one-sided differences.
        // They should be good to about half-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];
        
        System.out.println("Relative accuracy:");
        System.out.println("df/dy_1       df/dy_2");
        System.out.printf(" %.2fu        %.2fu\n", re[0]/u, re[1]/u);
        System.out.printf("(%.3e)  (%.3e)\n", re[0], re[1]);
    }
}

Output

     Numerical gradient:
           0            1    
0  10,722,141,696.00  60.48  

     Analytic gradient:
           0            1    
0  10,722,141,353.42  60.48  

Relative accuracy:
df/dy_1       df/dy_2
 2.14u        -0.00u
(3.195e-08)  (-1.175e-16)
Link to Java source.