Skip to content

Instantly share code, notes, and snippets.

@franktoffel
Last active July 18, 2022 02:33
Show Gist options
  • Save franktoffel/e65320e553a9c02870931bb288684ae8 to your computer and use it in GitHub Desktop.
Save franktoffel/e65320e553a9c02870931bb288684ae8 to your computer and use it in GitHub Desktop.
Complex-step derivative approximation (implemented in Python, NumPy, matplotlib)
# 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()
@franktoffel
Copy link
Author

function
error

@tstevye
Copy link

tstevye commented Mar 13, 2018

very interesting post here. Thanks for this.

@tstevye
Copy link

tstevye commented Mar 13, 2018

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment