pdefields.py

from scipy.sparse import linalg
import numpy as np
import map_utils

def dim2w(dim):
    return np.arange(-np.ceil((dim-1)/2.), np.floor((dim-1)/2.)+1)

class fractional_modified_laplacian(linalg.LinearOperator):
    
    def __init__(self, dims, kappa, alpha):
        # Remember the size of the grid, kappa, and alpha.
        self.dims = dims
        self.kappa = kappa
        self.alpha = alpha
        self.xsize = np.prod(dims)
        
        self.wsq = dim2w(dims[0])**2
        for d in self.dims[1:]:
            self.wsq = np.add.outer(self.wsq, dim2w(d)**2)
        self.wsq = np.fft.ifftshift(self.wsq)
        
        # let the LinearOperator class do whatever it needs to do to initialize.
        linalg.LinearOperator.__init__(self, (self.xsize,self.xsize), None, dtype='float')
    
    def __call__(self, x):
        return self.matvec(x.ravel()).reshape(self.dims)
    
    def matvec(self, x):
        "Matvec computes Lx for given input x."
        # x is coming in as a vector, reshape it to an n-dimensional array.
        x_ = x.reshape(self.dims)
        
        # Apply an n-dimensional fft.
        k = np.fft.fftn(x_)
        
        # Apply the fractional modified Laplacian in Fourier space.
        # The fractional modified Laplacian is (\kappa^2 - \Delta)^{\alpha/2}
        k_trans = (self.kappa**2+self.wsq)**(self.alpha/2.)*k

        # Apply an n-dimensional inverse fft, reshape to a vector and return.
        y_ = np.fft.ifftn(k_trans).real
        return y_.ravel()
        
def expand_matrix(op):
    out = np.empty((op.xsize, op.xsize))
    x = np.zeros(op.xsize)
    for i in xrange(op.xsize):
        x[i] = 1
        out[i] = op.matvec(x)
        x[i] = 0
    out[np.where(np.abs(out)<1.e-10)]=0
    return out

def genfield(L):
    """
    Solves the system Lx = W using the gmres method, see 
    http://docs.scipy.org/doc/scipy/reference/sparse.linalg.html#scipy.sparse.linalg.gmres
    """
    w = np.random.normal(size=L.dims)
    return linalg.lgmres(L, w.ravel())[0].reshape(L.dims)

if __name__ == '__main__':
    import euclidean
    n = 101
    L = euclidean.fractional_modified_laplacian((n,n), 10, 2)
    import pylab as pl
    # x = np.cos(np.linspace(-1,1,n)*np.pi)
    # x = np.add.outer(x,x)
    pl.clf()
    # pl.subplot(1,2,1)
    # pl.imshow(x,interpolation='nearest')
    # pl.colorbar()
    # pl.subplot(1,2,2)
    # Lx = L(x)
    # pl.imshow(Lx,interpolation='nearest')
    # pl.colorbar()
    
    # Get L as a full matrix
    # A = expand_matrix(L)
    # e,v = np.linalg.eigh(A)

    # Number of zero eigenvalues.
    # print (np.abs(e)<1.e-10).sum()
    
    f = genfield(L)
    pl.imshow(f,interpolation='nearest')

The new version seems to work, broadly speaking, in spite of the boundary effects. However, it's really slow for smooth fields: small kappa, large alpha, high resolution. Can the method be improved?