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 illustrates that the operator has one zero eigenvalue, which makes sense because constants are in the nullspace and satisfy the periodic boundary conditions. Also making sense, that zero eigenvalue goes away when kappa is nonzero.

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?