Question

I'm trying to calculate the laplacian of a 2d field with periodic boundary conditions (i.e. the divergence of the gradient) and I'm having trouble getting it to work.

I've tried two different methods. First I used an example from scipy.sparse:

def poisson2d(N,epsilon=1.0):
    """
    Return a sparse CSR matrix for the 2d poisson problem
    with standard 5-point finite difference stencil on a
    square N-by-N grid.
    """
    D = (2 + 2*epsilon)*ones(N*N)
    T =  -epsilon * ones(N*N)
    O =  -ones(N*N)
    T[N-1::N] = 0
    return spdiags([D,O,T,T,O],[0,-N,-1,1,N],N*N,N*N).tocoo().tocsr() #eliminate explicit zeros
dx=0.01
grid=arange(-1,1,dx)
N=len(grid)
D=poisson2d(N)
f=lambda x,y: x^2
phi=matrix([[f(x,y) for x in grid] for y in grid])
dphi2=D*phi.reshape((N**2,1))

And I calculated the laplacian using a filter:

dx=0.01
grid=arange(-1,1,dx)
N=len(grid)
f=lambda x,y: x**2
phi=array([[f(x,y) for x in grid] for y in grid])
stencil=[[0,1,0],[1,-4,1],[0,1,0]]
dphi2=filters.convolve(phi,stencil,mode='wrap')

In both cases I get the right answer in the middle of the field which for the function x**2 is just 2, but at the boundaries I get something totally wrong. Why isn't the 'wrap' mode working in the filter case?

Answer 1

i like this answer (click again to cancel)

0

i dont like this answer (click again to cancel)

The structure of the Poisson matrix is incorrect. Below is some code to set it up. Please also let me know where you found the scipy.sparse example so that I can fix it.

import numpy as np
import scipy.sparse as sp

def lapmat(N):
    """Model matrix for solving Laplace's equation,
    discretised on an NxN grid.

    Parameters
    ----------
    N : int
        Grid dimension.

    Returns
    -------
    A : sparse matrix
        Model matrix.

    """
    E = np.ones(N)
    T = sp.dia_matrix(((-E, E*2, -E), (-1, 0, 1)), shape=(N, N))
    I = sp.eye(N, N)

    return sp.kron(T, I) + sp.kron(I, T)

If you use matplotlib.pyplot.spy to visualise the two matrices, you'll notice gaps in the lapmat version that are not present in the example above.

permanent link | flag offensive

updated2010-03-05 07:06:06

asked 2010-03-05 07:04:57
	Stefan 31

add comment

Related questions