of extreme eigenvalues and/or eigenvectors can be computed.

>>> lap = LaplacianNd(grid_shape, boundary_conditions='periodic')
>>> lap.eigenvalues()
array([-4., -3., -3., -1., -1.,  0.])
>>> lap.eigenvalues()[-2:]
array([-1.,  0.])
>>> lap.eigenvalues(2)
array([-1.,  0.])
>>> lap.eigenvectors(1)
array([[0.40824829],
       [0.40824829],
       [0.40824829],
       [0.40824829],
       [0.40824829],
       [0.40824829]])
>>> lap.eigenvectors(2)
array([[ 0.5       ,  0.40824829],
       [ 0.        ,  0.40824829],
       [-0.5       ,  0.40824829],
       [-0.5       ,  0.40824829],
       [ 0.        ,  0.40824829],
       [ 0.5       ,  0.40824829]])
>>> lap.eigenvectors()
array([[ 0.40824829,  0.28867513,  0.28867513,  0.5       ,  0.5       ,
         0.40824829],
       [-0.40824829, -0.57735027, -0.57735027,  0.        ,  0.        ,
         0.40824829],
       [ 0.40824829,  0.28867513,  0.28867513, -0.5       , -0.5       ,
         0.40824829],
       [-0.40824829,  0.28867513,  0.28867513, -0.5       , -0.5       ,
         0.40824829],
       [ 0.40824829, -0.57735027, -0.57735027,  0.        ,  0.        ,
         0.40824829],
       [-0.40824829,  0.28867513,  0.28867513,  0.5       ,  0.5       ,
         0.40824829]])

The two-dimensional Laplacian is illustrated on a regular grid with
``grid_shape = (2, 3)`` points in each dimension.

>>> grid_shape = (2, 3)
>>> n = np.prod(grid_shape)

Numeration of grid points is as follows:

>>> np.arange(n).reshape(grid_shape + (-1,))
array([[[0],
        [1],
        [2]],
<BLANKLINE>
       [[3],
        [4],
        [5]]])

Each of the boundary conditions ``'dirichlet'``, ``'periodic'``, and
``'neumann'`` is illustrated separately; with ``'dirichlet'``

>>> lap = LaplacianNd(grid_shape, boundary_conditions='dirichlet')
>>> lap.tosparse()
<Compressed Sparse Row sparse array of dtype 'int8'
    with 20 stored elements and shape (6, 6)>
>>> lap.toarray()
array([[-4,  1,  0,  1,  0,  0],
       [ 1, -4,  1,  0,  1,  0],
       [ 0,  1, -4,  0,  0,  1],
       [ 1,  0,  0, -4,  1,  0],
       [ 0,  1,  0,  1, -4,  1],
       [ 0,  0,  1,  0,  1, -4]], dtype=int8)
>>> np.array_equal(lap.matmat(np.eye(n)), lap.toarray())
True
>>> np.array_equal(lap.tosparse().toarray(), lap.toarray())
True
>>> lap.eigenvalues()
array([-6.41421356, -5.        , -4.41421356, -3.58578644, -3.        ,
       -1.58578644])
>>> eigvals = eigvalsh(lap.toarray().astype(np.float64))
>>> np.allclose(lap.eigenvalues(), eigvals)
True
>>> np.allclose(lap.toarray() @ lap.eigenvectors(),
...             lap.eigenvectors() @ np.diag(lap.eigenvalues()))
True

with ``'periodic'``

>>> lap = LaplacianNd(grid_shape, boundary_conditions='periodic')
>>> lap.tosparse()
<Compressed Sparse Row sparse array of dtype 'int8'
    with 24 stored elements and shape (6, 6)>
>>> lap.toarray()
    array([[-4,  1,  1,  2,  0,  0],
           [ 1, -4,  1,  0,  2,  0],
           [ 1,  1, -4,  0,  0,  2],
           [ 2,  0,  0, -4,  1,  1],
           [ 0,  2,  0,  1, -4,  1],
           [ 0,  0,  2,  1,  1, -4]], dtype=int8)
>>> np.array_equal(lap.matmat(np.eye(n)), lap.toarray())
True
>>> np.array_equal(lap.tosparse().toarray(), lap.toarray())
True
>>> lap.eigenvalues()
array([-7., -7., -4., -3., -3.,  0.])
>>> eigvals = eigvalsh(lap.toarray().astype(np.float64))
>>> np.allclose(lap.eigenvalues(), eigvals)
True
>>> np.allclose(lap.toarray() @ lap.eigenvectors(),
...             lap.eigenvectors() @ np.diag(lap.eigenvalues()))
True

and with ``'neumann'``

>>> lap = LaplacianNd(grid_shape, boundary_conditions='neumann')
>>> lap.tosparse()
<Compressed Sparse Row sparse array of dtype 'int8'
    with 20 stored elements and shape (6, 6)>
>>> lap.toarray()
array([[-2,  1,  0,  1,  0,  0],
       [ 1, -3,  1,  0,  1,  0],
       [ 0,  1, -2,  0,  0,  1],
       [ 1,  0,  0, -2,  1,  0],
       [ 0,  1,  0,  1, -3,  1],
       [ 0,  0,  1,  0,  1, -2]], dtype=int8)
>>> np.array_equal(lap.matmat(np.eye(n)), lap.toarray())
True
>>> np.array_equal(lap.tosparse().toarray(), lap.toarray())
True
>>> lap.eigenvalues()
array([-5., -3., -3., -2., -1.,  0.])
>>> eigvals = eigvalsh(lap.toarray().astype(np.float64))
>>> np.allclose(lap.eigenvalues(), eigvals)
True
>>> np.allclose(lap.toarray() @ lap.eigenvectors(),
...             lap.eigenvectors() @ np.diag(lap.eigenvalues()))
True

Ú