maroba / findiff

Python package for numerical derivatives and partial differential equations in any number of dimensions.
MIT License
413 stars 59 forks source link

Accuracy not affecting stencil generation #31

Closed orebas closed 3 years ago

orebas commented 3 years ago

It seems like when I generate a Stencil, the accuracy is being ignored. Repro: x = np.linspace(0, 10, 11) y = np.linspace(0, 10, 11) X, Y = np.meshgrid(x, y, indexing='ij') u = X+Y d1x = FinDiff(0, x[1] - x[0],1,acc=10) stencil1 = d1x.stencil(u.shape) print(stencil1)

Output (this is the same no matter when I change accuracy too): ('L', 'L'): {(0, 0): -1.5, (1, 0): 2.0, (2, 0): -0.5} ('L', 'C'): {(0, 0): -1.5, (1, 0): 2.0, (2, 0): -0.5} ('L', 'H'): {(0, 0): -1.5, (1, 0): 2.0, (2, 0): -0.5} ('C', 'L'): {(-1, 0): -0.5, (0, 0): 0.0, (1, 0): 0.5} ('C', 'C'): {(-1, 0): -0.5, (0, 0): 0.0, (1, 0): 0.5} ('C', 'H'): {(-1, 0): -0.5, (0, 0): 0.0, (1, 0): 0.5} ('H', 'L'): {(-2, 0): 0.5, (-1, 0): -2.0, (0, 0): 1.5} ('H', 'C'): {(-2, 0): 0.5, (-1, 0): -2.0, (0, 0): 1.5} ('H', 'H'): {(-2, 0): 0.5, (-1, 0): -2.0, (0, 0): 1.5}

orebas commented 3 years ago

As a comment, I wanted to let you know that I am finding this package otherwise quite useful, indeed I have not found anything as user-friendly. Thank you.

maroba commented 3 years ago

Thanks for compliment!

Regarding your issue, the acc property doesn't seem to be passed down correctly to the stencil generator. However, until I fix this, you can alternatively attach the acc argument to the stencil call (instead of the FinDiff constructor) and write:

d1x = FinDiff(0, x[1] - x[0], 1)
stencil1 = d1x.stencil(u.shape, acc=4)
print(stencil1)

Output:

('L', 'L'): {(0, 0): -2.083333333333331, (1, 0): 3.9999999999999916, (2, 0): -2.999999999999989, (3, 0): 1.3333333333333268, (4, 0): -0.24999999999999858}
('L', 'C'): {(0, 0): -2.083333333333331, (1, 0): 3.9999999999999916, (2, 0): -2.999999999999989, (3, 0): 1.3333333333333268, (4, 0): -0.24999999999999858}
('L', 'H'): {(0, 0): -2.083333333333331, (1, 0): 3.9999999999999916, (2, 0): -2.999999999999989, (3, 0): 1.3333333333333268, (4, 0): -0.24999999999999858}
('C', 'L'): {(-2, 0): 0.08333333333333333, (-1, 0): -0.6666666666666666, (0, 0): 0.0, (1, 0): 0.6666666666666666, (2, 0): -0.08333333333333333}
('C', 'C'): {(-2, 0): 0.08333333333333333, (-1, 0): -0.6666666666666666, (0, 0): 0.0, (1, 0): 0.6666666666666666, (2, 0): -0.08333333333333333}
('C', 'H'): {(-2, 0): 0.08333333333333333, (-1, 0): -0.6666666666666666, (0, 0): 0.0, (1, 0): 0.6666666666666666, (2, 0): -0.08333333333333333}
('H', 'L'): {(-4, 0): 0.24999999999999958, (-3, 0): -1.3333333333333313, (-2, 0): 2.9999999999999956, (-1, 0): -3.999999999999996, (0, 0): 2.0833333333333317}
('H', 'C'): {(-4, 0): 0.24999999999999958, (-3, 0): -1.3333333333333313, (-2, 0): 2.9999999999999956, (-1, 0): -3.999999999999996, (0, 0): 2.0833333333333317}
('H', 'H'): {(-4, 0): 0.24999999999999958, (-3, 0): -1.3333333333333313, (-2, 0): 2.9999999999999956, (-1, 0): -3.999999999999996, (0, 0): 2.0833333333333317}
orebas commented 3 years ago

Cool! Thank you so much for taking a look.

Just to let you know, I was actually trying to calculate mixed partials. (Specifically d^2/dxdy on a 2d grid, but to higher order.) I ended up applying one stencil after another (manually) in my code.

Findiff doesn't seem to like it. I am not surprised because this is harder to implement. I'm just letting you know - as above my "acute" need is solved. If it would be helpful I can upload this to the github page, or maybe even try to work on it. I think you need to "tensor" the stencils somehow. CODE:

x = np.linspace(0, 10, 11) y = np.linspace(0, 10, 11) X, Y = np.meshgrid(x, y, indexing='ij') u = X+Y dx = x[1]-x[0] dy = y[1]-y[0] d1x = FinDiff((0, dx), (1, dy)) stencil1 = d1x.stencil(u.shape) mat1=d1x.matrix(u.shape) print(stencil1) print(mat1)

Error: File "C:/Users/oreba/PycharmProjects/untitled/HelloWorld.py", line 13, in main stencil1 = d1x.stencil(u.shape) File "C:\ProgramData\Anaconda3\envs\MainEnv\lib\site-packages\findiff\operators.py", line 99, in stencil return self.pds.stencil(shape, h, acc, old_stl) AttributeError: 'Mul' object has no attribute 'stencil' On Mon, Dec 21, 2020 at 2:53 AM Matthias Baer notifications@github.com wrote:

Thanks for compliment!

Regarding your issue, the acc property doesn't seem to be passed down correctly to the stencil generator. However, until I fix this, you can alternatively attach the acc argument to the stencil call (instead of the FinDiff constructor) and write:

d1x = FinDiff(0, x[1] - x[0], 1)stencil1 = d1x.stencil(u.shape, acc=4)print(stencil1)

Output:

('L', 'L'): {(0, 0): -2.083333333333331, (1, 0): 3.9999999999999916, (2, 0): -2.999999999999989, (3, 0): 1.3333333333333268, (4, 0): -0.24999999999999858} ('L', 'C'): {(0, 0): -2.083333333333331, (1, 0): 3.9999999999999916, (2, 0): -2.999999999999989, (3, 0): 1.3333333333333268, (4, 0): -0.24999999999999858} ('L', 'H'): {(0, 0): -2.083333333333331, (1, 0): 3.9999999999999916, (2, 0): -2.999999999999989, (3, 0): 1.3333333333333268, (4, 0): -0.24999999999999858} ('C', 'L'): {(-2, 0): 0.08333333333333333, (-1, 0): -0.6666666666666666, (0, 0): 0.0, (1, 0): 0.6666666666666666, (2, 0): -0.08333333333333333} ('C', 'C'): {(-2, 0): 0.08333333333333333, (-1, 0): -0.6666666666666666, (0, 0): 0.0, (1, 0): 0.6666666666666666, (2, 0): -0.08333333333333333} ('C', 'H'): {(-2, 0): 0.08333333333333333, (-1, 0): -0.6666666666666666, (0, 0): 0.0, (1, 0): 0.6666666666666666, (2, 0): -0.08333333333333333} ('H', 'L'): {(-4, 0): 0.24999999999999958, (-3, 0): -1.3333333333333313, (-2, 0): 2.9999999999999956, (-1, 0): -3.999999999999996, (0, 0): 2.0833333333333317} ('H', 'C'): {(-4, 0): 0.24999999999999958, (-3, 0): -1.3333333333333313, (-2, 0): 2.9999999999999956, (-1, 0): -3.999999999999996, (0, 0): 2.0833333333333317} ('H', 'H'): {(-4, 0): 0.24999999999999958, (-3, 0): -1.3333333333333313, (-2, 0): 2.9999999999999956, (-1, 0): -3.999999999999996, (0, 0): 2.0833333333333317}

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/maroba/findiff/issues/31#issuecomment-748825467, or unsubscribe https://github.com/notifications/unsubscribe-auth/AK37HISLTHMNWV7SWC7XTUTSV35GVANCNFSM4VDTHZTQ .

maroba commented 3 years ago

Yes, you are right. Stencils for mixed partial derivatives were not yet supported.

In the meantime I have fixed the issue that you reported and also implemented the stencils for mixed partial derivatives. Both is uploaded on github and pypi as version v0.8.6.

By the way, do you explicitly need the stencils or do you just want to apply the derivatives?

orebas commented 3 years ago

I am using the stencils - coding some finite difference stuff in C++. Believe it or not I cannot find a similar library in C++. (I don't want the whole sparse matrix, my data is too large for that and I am already memory constrained.)

Thanks for the incredibly prompt response, I am very grateful. --Oren Bassik

On Tue, Dec 22, 2020 at 9:55 AM Matthias Baer notifications@github.com wrote:

Yes, you are right. Stencils for mixed partial derivatives were not yet supported.

In the meantime I have fixed the issue that you reported and also implemented the stencils for mixed partial derivatives. Both is uploaded on github and pypi as version v0.8.6.

By the way, do you explicitly need the stencils or do you just want to apply the derivatives?

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/maroba/findiff/issues/31#issuecomment-749580103, or unsubscribe https://github.com/notifications/unsubscribe-auth/AK37HIVTKTGNUFQYXUUQCMLSWCXNTANCNFSM4VDTHZTQ .