• An Effective Algorithm for Computing the Numerical Range by Carl C. Cowen and Elad Harel
• Field of values and numerical abscissa by Nick Trefethen, November

import numpy as np
import matplotlib.mlab as mlab
 
 
def numerical_range(A,resolution=0.01):
    """
    Function implements algorithm for calculation of numerical range
    http://www.math.iupui.edu/~ccowen/Downloads/33NumRange.html
    """
    A=np.asmatrix(A)
    th=np.arange(0,2*np.pi+resolution,resolution)
    k=0
    w=[]
    for j in th:
        Ath=np.exp(1j*-j)*A
        Hth=(Ath+Ath.H)/2
        e,r=np.linalg.eigh(Hth)
        r=np.matrix(r)
        e=np.real(e)
        m=e.max()
        s=mlab.find(e==m)
        if np.size(s)==1:
            w.append(np.matrix.item(r[:,s].H*A*r[:,s]))
        else:
            Kth=1j*(Hth-Ath)
            pKp=r[:,s].H*Kth*r[:,s]
            ee,rr=np.linalg.eigh(pKp)
            rr=np.matrix(rr)
            ee=np.real(ee)
            mm=ee.min()
            sm=mlab.find(ee==mm)
            temp=rr[:,sm[0]].H*r[:,s].H*A*r[:,s]*rr[:,sm[0]]
            w.append(temp[0,0])
            k+=1
            mM=ee.max()
            sM=mlab.find(ee==mM)
            temp=rr[:,sM[0]].H*r[:,s].H*A*r[:,s]*rr[:,sM[0]]
            w.append(temp[0,0])
        k+=1
    return w