From 81d22c3b565ae3d78bd576d57717294573a75f6b Mon Sep 17 00:00:00 2001 From: Michael BAUDIN Date: Tue, 10 Dec 2024 23:07:28 +0100 Subject: [PATCH] Access the threshold in the algorithms. Clarify the examples. Plot the condition error in Gill, Murray, Saunders and Wright --- README.md | 4 + doc/examples/plot_general_fd.py | 1 + .../plot_gill_murray_saunders_wright.py | 197 +++++++++++++----- doc/examples/plot_stepleman_winarsky.py | 62 ++++-- doc/examples/plot_stepleman_winarsky_plots.py | 90 +++++--- numericalderivative/_DumontetVignes.py | 130 +++++++++++- .../_GillMurraySaundersWright.py | 78 ++++++- numericalderivative/_SteplemanWinarsky.py | 37 ++++ 8 files changed, 480 insertions(+), 119 deletions(-) diff --git a/README.md b/README.md index afd218f..64577a1 100644 --- a/README.md +++ b/README.md @@ -68,6 +68,10 @@ pip install numericalderivative RAIRO. Analyse numérique, 11 (1), 13-25. ## Roadmap +- Evaluate and plot the criteria c(hphi) for Gill, Murray, Saunders & Wright. +- Evaluate and plot the criteria ell(h) for Dumontet & Vignes. +- See why the Dumontet & Vignes method fails on the polynomial function. +- See why the Gill, Murray, Saunders & Wright method fails on the sin function. - Implement the method of: Shi, H. J. M., Xie, Y., Xuan, M. Q., & Nocedal, J. (2022). Adaptive finite-difference interval estimation for noisy derivative-free optimization. _SIAM Journal on Scientific Computing_, _44_(4), A2302-A2321. diff --git a/doc/examples/plot_general_fd.py b/doc/examples/plot_general_fd.py index 4ee5ed6..d6ab0d7 100644 --- a/doc/examples/plot_general_fd.py +++ b/doc/examples/plot_general_fd.py @@ -23,6 +23,7 @@ # Compute the first derivative using forward F.D. formula # ------------------------------------------------------- + # %% # This is the function we want to compute the derivative of. def scaled_exp(x): diff --git a/doc/examples/plot_gill_murray_saunders_wright.py b/doc/examples/plot_gill_murray_saunders_wright.py index 211c348..8374cf1 100644 --- a/doc/examples/plot_gill_murray_saunders_wright.py +++ b/doc/examples/plot_gill_murray_saunders_wright.py @@ -21,9 +21,9 @@ # %% def compute_first_derivative_GMSW( - f, + function, x, - f_prime, + first_derivative, kmin, kmax, verbose=False, @@ -33,11 +33,11 @@ def compute_first_derivative_GMSW( Parameters ---------- - f : function + function : function The function. x : float The point where the derivative is to be evaluated - f_prime : function + first_derivative : function The exact first derivative of the function. kmin : float, > 0 The minimum step k for the second derivative. @@ -55,11 +55,11 @@ def compute_first_derivative_GMSW( feval : int The number of function evaluations. """ - algorithm = nd.GillMurraySaundersWright(f, x, verbose=verbose) + algorithm = nd.GillMurraySaundersWright(function, x, verbose=verbose) step, number_of_iterations = algorithm.compute_step(kmin, kmax) f_prime_approx = algorithm.compute_first_derivative(step) feval = algorithm.get_number_of_function_evaluations() - f_prime_exact = f_prime(x) + f_prime_exact = first_derivative(x) if verbose: print(f"Computed step = {step:.3e}") print(f"Number of iterations = {number_of_iterations}") @@ -74,8 +74,8 @@ def compute_first_derivative_GMSW( kmin = 1.0e-15 kmax = 1.0e1 x = 1.0 -benchmark = nd.ExponentialProblem() -second_derivative_value = benchmark.second_derivative(x) +problem = nd.ExponentialProblem() +second_derivative_value = problem.second_derivative(x) optimal_step, absolute_error = nd.FirstDerivativeForward.compute_step( second_derivative_value ) @@ -84,9 +84,9 @@ def compute_first_derivative_GMSW( absolute_error, number_of_function_evaluations, ) = compute_first_derivative_GMSW( - benchmark.function, + problem.function, x, - benchmark.first_derivative, + problem.first_derivative, kmin, kmax, verbose=True, @@ -101,8 +101,9 @@ def compute_first_derivative_GMSW( kmin = 1.0e-9 kmax = 1.0e8 x = 1.0 -benchmark = nd.ScaledExponentialProblem() -second_derivative_value = benchmark.second_derivative(x) +problem = nd.ScaledExponentialProblem() +second_derivative = problem.get_second_derivative() +second_derivative_value = second_derivative(x) optimal_step, absolute_error = nd.FirstDerivativeForward.compute_step( second_derivative_value ) @@ -111,9 +112,9 @@ def compute_first_derivative_GMSW( absolute_error, number_of_function_evaluations, ) = compute_first_derivative_GMSW( - benchmark.function, + problem.get_function(), x, - benchmark.first_derivative, + problem.get_first_derivative(), kmin, kmax, verbose=True, @@ -126,7 +127,7 @@ def compute_first_derivative_GMSW( # %% def benchmark_method( - function, derivative_function, test_points, kmin, kmax, verbose=False + function, first_derivative, test_points, kmin, kmax, verbose=False ): """ Apply Gill, Murray, Saunders & Wright method to compute the approximate first @@ -136,10 +137,10 @@ def benchmark_method( ---------- f : function The function. - derivative_function : function + first_derivative : function The exact first derivative of the function test_points : list(float) - The list of x points where the benchmark must be performed. + The list of x points where the problem must be performed. kmin : float, > 0 The minimum step k for the second derivative. kmax : float, > kmin @@ -166,9 +167,9 @@ def benchmark_method( absolute_error, number_of_function_evaluations, ) = compute_first_derivative_GMSW( - function, x, derivative_function, kmin, kmax, verbose + function, x, first_derivative, kmin, kmax, verbose ) - relative_error = absolute_error / abs(derivative_function(x)) + relative_error = absolute_error / abs(first_derivative(x)) if verbose: print( "x = %.3f, abs. error = %.3e, rel. error = %.3e, Func. eval. = %d" @@ -188,30 +189,98 @@ def benchmark_method( # %% print("+ Benchmark on several points") number_of_test_points = 100 -test_points = np.linspace(0.01, 12.2, number_of_test_points) -kmin = 1.0e-13 -kmax = 1.0e-1 -benchmark = nd.ExponentialProblem() +problem = nd.ExponentialProblem() +interval = problem.get_interval() +test_points = np.linspace(interval[0], interval[1], number_of_test_points) +kmin = 1.0e-12 +kmax = 1.0e1 average_error, average_feval = benchmark_method( - benchmark.function, benchmark.first_derivative, test_points, kmin, kmax, True + problem.get_function(), + problem.get_fifth_derivative(), + test_points, + kmin, + kmax, + True, ) +# %% +# Plot the condition error depending on k. +number_of_points = 1000 +problem = nd.ScaledExponentialProblem() +x = problem.get_x() +name = problem.get_name() +function = problem.get_function() +kmin = 1.0e-10 +kmax = 1.0e5 +k_array = np.logspace(np.log10(kmin), np.log10(kmax), number_of_points) +algorithm = nd.GillMurraySaundersWright(function, x) +c_min, c_max = algorithm.get_threshold_min_max() +condition_array = np.zeros((number_of_points)) +for i in range(number_of_points): + condition_array[i] = algorithm.compute_condition(k_array[i]) +# %% +# The next plot presents the condition error :math:`c(h_\Phi)` depending +# on :math:`h_\Phi`. +# The two horizontal lines represent the minimum and maximum threshold +# values. +# We search for the value of :math:`h_\Phi` such that the condition +# error is between these two limits. + +# %% +pl.figure(figsize=(4.0, 3.0)) +pl.title(f"Condition error of {name} at x = {x}") +pl.plot(k_array, condition_array) +pl.plot([kmin, kmax], [c_min] * 2, label=r"$c_{min}$") +pl.plot([kmin, kmax], [c_max] * 2, label=r"$c_{max}$") +pl.xlabel(r"$h_\Phi$") +pl.ylabel(r"$c(h_\Phi$)") +pl.xscale("log") +pl.yscale("log") +pl.legend(bbox_to_anchor=(1.0, 1.0)) +pl.tight_layout() + +# %% +# The previous plot shows that the condition error is a decreasing function +# of :math:`h_\Phi`. + +# %% # For each function, at point x = 1, plot the error vs the step computed # by the method # %% def plot_error_vs_h_with_GMSW_steps( - name, function, function_derivative, x, h_array, kmin, kmax, verbose=False + name, function, first_derivative, x, step_array, kmin, kmax, verbose=False ): + """ + Plot the computed error depending on the step for an array of F.D. steps + + Parameters + ---------- + name : str + The name of the problem + function : function + The function. + first_derivative : function + The exact first derivative of the function + x : float + The input point where the test is done + step_array : list(float) + The array of finite difference steps + kmin : float, > 0 + The minimum step k for the second derivative. + kmax : float, > kmin + The maximum step k for the second derivative. + verbose : bool, optional + Set to True to print intermediate messages. The default is False. + """ algorithm = nd.GillMurraySaundersWright(function, x) - number_of_points = len(h_array) + number_of_points = len(step_array) error_array = np.zeros((number_of_points)) for i in range(number_of_points): - h = h_array[i] - f_prime_approx = algorithm.compute_first_derivative(h) - error_array[i] = abs(f_prime_approx - function_derivative(x)) + f_prime_approx = algorithm.compute_first_derivative(step_array[i]) + error_array[i] = abs(f_prime_approx - first_derivative(x)) step, number_of_iterations = algorithm.compute_step(kmin, kmax) @@ -223,7 +292,7 @@ def plot_error_vs_h_with_GMSW_steps( maximum_error = np.nanmax(error_array) pl.figure() - pl.plot(h_array, error_array) + pl.plot(step_array, error_array) pl.plot( [step] * 2, [minimum_error, maximum_error], @@ -241,13 +310,31 @@ def plot_error_vs_h_with_GMSW_steps( # %% -def plot_error_vs_h_benchmark(benchmark, x, h_array, kmin, kmax, verbose=False): +def plot_error_vs_h_benchmark(problem, x, step_array, kmin, kmax, verbose=False): + """ + Plot the computed error depending on the step for an array of F.D. steps + + Parameters + ---------- + problem : nd.BenchmarkProblem + The problem + x : float + The input point where the test is done + step_array : list(float) + The array of finite difference steps + kmin : float, > 0 + The minimum step k for the second derivative. + kmax : float, > kmin + The maximum step k for the second derivative. + verbose : bool, optional + Set to True to print intermediate messages. The default is False. + """ plot_error_vs_h_with_GMSW_steps( - benchmark.name, - benchmark.function, - benchmark.first_derivative, + problem.get_name(), + problem.get_function(), + problem.get_first_derivative(), x, - h_array, + step_array, kmin, kmax, verbose, @@ -255,57 +342,57 @@ def plot_error_vs_h_benchmark(benchmark, x, h_array, kmin, kmax, verbose=False): # %% -benchmark = nd.ExponentialProblem() +problem = nd.ExponentialProblem() x = 1.0 number_of_points = 1000 -h_array = np.logspace(-15.0, -1.0, number_of_points) +step_array = np.logspace(-15.0, -1.0, number_of_points) kmin = 1.0e-15 kmax = 1.0e-1 -plot_error_vs_h_benchmark(benchmark, x, h_array, kmin, kmax, True) +plot_error_vs_h_benchmark(problem, x, step_array, kmin, kmax, True) # %% x = 12.0 -h_array = np.logspace(-15.0, -1.0, number_of_points) -plot_error_vs_h_benchmark(benchmark, x, h_array, kmin, kmax) +step_array = np.logspace(-15.0, -1.0, number_of_points) +plot_error_vs_h_benchmark(problem, x, step_array, kmin, kmax) # %% -benchmark = nd.ScaledExponentialProblem() +problem = nd.ScaledExponentialProblem() x = 1.0 kmin = 1.0e-10 kmax = 1.0e8 -h_array = np.logspace(-10.0, 8.0, number_of_points) -plot_error_vs_h_benchmark(benchmark, x, h_array, kmin, kmax) +step_array = np.logspace(-10.0, 8.0, number_of_points) +plot_error_vs_h_benchmark(problem, x, step_array, kmin, kmax) # %% -benchmark = nd.LogarithmicProblem() +problem = nd.LogarithmicProblem() x = 1.1 kmin = 1.0e-14 kmax = 1.0e-1 -h_array = np.logspace(-15.0, -1.0, number_of_points) -plot_error_vs_h_benchmark(benchmark, x, h_array, kmin, kmax, True) +step_array = np.logspace(-15.0, -1.0, number_of_points) +plot_error_vs_h_benchmark(problem, x, step_array, kmin, kmax, True) # %% -benchmark = nd.SinProblem() +problem = nd.SinProblem() x = 1.0 kmin = 1.0e-15 kmax = 1.0e-1 -h_array = np.logspace(-15.0, -1.0, number_of_points) -plot_error_vs_h_benchmark(benchmark, x, h_array, kmin, kmax) +step_array = np.logspace(-15.0, -1.0, number_of_points) +plot_error_vs_h_benchmark(problem, x, step_array, kmin, kmax) # %% -benchmark = nd.SquareRootProblem() +problem = nd.SquareRootProblem() x = 1.0 kmin = 1.0e-15 kmax = 1.0e-1 -h_array = np.logspace(-15.0, -1.0, number_of_points) -plot_error_vs_h_benchmark(benchmark, x, h_array, kmin, kmax, True) +step_array = np.logspace(-15.0, -1.0, number_of_points) +plot_error_vs_h_benchmark(problem, x, step_array, kmin, kmax, True) # %% -benchmark = nd.AtanProblem() +problem = nd.AtanProblem() x = 1.1 kmin = 1.0e-15 kmax = 1.0e-1 -h_array = np.logspace(-15.0, -1.0, number_of_points) -plot_error_vs_h_benchmark(benchmark, x, h_array, kmin, kmax) +step_array = np.logspace(-15.0, -1.0, number_of_points) +plot_error_vs_h_benchmark(problem, x, step_array, kmin, kmax) # %% diff --git a/doc/examples/plot_stepleman_winarsky.py b/doc/examples/plot_stepleman_winarsky.py index ee20f6a..28a4e28 100644 --- a/doc/examples/plot_stepleman_winarsky.py +++ b/doc/examples/plot_stepleman_winarsky.py @@ -41,16 +41,16 @@ x = 1.0 finite_difference = nd.FirstDerivativeCentral(benchmark.function, x) number_of_points = 1000 -h_array = np.logspace(-7.0, 5.0, number_of_points) +step_array = np.logspace(-7.0, 5.0, number_of_points) error_array = np.zeros((number_of_points)) for i in range(number_of_points): - h = h_array[i] + h = step_array[i] f_prime_approx = finite_difference.compute(h) error_array[i] = abs(f_prime_approx - benchmark.first_derivative(x)) # %% pl.figure() -pl.plot(h_array, error_array) +pl.plot(step_array, error_array) pl.plot([optimum_step] * 2, [min(error_array), max(error_array)], label=r"$h^*$") pl.title("Central finite difference") pl.xlabel("h") @@ -77,8 +77,30 @@ # %% -def fd_difference(h1, h2, f, x): - finite_difference = nd.FirstDerivativeCentral(f, x) +def fd_difference(h1, h2, function, x): + """ + Compute the difference of central difference approx. for different step sizes + + This function computes the absolute value of the difference of approximations + evaluated at two different steps h1 and h2: + + d = abs(FD(h1) - FD(h2)) + + where FD(h) is the approximation from the finite difference formula + evaluated from the step h. + + Parameters + ---------- + h1 : float, > 0 + The first step + h2 : float, > 0 + The second step + function : function + The function + x : float + The input point where the derivative is approximated. + """ + finite_difference = nd.FirstDerivativeCentral(function, x) f_prime_approx_1 = finite_difference.compute(h1) f_prime_approx_2 = finite_difference.compute(h2) diff_current = abs(f_prime_approx_1 - f_prime_approx_2) @@ -88,15 +110,15 @@ def fd_difference(h1, h2, f, x): # %% # 3. Plot the evolution of | FD(h) - FD(h / 2) | for different values of h number_of_points = 1000 -h_array = np.logspace(-7.0, 5.0, number_of_points) +step_array = np.logspace(-7.0, 5.0, number_of_points) diff_array = np.zeros((number_of_points)) for i in range(number_of_points): - h = h_array[i] + h = step_array[i] diff_array[i] = fd_difference(h, h / 2, benchmark.function, x) # %% pl.figure() -pl.plot(h_array, diff_array) +pl.plot(step_array, diff_array) pl.title("F.D. difference") pl.xlabel("h") pl.ylabel(r"$|\operatorname{FD}(h) - \operatorname{FD}(h / 2) |$") @@ -110,19 +132,21 @@ def fd_difference(h1, h2, f, x): number_of_points = 20 h_initial = 1.0e5 beta = 4.0 -h_array = np.zeros((number_of_points)) +step_array = np.zeros((number_of_points)) diff_array = np.zeros((number_of_points)) for i in range(number_of_points): if i == 0: - h_array[i] = h_initial / beta - diff_array[i] = fd_difference(h_array[i], h_initial, benchmark.function, x) + step_array[i] = h_initial / beta + diff_array[i] = fd_difference(step_array[i], h_initial, benchmark.function, x) else: - h_array[i] = h_array[i - 1] / beta - diff_array[i] = fd_difference(h_array[i], h_array[i - 1], benchmark.function, x) + step_array[i] = step_array[i - 1] / beta + diff_array[i] = fd_difference( + step_array[i], step_array[i - 1], benchmark.function, x + ) # %% pl.figure() -pl.plot(h_array, diff_array, "o") +pl.plot(step_array, diff_array, "o") pl.title("F.D. difference") pl.xlabel("h") pl.ylabel(r"$|\operatorname{FD}(h) - \operatorname{FD}(h / 2) |$") @@ -159,21 +183,21 @@ def fd_difference(h1, h2, f, x): # %% number_of_points = 1000 -h_array = np.logspace(-7.0, 7.0, number_of_points) +step_array = np.logspace(-7.0, 7.0, number_of_points) n_digits_array = np.zeros((number_of_points)) for i in range(number_of_points): - h = h_array[i] + h = step_array[i] n_digits_array[i] = algorithm.number_of_lost_digits(h) # %% y_max = algorithm.number_of_lost_digits(h_reference) pl.figure() -pl.plot(h_array, n_digits_array, label="$N(h)$") +pl.plot(step_array, n_digits_array, label="$N(h)$") pl.plot([h_reference] * 2, [0.0, y_max], "--", label=r"$h_{ref}$") pl.plot([step_zero] * 2, [0.0, y_max], "--", label=r"$h^{(0)}$") pl.plot([estim_step] * 2, [0.0, y_max], "--", label=r"$h^\star$") pl.plot( - h_array, + step_array, np.array([threshold] * number_of_points), ":", label=r"$T$", @@ -188,7 +212,7 @@ def fd_difference(h1, h2, f, x): # %% pl.figure() -pl.plot(h_array, error_array) +pl.plot(step_array, error_array) pl.plot([step_zero] * 2, [0.0, 1.0e-9], "--", label=r"$h^{(0)}$") pl.plot([estim_step] * 2, [0.0, 1.0e-9], "--", label=r"$h^\star$") pl.title("Finite difference") diff --git a/doc/examples/plot_stepleman_winarsky_plots.py b/doc/examples/plot_stepleman_winarsky_plots.py index 7ec4ef8..143663b 100644 --- a/doc/examples/plot_stepleman_winarsky_plots.py +++ b/doc/examples/plot_stepleman_winarsky_plots.py @@ -26,15 +26,15 @@ # Plot the number of lost digits for exp number_of_points = 100 x = 1.0 -h_array = np.logspace(-15.0, 1.0, number_of_points) +step_array = np.logspace(-15.0, 1.0, number_of_points) n_digits_array = np.zeros((number_of_points)) algorithm = nd.SteplemanWinarsky(np.exp, x) for i in range(number_of_points): - h = h_array[i] + h = step_array[i] n_digits_array[i] = algorithm.number_of_lost_digits(h) pl.figure() -pl.plot(h_array, n_digits_array) +pl.plot(step_array, n_digits_array) pl.title(r"Number of digits lost by F.D.. $f(x) = \exp(x)$") pl.xlabel("h") pl.ylabel("$N(h)$") @@ -43,16 +43,16 @@ # %% # Plot the number of lost digits for sin x = 1.0 -h_array = np.logspace(-7.0, 2.0, number_of_points) +step_array = np.logspace(-7.0, 2.0, number_of_points) n_digits_array = np.zeros((number_of_points)) algorithm = nd.SteplemanWinarsky(np.sin, x) for i in range(number_of_points): - h = h_array[i] + h = step_array[i] n_digits_array[i] = algorithm.number_of_lost_digits(h) # %% pl.figure() -pl.plot(h_array, n_digits_array) +pl.plot(step_array, n_digits_array) pl.title(r"Number of digits lost by F.D.. $f(x) = \sin(x)$") pl.xlabel("h") pl.ylabel("$N(h)$") @@ -65,13 +65,35 @@ # %% def plot_error_vs_h_with_SW_steps( - name, function, function_derivative, x, h_array, h_min, h_max, verbose=False + name, function, function_derivative, x, step_array, h_min, h_max, verbose=False ): + """ + Plot the computed error depending on the step for an array of F.D. steps + + Parameters + ---------- + name : str + The name of the problem + function : function + The function. + first_derivative : function + The exact first derivative of the function + x : float + The input point where the test is done + step_array : list(float) + The array of finite difference steps + h_min : float, > 0 + The lower bound to bracket the initial differentiation step. + h_max : float, > kmin + The upper bound to bracket the initial differentiation step. + verbose : bool, optional + Set to True to print intermediate messages. The default is False. + """ algorithm = nd.SteplemanWinarsky(function, x) - number_of_points = len(h_array) + number_of_points = len(step_array) error_array = np.zeros((number_of_points)) for i in range(number_of_points): - h = h_array[i] + h = step_array[i] f_prime_approx = algorithm.compute_first_derivative(h) error_array[i] = abs(f_prime_approx - function_derivative(x)) @@ -93,7 +115,7 @@ def plot_error_vs_h_with_SW_steps( maximum_error = np.nanmax(error_array) pl.figure() - pl.plot(h_array, error_array) + pl.plot(step_array, error_array) pl.plot( [h_min] * 2, [minimum_error, maximum_error], @@ -124,13 +146,31 @@ def plot_error_vs_h_with_SW_steps( # %% -def plot_error_vs_h_benchmark(benchmark, x, h_array, h_min, h_max, verbose=False): +def plot_error_vs_h_benchmark(benchmark, x, step_array, h_min, h_max, verbose=False): + """ + Plot the computed error depending on the step for an array of F.D. steps + + Parameters + ---------- + problem : nd.BenchmarkProblem + The problem + x : float + The input point where the test is done + step_array : list(float) + The array of finite difference steps + kmin : float, > 0 + The minimum step k for the second derivative. + kmax : float, > kmin + The maximum step k for the second derivative. + verbose : bool, optional + Set to True to print intermediate messages. The default is False. + """ plot_error_vs_h_with_SW_steps( benchmark.name, benchmark.function, benchmark.first_derivative, x, - h_array, + step_array, h_min, h_max, True, @@ -141,42 +181,42 @@ def plot_error_vs_h_benchmark(benchmark, x, h_array, h_min, h_max, verbose=False benchmark = nd.ExponentialProblem() x = 1.0 number_of_points = 1000 -h_array = np.logspace(-15.0, 1.0, number_of_points) -plot_error_vs_h_benchmark(benchmark, x, h_array, 1.0e-10, 1.0e0, True) +step_array = np.logspace(-15.0, 1.0, number_of_points) +plot_error_vs_h_benchmark(benchmark, x, step_array, 1.0e-10, 1.0e0, True) # %% x = 12.0 -h_array = np.logspace(-15.0, 1.0, number_of_points) -plot_error_vs_h_benchmark(benchmark, x, h_array, 1.0e-10, 1.0e0) +step_array = np.logspace(-15.0, 1.0, number_of_points) +plot_error_vs_h_benchmark(benchmark, x, step_array, 1.0e-10, 1.0e0) if False: benchmark = nd.LogarithmicProblem() x = 1.0 - plot_error_vs_h_benchmark(benchmark, x, h_array, 1.0e-15, 1.0e0, True) + plot_error_vs_h_benchmark(benchmark, x, step_array, 1.0e-15, 1.0e0, True) # %% benchmark = nd.LogarithmicProblem() x = 1.1 -h_array = np.logspace(-15.0, -1.0, number_of_points) -plot_error_vs_h_benchmark(benchmark, x, h_array, 1.0e-14, 1.0e-1, True) +step_array = np.logspace(-15.0, -1.0, number_of_points) +plot_error_vs_h_benchmark(benchmark, x, step_array, 1.0e-14, 1.0e-1, True) # %% benchmark = nd.SinProblem() x = 1.0 -h_array = np.logspace(-15.0, 0.0, number_of_points) -plot_error_vs_h_benchmark(benchmark, x, h_array, 1.0e-15, 1.0e-0) +step_array = np.logspace(-15.0, 0.0, number_of_points) +plot_error_vs_h_benchmark(benchmark, x, step_array, 1.0e-15, 1.0e-0) # %% benchmark = nd.SquareRootProblem() x = 1.0 -h_array = np.logspace(-15.0, 0.0, number_of_points) -plot_error_vs_h_benchmark(benchmark, x, h_array, 1.0e-15, 1.0e-0, True) +step_array = np.logspace(-15.0, 0.0, number_of_points) +plot_error_vs_h_benchmark(benchmark, x, step_array, 1.0e-15, 1.0e-0, True) # %% benchmark = nd.AtanProblem() x = 1.0 -h_array = np.logspace(-15.0, 0.0, number_of_points) -plot_error_vs_h_benchmark(benchmark, x, h_array, 1.0e-15, 1.0e-0) +step_array = np.logspace(-15.0, 0.0, number_of_points) +plot_error_vs_h_benchmark(benchmark, x, step_array, 1.0e-15, 1.0e-0) # %% benchmark = nd.ExponentialProblem() diff --git a/numericalderivative/_DumontetVignes.py b/numericalderivative/_DumontetVignes.py index 604b425..e65967a 100644 --- a/numericalderivative/_DumontetVignes.py +++ b/numericalderivative/_DumontetVignes.py @@ -9,13 +9,102 @@ class DumontetVignes: - """ + r""" Compute an approximately optimal step for the central F.D. formula for the first derivative The method is based on computing the third derivative. Then the optimal step for the central formula for the first derivative is computed from the third derivative. + The goal of the method is to compute the step of the central finite + difference approximation of the first derivative (see + (Dumontet & Vignes, 1977) eq. 2 page 13): + + .. math:: + + f'(x) \approx \frac{f(x + h) - f(x - h)}{2h} + + where :math:`f` is the function, :math:`x \in \mathbb{R}` is the + input point and :math:`h > 0` is the step. + The optimal step is (see (Dumontet & Vignes, 1977) eq. 24 page 18): + + .. math:: + + h^\star + = \left( \frac{3 \epsilon_f \left|f(x)\right|}{\left|f^{(3)}(x)\right|} \right)^{1/3} + + where :math:`\epsilon_f > 0` is the absolute error of the function evaluation. + The goal of the method is to compute :math:`h^\star` using + function evaluations only. + + The third derivative is approximated using the central finite difference formula + (see (Dumontet & Vignes, 1977) eq. 25 page 18): + + .. math:: + + f^{(3)}_k(x) + = \frac{f(x + 2k) - f(x - 2k) - 2 (f(x + k) - f(x - k))}{2k^3} + + where :math:`k > 0` is the step used for the third derivative. + The method introduces :math:`f^{(3)}_{inf}(x)` and :math:`f^{(3)}_{sup}(x)` + such that: + + .. math:: + + f^{(3)}_{inf}(x) \leq f^{(3)}_k(x_0) \leq f^{(3)}_{sup}(x). + + We evaluate the function on 4 points (see (Dumontet & Vignes, 1977) eq. 28 + page 19, with 2 errors corrected with respect to the original paper): + + .. math:: + + & T_1 = f(x + 2k), \qquad T_2 = -f(x - 2k), \\ + & T_3 = -2f(x + k) , \qquad T_4 = 2f(x - k). + + Let :math:`A` and :math:`B` defined by (see (Dumontet & Vignes, 1977) page 19): + + .. math:: + + A = \sum_{i = 1}^4 \max(T_i, 0), \quad + B = \sum_{i = 1}^4 \min(T_i, 0). \tag{16} + + The lower and upper bounds of :math:`f^{(3)}_k(x)` are computed + from (see (Dumontet & Vignes, 1977) eq. 30 page 20): + + .. math:: + + f^{(3)}_{inf}(x_0) + = \frac{\frac{A}{1 + \epsilon_f} + \frac{B}{1 - \epsilon_f}}{2 k^3}, \qquad + f^{(3)}_{sup}(x_0) + = \frac{\frac{A}{1 - \epsilon_f} + \frac{B}{1 + \epsilon_f}}{2 k^3}. + + We introduce the ratio (see (Dumontet & Vignes, 1977) eq. 32 page 20): + + .. math:: + + L(k) = \frac{f^{(3)}_{sup}(x)}{f^{(3)}_{inf}(x)} \geq 1. + + We search for :math:`k` such that the ratio :math:`L` is: + + - neither too close to 1 because it would mean that :math:`k` is too large + meaning that the truncation error dominates, + - nor too far away from 1 because it would mean that :math:`k` is too small + meaning that the rounding error dominates. + + Let :math:`k_{inf}` and :math:`k_{sup}` two real numbers representing the + minimum and maximum bounds for :math:`k`. + We search for :math:`k \in [k_{inf}, k_{sup}]` such that (see + (Dumontet & Vignes, 1977) eq. 33 page 20): + + .. math:: + + \ell(k) \in [L_1, L_2] \cup [L_3, L_4] + + where: + + - :math:`L_3 = 2` and :math:`L_2 = \frac{1}{L_3}`, + - :math:`L_4 = 15` and :math:`L_1 = \frac{1}{L_4}`. + Parameters ---------- function : function @@ -26,9 +115,9 @@ class DumontetVignes: The relative precision of evaluation of f. number_of_digits : int The maximum number of digits of the floating point system. - ell_1 : float + ell_3 : float The minimum bound of the L ratio. - ell_2 : float, > ell_1 + ell_4 : float, > ell_1 The maximum bound of the L ratio. args : list A list of optional arguments that the function takes as inputs. @@ -70,8 +159,8 @@ def __init__( x, relative_precision=1.0e-14, number_of_digits=53, - ell_1=1.0 / 15.0, - ell_2=1.0 / 2.0, + ell_3=2.0, + ell_4=15.0, args=None, verbose=False, ): @@ -82,20 +171,39 @@ def __init__( ) self.relative_precision = relative_precision self.number_of_digits = number_of_digits - if ell_2 <= ell_1: + if ell_4 <= ell_3: raise ValueError( - f"We must have ell_2 > ell_1, but ell_1 = {ell_1} and ell_2 = {ell_2}" + f"We must have ell_4 > ell_3, but ell_4 = {ell_4} and ell_3 = {ell_3}" ) # Eq. 34, fixed - self.ell_1 = ell_1 - self.ell_2 = ell_2 - self.ell_3 = 1.0 / ell_2 - self.ell_4 = 1.0 / ell_1 + self.ell_3 = ell_3 + self.ell_4 = ell_4 + self.ell_1 = 1.0 / ell_4 + self.ell_2 = 1.0 / ell_3 self.verbose = verbose self.first_derivative_central = nd.FirstDerivativeCentral(function, x, args) self.function = nd.FunctionWithArguments(function, args) self.x = x + def get_ell_min_max(self): + """ + Return the minimum and maximum of the L ratio + + The parameters L1 and L2 can be computed from the equations: + + .. math:: + + L_2 = \frac{1}{L_3}, \qquad L_1 = \frac{1}{L_4}. + + Returns + ------- + ell_3 : float, > 0 + The lower bound of the L ratio. + ell_4 : float, > 0 + The upper bound of the L ratio. + """ + return [self.ell_3, self.ell_4] + def compute_ell(self, k): """ Compute the L ratio depending on k. diff --git a/numericalderivative/_GillMurraySaundersWright.py b/numericalderivative/_GillMurraySaundersWright.py index 12d0c1d..828bcef 100644 --- a/numericalderivative/_GillMurraySaundersWright.py +++ b/numericalderivative/_GillMurraySaundersWright.py @@ -9,18 +9,66 @@ class GillMurraySaundersWright: - """ + r""" Compute an approximately optimal step for the forward F.D. formula of the first derivative The method is based on three steps: - - compute an approximate optimal step for the second derivative using central finite difference formula, - - compute the approximate second derivative using central finite difference formula, - - compute the approximate optimal step for the first derivative using the forward finite difference formula. + - compute an approximate optimal step :math:`h_\Phi` for the second + derivative using central finite difference formula, + - compute the approximate second derivative using central finite + difference formula, + - compute the approximate optimal step for the first derivative using the + forward finite difference formula. Finally, this approximately optimal step can be use to compute the first derivative using the forward finite difference formula. + The goal of the method is to compute the approximation of :math:`f'(x)` + using the forward finite difference formula (see (G, M, S & W, 1983) eq. 1 page 311): + + .. math:: + + f'(x) \approx \frac{f(x + h) - f(x)}{h} + + where :math:`f` is the function, :math:`x \in \mathbb{R}` is the + input point and :math:`h > 0` is the step. + If :math:`f''(x) \neq 0`, then the step which minimizes the total error is: + + .. math:: + + h^\star = 2 \sqrt{\frac{\epsilon_f}{\left|f''(x)\right|}} + + where :math:`\epsilon_f > 0` is the absolute error of the function evaluation. + The goal of the method is to compute :math:`h^\star` using + function evaluations only. + An approximate value of the second derivative can be computed from the + central finite difference formula (see (G, M, S & W, 1983) eq. 8 page 314): + + .. math:: + + f''(x) \approx \Phi(h_\Phi) + = \frac{f(x + h_\Phi) - 2 f(x) + f(x - h_\Phi)}{h_\Phi^2}. + + where :math:`\Phi` is the approximation of :math:`f''(x)` from the + central finite difference formula and :math:`h_\Phi > 0` is the step of + the second derivative finite difference formula. + The method is based on the condition error (see (G, M, S & W, 1983) eq. 1 page 315): + + .. math:: + + c(h_\Phi) = \frac{4 \epsilon_f}{h_\Phi^2 |\Phi|}. + + The condition error is a decreasing function of :math:`h_\Phi`. + The algorithm searches for a step :math:`h_\Phi` such that: + + .. math:: + + c_{\min} \leq c(h_\Phi) \leq c_{\max} + + where :math:`c_{\min}` and :math:`c_{\max}` are thresholds defined by the + user. + This algorithm is a simplified version of the algorithm published in (Gill, Murray, Saunders & Wright, 1983) section 3.2 page 316. While (Gill, Murray, Saunders & Wright, 1983) simultaneously considers @@ -110,18 +158,30 @@ def __init__( self.absolute_precision = abs(relative_precision * self.y) self.first_derivative_forward = nd.FirstDerivativeForward(function, x, args) + def get_threshold_min_max(self): + """ + Return the threshold min and max of the condition error + + Returns + ------- + c_threshold_min : float, > 0 + The minimum value of the threshold of the condition error. + c_threshold_max : float, > 0 + The maxiimum value of the threshold of the condition error. + """ + return [self.c_threshold_min, self.c_threshold_max] + def compute_condition(self, k): r""" Compute the condition error for given step k. - This is the condition error of the finite difference formula - of the second derivative finite difference : + This function evaluates the condition error :math:`c(h_\Phi)` of the + finite difference formula of the second derivative finite difference + depending on the step :math:`h_\Phi`: .. math:: - f''(x) \approx \frac{f(x + k) - 2 f(x) + f(x - k)}{k^2} - - The condition error is a decreasing function of k. + c(h_\Phi) = \frac{4 \epsilon_f}{h_\Phi^2 |\Phi|}. Parameters ---------- diff --git a/numericalderivative/_SteplemanWinarsky.py b/numericalderivative/_SteplemanWinarsky.py index d3b399b..049e1dd 100644 --- a/numericalderivative/_SteplemanWinarsky.py +++ b/numericalderivative/_SteplemanWinarsky.py @@ -21,6 +21,43 @@ class SteplemanWinarsky: f'(x) \approx \frac{f(x + h) - f(x - h)}{2 h} + where :math:`f` is the function, :math:`x \in \mathbb{R}` is the + input point and :math:`h > 0` is the step. + The goal of the method is to compute an approximately optimal + :math:`h^\star` for the central F.D. formula using function evaluations + only. + + The method introduces the function :math:`g` defined by: + + .. math:: + + g(t) = f(x + t) - f(x - t) + + for any :math:`t \in \mathbb{R}`. + We introduce the monotic sequence of step sizes :math:`\{h_i\}_{i \geq 0}` defined + by the equation: + + .. math:: + + h_{i + 1} = \frac{h_i}{\beta}, \quad i=0,1,2,... + + Therefore, under some smoothness hypotheses on :math:`g`, + there exists :math:`N \in \mathbb{N}` such that for any + :math:`i \geq N`, we have: + + .. math:: + + \left|d(h_{i + 1}) - d(h_i)\right| + \leq \left|d(h_{i}) - d(h_{i - 1})\right|. + + The previous theorem states that the sequence + :math:`\left\{\left|d(h_{i}) - d(h_{i - 1})\right|\right\}_{i \geq 0}` + is monotonic in exact arithmetic. + + The method starts from an initial step :math:`h_0 > 0`. + It then reduces this step iteratively until the monotonicity property + is broken. + Parameters ---------- function : function