Complex-step derivative approximation (implemented in Python, NumPy, matplotlib)

complex_step_derivative.py

# coding: utf-8

'''The following code reproduces an example of the paper:
'The Complex-Step Derivative Approximation'
by Joaquim R. R. A. MARTINS, Peter STURDZA and Juan J. Alonso published in 2003.

License: MIT
Author: FJ Navarro Brull
'''



import numpy as np
import matplotlib.pyplot as plt

# Define function to be tested
def F(x):
    return np.exp(x)/((np.cos(x))**3 + np.sin(x)**3)

# Number of evaluation points
N = 1000

x = np.linspace(-np.pi/4,np.pi/2, N, endpoint=True) 
y =  F(x)

# Where the derivative will be calculated
point_x = np.pi/4
point_y = F(point_x)

# Take a look at the function being analysed
plt.figure()
plt.plot(x,y)
plt.plot(point_x, point_y, 'o')

plt.xlim((-np.pi/4, np.pi/2))
plt.ylim((0, 6))

plt.ylabel('F(x)')
plt.xlabel('x')
plt.annotate(xy=(-0.5,4),s=r'$F(x)=\frac{e^x}{\cos(x)^3+\sin(x)^3}$',fontsize=20)

plt.plot()
plt.show()


# Define derivatives

def dF_complex(x,h):
    '''Complex step aproximation'''
    #x =  np.array(x, dtype='cfloat')
    return np.imag(F(x + 1j*h))/h

def dF_ff(x,h):
    '''Finite forward-difference approximation'''
    return (F(x+h) - F(x))/h

def dF_cf(x,h):
    '''Finite central-difference approximation'''
    return (F(x+h) - F(x-h))/(2*h)

def dF_analytic(x):
    '''Analytic derivative (to obtain the exact value)'''
    return ((np.exp(x)*(np.cos(3*x) + np.sin(3*x)/2 + (3*np.sin(x))/2)) / 
            (np.cos(x)**3 + np.sin(x)**3)**2)


# Test derivatives as a function of derivative step size

point_x = np.pi/4
h_values = np.logspace(-15,-1,15, endpoint=True)

# Preallocate values in memory
error_complex = np.zeros(h_values.shape)
error_cf = np.zeros(h_values.shape)
error_ff = np.zeros(h_values.shape)

for i, h in enumerate(h_values):
    error_complex[i] = np.abs(dF_complex(point_x, h) - dF_analytic(point_x))/np.abs(dF_analytic(point_x))
    error_ff[i] = np.abs(dF_ff(point_x, h) - dF_analytic(point_x))/np.abs(dF_analytic(point_x))
    error_cf[i] = np.abs(dF_cf(point_x, h) - dF_analytic(point_x))/np.abs(dF_analytic(point_x))

# Plot the results

plt.figure()
plt.loglog(h_values, error_complex, label='Complex step')
plt.loglog(h_values, error_ff, label='Forward difference')
plt.loglog(h_values, error_cf, label='Central-difference')

plt.gca().invert_xaxis()
plt.legend()
plt.ylabel('Normalized error')
plt.xlabel('Step size (h)')

plt.show()

I am currently trying to apply the complex step method to the heat equation. Wish to have some ideas on how to go about it, Thanks