.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_example/plot_gill_murray_saunders_wright_benchmark.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_example_plot_gill_murray_saunders_wright_benchmark.py: Benchmark Gill, Murray, Saunders and Wright method ================================================== The goal of this example is to benchmark the :class:`~numericalderivative.GillMurraySaundersWright` on a collection of test problems. These problems are created by the :meth:`~numericalderivative.build_benchmark()` static method, which returns a list of problems. References ---------- - Gill, P. E., Murray, W., Saunders, M. A., & Wright, M. H. (1983). Computing forward-difference intervals for numerical optimization. SIAM Journal on Scientific and Statistical Computing, 4(2), 310-321. .. GENERATED FROM PYTHON SOURCE LINES 19-24 .. code-block:: Python import numpy as np import pylab as pl import tabulate import numericalderivative as nd .. GENERATED FROM PYTHON SOURCE LINES 25-27 The next function is an oracle which returns the absolute precision of the value of the function. .. GENERATED FROM PYTHON SOURCE LINES 30-116 .. code-block:: Python def absolute_precision_oracle(function, x, relative_precision): """ Return the absolute precision of the function value This oracle may fail if the function value is zero. Parameters ---------- function : function The function x : float The test point relative_precision : float, > 0, optional The relative precision of evaluation of f. Returns ------- absolute_precision : float, >= 0 The absolute precision """ function_value = function(x) if function_value == 0.0: raise ValueError( "The function value is zero: " "cannot compute the absolute precision" ) absolute_precision = relative_precision * abs(function_value) return absolute_precision class GillMurraySaundersWrightMethod: def __init__(self, kmin, kmax, relative_precision): """ Create a GillMurraySaundersWright method to compute the approximate first derivative Notice that the algorithm is parametrized here based on the relative precision of the value of the function f. Then an oracle computes the absolute precision depending on the function, the point x and the relative precision. Parameters ---------- kmin : float, kmin > 0 A minimum bound for the finite difference step of the third derivative. If no value is provided, the default is to compute the smallest possible kmin using number_of_digits and x. kmax : float, kmax > kmin > 0 A maximum bound for the finite difference step of the third derivative. If no value is provided, the default is to compute the largest possible kmax using number_of_digits and x. relative_precision : float, > 0, optional The relative precision of evaluation of f. """ self.kmin = kmin self.kmax = kmax self.relative_precision = relative_precision def compute_first_derivative(self, function, x): """ Compute the first derivative using GillMurraySaundersWright Parameters ---------- function : function The function x : float The test point Returns ------- f_prime_approx : float The approximate value of the first derivative of the function at point x number_of_function_evaluations : int The number of function evaluations. """ absolute_precision = absolute_precision_oracle( function, x, self.relative_precision ) algorithm = nd.GillMurraySaundersWright( function, x, absolute_precision=absolute_precision ) step, _ = algorithm.find_step(kmin, kmax) f_prime_approx = algorithm.compute_first_derivative(step) number_of_function_evaluations = algorithm.get_number_of_function_evaluations() return f_prime_approx, number_of_function_evaluations .. GENERATED FROM PYTHON SOURCE LINES 117-119 The next example computes the approximate derivative on the :class:`~numericalderivative.ExponentialProblem`. .. GENERATED FROM PYTHON SOURCE LINES 121-143 .. code-block:: Python print("+ Benchmark on several points") number_of_test_points = 20 kmin = 1.0e-16 kmax = 1.0e-1 problem = nd.ExponentialProblem() print(problem) interval = problem.get_interval() test_points = np.linspace(interval[0], interval[1], number_of_test_points) relative_precision = 1.0e-15 method = GillMurraySaundersWrightMethod(kmin, kmax, relative_precision) average_relative_error, average_feval, data = nd.benchmark_method( problem.get_function(), problem.get_first_derivative(), test_points, method.compute_first_derivative, True, ) print("Average relative error =", average_relative_error) print("Average number of function evaluations =", average_feval) tabulate.tabulate(data, headers=["x", "Rel. err.", "F. Eval."], tablefmt="html") .. rst-class:: sphx-glr-script-out .. code-block:: none + Benchmark on several points DerivativeBenchmarkProblem name = exp x = 1.0 f(x) = 2.718281828459045 f'(x) = 2.718281828459045 f''(x) = 2.718281828459045 f^(3)(x) = 2.718281828459045 f^(4)(x) = 2.718281828459045 f^(5)(x) = 2.718281828459045 x = 0.000, abs. error = 3.252e-08, rel. error = 3.252e-08, Func. eval. = 16 x = 0.632, abs. error = 5.986e-08, rel. error = 3.183e-08, Func. eval. = 16 x = 1.263, abs. error = 1.075e-07, rel. error = 3.040e-08, Func. eval. = 16 x = 1.895, abs. error = 2.074e-07, rel. error = 3.118e-08, Func. eval. = 16 x = 2.526, abs. error = 3.938e-07, rel. error = 3.148e-08, Func. eval. = 16 x = 3.158, abs. error = 7.683e-07, rel. error = 3.267e-08, Func. eval. = 16 x = 3.789, abs. error = 1.399e-06, rel. error = 3.163e-08, Func. eval. = 16 x = 4.421, abs. error = 2.819e-06, rel. error = 3.389e-08, Func. eval. = 16 x = 5.053, abs. error = 4.939e-06, rel. error = 3.157e-08, Func. eval. = 16 x = 5.684, abs. error = 9.331e-06, rel. error = 3.172e-08, Func. eval. = 16 x = 6.316, abs. error = 1.710e-05, rel. error = 3.090e-08, Func. eval. = 16 x = 6.947, abs. error = 3.371e-05, rel. error = 3.240e-08, Func. eval. = 16 x = 7.579, abs. error = 6.325e-05, rel. error = 3.233e-08, Func. eval. = 16 x = 8.211, abs. error = 1.126e-04, rel. error = 3.060e-08, Func. eval. = 16 x = 8.842, abs. error = 2.160e-04, rel. error = 3.122e-08, Func. eval. = 16 x = 9.474, abs. error = 4.027e-04, rel. error = 3.095e-08, Func. eval. = 16 x = 10.105, abs. error = 7.612e-04, rel. error = 3.111e-08, Func. eval. = 16 x = 10.737, abs. error = 1.391e-03, rel. error = 3.023e-08, Func. eval. = 16 x = 11.368, abs. error = 2.736e-03, rel. error = 3.161e-08, Func. eval. = 16 x = 12.000, abs. error = 5.389e-03, rel. error = 3.311e-08, Func. eval. = 16 Average relative error = 3.1666900098067314e-08 Average number of function evaluations = 16.0 .. raw:: html
x Rel. err. F. Eval.
0 3.25181e-08 16
0.6315793.18296e-08 16
1.26316 3.03995e-08 16
1.89474 3.11816e-08 16
2.52632 3.14818e-08 16
3.15789 3.26655e-08 16
3.78947 3.16328e-08 16
4.42105 3.38876e-08 16
5.05263 3.15726e-08 16
5.68421 3.17182e-08 16
6.31579 3.08999e-08 16
6.94737 3.24025e-08 16
7.57895 3.23281e-08 16
8.21053 3.06016e-08 16
8.84211 3.1215e-08 16
9.47368 3.09499e-08 16
10.1053 3.11055e-08 16
10.7368 3.02255e-08 16
11.3684 3.16141e-08 16
12 3.31086e-08 16


.. GENERATED FROM PYTHON SOURCE LINES 144-145 Map from the problem name to kmax .. GENERATED FROM PYTHON SOURCE LINES 147-166 .. code-block:: Python kmax_map = { "polynomial": 1.0, "inverse": 1.0e0, "exp": 1.0e-1, "log": 1.0e-3, # x > 0 "sqrt": 1.0e-3, # x > 0 "atan": 1.0e0, "sin": 1.0e0, "scaled exp": 1.0e5, "GMSW": 1.0e0, "SXXN1": 1.0e0, "SXXN2": 1.0e0, # Fails "SXXN3": 1.0e0, "SXXN4": 1.0e0, "Oliver1": 1.0e0, "Oliver2": 1.0e0, "Oliver3": 1.0e-3, } .. GENERATED FROM PYTHON SOURCE LINES 167-169 Benchmark the :class:`~numericalderivative.GillMurraySaundersWright` class on a collection of problems. .. GENERATED FROM PYTHON SOURCE LINES 171-227 .. code-block:: Python number_of_test_points = 100 # This value can significantly change the results relative_precision = 1.0e-15 delta_x = 1.0e-9 data = [] function_list = nd.build_benchmark() number_of_functions = len(function_list) average_relative_error_list = [] average_feval_list = [] for i in range(number_of_functions): problem = function_list[i] function = problem.get_function() first_derivative = problem.get_first_derivative() name = problem.get_name() kmax = kmax_map[name] kmin = 1.0e-16 * kmax lower_x_bound, upper_x_bound = problem.get_interval() if name == "sin": # Change the lower and upper bound so that the points +/-pi # are excluded (see below for details). lower_x_bound += delta_x upper_x_bound -= delta_x test_points = np.linspace(lower_x_bound, upper_x_bound, number_of_test_points) print(f"Function #{i}, {name}") method = GillMurraySaundersWrightMethod(kmin, kmax, relative_precision) average_relative_error, average_feval, _ = nd.benchmark_method( function, first_derivative, test_points, method.compute_first_derivative, ) average_relative_error_list.append(average_relative_error) average_feval_list.append(average_feval) data.append( ( name, kmin, kmax, average_relative_error, average_feval, ) ) data.append( [ "Average", "-", "-", np.nanmean(average_relative_error_list), np.nanmean(average_feval_list), ] ) tabulate.tabulate( data, headers=["Name", "kmin", "kmax", "Average error", "Average func. eval"], tablefmt="html", ) .. rst-class:: sphx-glr-script-out .. code-block:: none Function #0, polynomial Function #1, inverse Function #2, exp Function #3, log Function #4, sqrt Function #5, atan Function #6, sin Function #7, scaled exp Function #8, GMSW Function #9, SXXN1 Function #10, SXXN2 Function #11, SXXN3 Function #12, SXXN4 Function #13, Oliver1 Function #14, Oliver2 Function #15, Oliver3 .. raw:: html
Name kmin kmax Average error Average func. eval
polynomial1e-16 1.0 2.23087e-08 16.96
inverse 1e-16 1.0 4.48797e-08 16.88
exp 1e-17 0.1 3.16536e-08 18
log 1e-19 0.001 3.94827e-08 16.34
sqrt 1e-19 0.001 3.18116e-08 16.38
atan 1e-16 1.0 1.22055e-07 17.48
sin 1e-16 1.0 1.09185e-07 16
scaled exp1e-11 100000.0 3.16955e-08 18
GMSW 1e-16 1.0 3.10702e-08 18.28
SXXN1 1e-16 1.0 8.20049e-07 17.22
SXXN2 1e-16 1.0 5.08006e-08 12
SXXN3 1e-16 1.0 2.93358e-08 16.12
SXXN4 1e-16 1.0 2.597e-08 16.2
Oliver1 1e-16 1.0 3.14977e-08 18
Oliver2 1e-16 1.0 5.77814e-08 16.78
Oliver3 1e-19 0.001 3.47391e-08 17.08
Average - - 9.46447e-08 16.7325


.. GENERATED FROM PYTHON SOURCE LINES 228-238 Notice that the method cannot perform correctly for the sin function at the point Indeed, this function is such that :math:`f''(x) = 0` if :math:`x = \pm \pi`. In this case, the condition error is infinite and the method cannot work. Therefore, we make so that the points :math:`\pm \pi` are excluded from the benchmark. The same problem appears at the point :math:`x = 0`. This point is not included in the test set if the number of points is even (e.g. with `number_of_test_points = 100`), but it might appear if the number of test points is odd. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.065 seconds) .. _sphx_glr_download_auto_example_plot_gill_murray_saunders_wright_benchmark.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_gill_murray_saunders_wright_benchmark.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_gill_murray_saunders_wright_benchmark.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_gill_murray_saunders_wright_benchmark.zip `