Problems with the stencil matrix

[sotakao]

Hi, I've been trying to produce the stencil matrix for the Helmholtz operator

Helmholtz = Identity() - FinDiff(0, dx, 2)

I've noticed that general operators such as Minus does not have the stencil attribute to be produced. Would it be possible to include this in future version?

I've also found an error in the following code snippet of README:

import numpy as np
from findiff import FinDiff

x, y = np.linspace(0, 1, 100)  #Supposed to be x, y = np.linspace(0, 1, 100), np.linspace(0, 1, 100)?
X, Y = np.meshgrid(x, y, indexing='ij')
u = X**3 + Y**3
laplace_2d = FinDiff(0, x[1]-x[0], 2) + FinDiff(1, y[1]-y[0], 2)

stencil = laplace_2d.stencil(u.shape)

print(stencil)  #Supposed to be print(stencil.data)?

does not yield the desired result. Instead, I get:

{('L', 'L'): {(0, 0): 39203.99999999999, (0, 1): -49004.99999999999, (0, 2): 39203.99999999999, (0, 3): -9800.999999999998, (1, 0): -49004.99999999999, (2, 0): 39203.99999999999, (3, 0): -9800.999999999998}, ('L', 'C'): {(0, -1): 9800.999999999998, (0, 1): 9800.999999999998, (1, 0): -49004.99999999999, (2, 0): 39203.99999999999, (3, 0): -9800.999999999998}, ('L', 'H'): {(0, -3): -9800.999999999996, (0, -2): 39203.999999999985, (0, -1): -49004.99999999998, (0, 0): 39203.999999999985, (1, 0): -49004.99999999999, (2, 0): 39203.99999999999, (3, 0): -9800.999999999998}, ('C', 'L'): {(-1, 0): 9800.999999999998, (0, 1): -49004.99999999999, (0, 2): 39203.99999999999, (0, 3): -9800.999999999998, (1, 0): 9800.999999999998}, ('C', 'C'): {(-1, 0): 9800.999999999998, (0, -1): 9800.999999999998, (0, 0): -39203.99999999999, (0, 1): 9800.999999999998, (1, 0): 9800.999999999998}, ('C', 'H'): {(-1, 0): 9800.999999999998, (0, -3): -9800.999999999996, (0, -2): 39203.999999999985, (0, -1): -49004.99999999998, (0, 0): -3.637978807091713e-12, (1, 0): 9800.999999999998}, ('H', 'L'): {(-3, 0): -9800.999999999996, (-2, 0): 39203.999999999985, (-1, 0): -49004.99999999998, (0, 0): 39203.999999999985, (0, 1): -49004.99999999999, (0, 2): 39203.99999999999, (0, 3): -9800.999999999998}, ('H', 'C'): {(-3, 0): -9800.999999999996, (-2, 0): 39203.999999999985, (-1, 0): -49004.99999999998, (0, -1): 9800.999999999998, (0, 0): -3.637978807091713e-12, (0, 1): 9800.999999999998}, ('H', 'H'): {(-3, 0): -9800.999999999996, (-2, 0): 39203.999999999985, (-1, 0): -49004.99999999998, (0, -3): -9800.999999999996, (0, -2): 39203.999999999985, (0, -1): -49004.99999999998, (0, 0): 39203.999999999985}}

I believe x, y = np.linspace(0, 1, 100), np.linspace(0, 1, 100) should be x, y = np.linspace(0, 99, 100), np.linspace(0, 99, 100) to get the desired result.

Thank you for your time and for developing this nice package.

Update: I've edited this and now it is not so much a 'problem' with the stencil matrix but rather a feature request.

[maroba]

Hi,

I'm not sure if you mean stencil or matrix representation. Anyways, let's look at both.

First of all, yes it is

x, y = np.linspace(0, 1, 100), np.linspace(0, 1, 100)

or simply

x = y = np.linspace(0, 1, 100)

For the stencil, consider the simplified version with only a 5x5 grid of grid spacing 1:

laplace_2d = FinDiff(0, 1, 2) + FinDiff(1, 1, 2)
stencil = laplace_2d.stencil((5, 5))
print(stencil)  # as of version v0.9.1 (just released) it is just stencil an no longer stencil.data

which yields

{('L', 'L'): {(0, 0): 4.0, (0, 1): -5.0, (0, 2): 4.0, (0, 3): -1.0, (1, 0): -5.0, (2, 0): 4.0, (3, 0): -1.0},
 ('L', 'C'): {(0, -1): 1.0, (0, 1): 1.0, (1, 0): -5.0, (2, 0): 4.0, (3, 0): -1.0},
 ('L', 'H'): {(0, -3): -1.0, (0, -2): 4.0, (0, -1): -5.0, (0, 0): 4.0, (1, 0): -5.0,
              (2, 0): 4.0, (3, 0): -1.0},
 ('C', 'L'): {(-1, 0): 1.0, (0, 1): -5.0, (0, 2): 4.0, (0, 3): -1.0, (1, 0): 1.0},
 ('C', 'C'): {(-1, 0): 1.0, (0, -1): 1.0, (0, 0): -4.0, (0, 1): 1.0, (1, 0): 1.0},
 ('C', 'H'): {(-1, 0): 1.0, (0, -3): -1.0, (0, -2): 4.0, (0, -1): -5.0, (0, 0): 8.881784197001252e-16,
              (1, 0): 1.0},
 ('H', 'L'): {(-3, 0): -1.0, (-2, 0): 4.0, (-1, 0): -5.0, (0, 0): 4.0, (0, 1): -5.0,
              (0, 2): 4.0, (0, 3): -1.0},
 ('H', 'C'): {(-3, 0): -1.0, (-2, 0): 4.0, (-1, 0): -5.0, (0, -1): 1.0, (0, 0): 8.881784197001252e-16,
              (0, 1): 1.0},
 ('H', 'H'): {(-3, 0): -1.0, (-2, 0): 4.0, (-1, 0): -5.0, (0, -3): -1.0, (0, -2): 4.0,
              (0, -1): -5.0, (0, 0): 4.0}}

which looks ok. In the interior of the grid, the (C,C) stencil applies, near the boundary one of the others. L means low index boundary region, H means high index boundary region. The stencil in your case looks strange because it contains 1 / grid_spacing **2 terms, which is just 1 in my example.

If you mean matrix, then this is how you can get it:

laplace_2d = FinDiff(0, dx, 2) + FinDiff(1, dy, 2)
mat = laplace_2d.matrix(u.shape)   #  returns a scipy sparse matrix

# apply matrix representation to u:
mat.dot(u.reshape(-1)).reshape(u.shape)

Regarding the Helmholtz operator in 1D, the following works:

from findiff import *
dx = 1
H = Identity() - FinDiff(0, dx, 2)
print(H.matrix((10,)).toarray())

which gives

[[-1.  5. -4.  1.  0.  0.  0.  0.  0.  0.]
 [-1.  3. -1.  0.  0.  0.  0.  0.  0.  0.]
 [ 0. -1.  3. -1.  0.  0.  0.  0.  0.  0.]
 [ 0.  0. -1.  3. -1.  0.  0.  0.  0.  0.]
 [ 0.  0.  0. -1.  3. -1.  0.  0.  0.  0.]
 [ 0.  0.  0.  0. -1.  3. -1.  0.  0.  0.]
 [ 0.  0.  0.  0.  0. -1.  3. -1.  0.  0.]
 [ 0.  0.  0.  0.  0.  0. -1.  3. -1.  0.]
 [ 0.  0.  0.  0.  0.  0.  0. -1.  3. -1.]
 [ 0.  0.  0.  0.  0.  0.  1. -4.  5. -1.]]

which looks correct. For actual computation you need a reasonable grid spacing, of course.

The stencil method for Helmholtz example really gives an error. I will fix it in the next days as version 0.9.2.

Edit: What part is a feature request?

[maroba]

Ok, I just released v0.9.2 with the necessary fix.

[sotakao]

Thank you very much! I just checked that I can create the stencil for the Helmholtz operator with v0.9.2. Have a nice day.

[maroba]

Thanks for checking. Have a nice day, too!