Ctrl+K

Adaptation of step size

3.2. Adaptation of step size#

3.2.1. Estimation of nuemrical errors#

In the previous section, we are able to see good values of \(h\) because we know the exact derivative. In practice, we don’t know the actual error and therefore we cannot use a plot like the right plot of Example 2.1.1 to find an appropriate value for \(h\). Equation eq}eq_mean_diff suggests that the leading error is given by \(\frac{|f^{(3)}|}{3!} h^2\). Unless we know the third order derivative, we still don’t know the magnitude of error. However, we can estimate it numerically. Suppose that we calculate the derivative using two different vales of step sizes, \(h\) and \(h'=h/2\). We expect that the output is more accurate with \(h'\) than \(h\). Using Eq. (3.1), the change of the mean finite difference is given by

\[ \left | \Delta_\text{M} f(x) - \Delta'_\text{M} f(x) \right | \approx \frac{|f^{(3)}|}{3!} \frac{3}{4} h^2 = \frac{3}{4} \times \left [ \text{error in } \Delta_\text{M} f(x) \right ] \]

which suggests that the error is estimated by

(3.2)#\[ \text{Error in }\Delta_M f(x) \approx \left | \Delta_\text{M} f(x) - \Delta'_\text{M} f(x) \right | \]

Note that the exact value of \(f'(x)\) is not needed to find the error. The right handside can be interpreted as the improvement of the nemerical derivative due to a smaller step size. If there is not much improvement, there is no point to use a further smaller step size.

Based on this error estimate, we can write a code that automatically finds numerical derivative with a desired accuracy using the error estimate.

3.2.2. Adaptation of step size#

What is the best step size? As we saw in the previous section, the result is not accurate if \(h\) is too large. On the other hand, if \(h\) is too small, the result is again bad. We can find an optimal step size based on the estimated error (3.2).

Alrgorithm 2.2.1 Automatic adaptation of step size

  1. Set a value of \textit{tolerance} (allowed error).

  2. Set a reference step size \(h_1\) and evaluate a reference derivative \(g_1 = \Delta_\text{M} f(x)\) using \(h_1\).

  3. Set a new step size \(h_2=h_1/2\) and evaluate a new derivative \(g_2 = \Delta_\text{M} f(x)\) using \(h_2\). 4.Evaluate error \(\delta = |g_2-g_1|\).

  4. If \(\delta < \text{tolerance}\), \(g_2\) is the desired result. Stop the loop.

  5. If not, let \(g_1=g_2\) and \(h_1=h_2\) (previous \(g_2\) and \(h_2\) are now the reference).

  6. Go back to Step 3.

Here the absolute error is used. We can also use a relative error \(\left |\displaystyle \frac{g_2 - g_1}{g_1} \right |\) instead. Then, the tolerance specifies a desired relative error.

Example 2.2.1

Let us calculat the same problem as Example 2.1.1. We numerically evaluate the derivative of \(f(x)=\displaystyle\frac{1}{3}x^3\) at \(x=1\) again. This time, we do not specify the step length \(h\). Instead we specify a tolerance and the program will automatically find an appropriate \(h\) for the given tolerance. Let the absolute error tolerance be 0.001.

import numpy as np
import matplotlib.pyplot as plt
           
def func(x):  # define a function
    return x**3/3
    
tol=0.001

x=1.0
h=1.0   # initial interval
diff_old=(func(x+h)-func(x-h))/(2.0*h)  #  derivative first try
delta=np.finfo(float).max  # any value bigger than tol is OK.
    
print("{0:^10} {1:^16} {2:^16}".format('h','derivative','error'))
while (delta>tol):
    h=h/2.0
    diff_new=(func(x+h)-func(x-h))/(2.0*h)  # improved derivative
    delta=np.abs(diff_new-diff_old)
    print("{0:10.3e} {1:16.10e} {2:16.10e}".format(h,diff_new,delta))
    diff_old=diff_new
            
print("Tolerance is OK.")

    h         derivative         error      
 5.000e-01 1.0833333333e+00 2.5000000000e-01
 2.500e-01 1.0208333333e+00 6.2500000000e-02
 1.250e-01 1.0052083333e+00 1.5625000000e-02
 6.250e-02 1.0013020833e+00 3.9062500000e-03
 3.125e-02 1.0003255208e+00 9.7656250000e-04
Tolerance is OK.

The calculation found that \(h=3.125 \times 10^{-2}\) produces the derivative with four significant figures.

Exercise Solve the same problem with the relative error with tolerance 0.0001.

Last modified on 2/14/2024 by R. Kawai

Contents