Here's a simple convergence study using

`cos()`

, whose second derivative is `-cos()`

. We'll use python to loop over an array with the number of points per wavelength we want to use in our approximation and then calculate a residual, or error, between the analytical solution and our approximation.for ni in n:

x = 2*pi*arange(0,ni)/(ni-1)

h = x[1]-x[0]

f = cos(x)

analytical_d2f = -cos(x)

# take a positive imaginary step

x_step = x + complex(0,h)

f_pc = cos(x_step)

# take a negative imaginary step

x_step = x - complex(0,h)

f_mc = cos(x_step)

# take a positive real step

x_step = x + h

f_ph = cos(x_step)

# take a negative real step

x_step = x - h

f_mh = cos(x_step)

approx_d2f_c = ( f_ph + f_mh - f_pc.real - f_mc.real)/(2*h*h)

approx_d2f_h = ( f_ph + f_mh - 2*f )/(h*h)

resid_c[i] = sqrt( sum( (analytical_d2f - approx_d2f_c)**2 ) )

resid_h[i] = sqrt( sum( (analytical_d2f - approx_d2f_h)**2 ) )

i += 1

The approximated and analytical second derivatives are shown below for 10 points per wavelength using both a standard central difference with only a real step as well as the complex step formula.

The convergence behaviour for the two approximations is shown below.

Obviously we're getting much better convergence with our complex step formula.