1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
|
import numpy as np
class FiniteDifferenceUniform(object):
"""
Finite difference operator on a uniform grid.
:param int diff_order: Order of the derivative to take. Currently only 1 and
2 are supported.
:param int acc_order: Accuracy order of the operator. Derivatives of
polynomials up to this order will be computed
exactly.
"""
# public
diff_order = None
acc_order = None
"""
For each source point, this number of points in each direction is needed
to compute the derivative to desired order. Points at the boundary are
computed with lower order operators.
"""
stencil = None
# private
_coeffs = None
_div = None
_diff_table = [
{
2 : {
'stencil' : 1,
'coeffs' : (-1., 0., 1.),
'div' : 2.,
},
4 : {
'stencil' : 2,
'coeffs' : (1., -8., 0., 8., -1.),
'div' : 12.,
},
6 : {
'stencil' : 3,
'coeffs' : (-1., 9., -45., 0., 45., -9., 1.),
'div' : 60.,
},
8 : {
'stencil' : 4,
'coeffs' : (3., -32., 168., -672., 0., 672., -168., 32., -3.),
'div' : 840.,
},
},
{
2 : {
'stencil' : 1,
'coeffs' : (1., -2., 1.),
'div' : 1.,
},
4 : {
'stencil' : 2,
'coeffs' : (-1., 16., -30., 16., -1.),
'div' : 12.,
},
6 : {
'stencil' : 3,
'coeffs' : (2., -27., 270., -490., 270., -27., 2.),
'div' : 180.,
},
8 : {
'stencil' : 4,
'coeffs' : (-9., 128., -1008., 8064., -14350., 8064., -1008., 128., -9.),
'div' : 5040.,
},
},
]
def __init__(self, diff_order, acc_order):
self.diff_order = diff_order
self.acc_order = acc_order
self._coeffs = self._diff_table[diff_order - 1][acc_order]['coeffs']
self._div = self._diff_table[diff_order - 1][acc_order]['div']
self.stencil = self._diff_table[diff_order - 1][acc_order]['stencil']
def __call__(self, arr, axis, dx, bound = 'both'):
"""
Evaluate the derivative for an array.
:param array_like arr: Array to differentiate.
:param int axis: Axis along which to differentiate.
:param float dx: Spacing between array elements.
:return: derivative of arr
:rtype: Numpy array, same shape and type as arr.
"""
coeff = self._coeffs
div = self._div
stencil = self.stencil
bound_lower = bound == 'lower' or bound == 'both'
bound_upper = bound == 'upper' or bound == 'both'
slices_src = []
for i in range(stencil * 2 + 1):
sl = [slice(None) for _ in arr.shape]
start = i
end = -stencil * 2 + i
if end == 0:
end = None
sl[axis] = slice(start, end)
slices_src.append(tuple(sl))
slicelist_ret = [slice(None)] * len(arr.shape)
slicelist_ret[axis] = slice(stencil, -stencil)
ret = np.zeros_like(arr)
for s, c in zip(slices_src, coeff):
ret[tuple(slicelist_ret)] += c * arr[tuple(s)]
ret[tuple(slicelist_ret)] /= (div * (dx ** self.diff_order))
if bound_lower:
if self.acc_order > 2:
fd_lower = FiniteDifferenceUniform(self.diff_order, self.acc_order - 2)
slicelist_edge_src = [slice(None) for _ in arr.shape]
slicelist_edge_src[axis] = slice(0, 2 * stencil)
diff_edge = fd_lower(arr[tuple(slicelist_edge_src)], axis, dx, bound = 'lower')
slicelist_edge_dst = [slice(None) for _ in arr.shape]
slicelist_edge_dst[axis] = slice(0, stencil)
ret[tuple(slicelist_edge_dst)] = diff_edge[tuple(slicelist_edge_dst)]
else:
slicelist_edge_dst = [slice(None) for _ in arr.shape]
slicelist_edge_dst[axis] = slice(0, 1)
slicelist_edge_p1 = [slice(None) for _ in arr.shape]
slicelist_edge_p1[axis] = slice(1, 2)
slicelist_edge_m1 = [slice(None) for _ in arr.shape]
slicelist_edge_m1[axis] = slice(0, 1)
ret[tuple(slicelist_edge_dst)] = (arr[tuple(slicelist_edge_p1)] - arr[tuple(slicelist_edge_m1)]) / dx
if bound_upper:
if self.acc_order > 2:
fd_lower = FiniteDifferenceUniform(self.diff_order, self.acc_order - 2)
slicelist_edge_src = [slice(None) for _ in arr.shape]
slicelist_edge_src[axis] = slice(-2 * stencil, None)
diff_edge = fd_lower(arr[tuple(slicelist_edge_src)], axis, dx, bound = 'upper')
slicelist_edge_dst = [slice(None) for _ in arr.shape]
slicelist_edge_dst[axis] = slice(-stencil, None)
ret[tuple(slicelist_edge_dst)] = diff_edge[tuple(slicelist_edge_dst)]
else:
slicelist_edge_dst = [slice(None) for _ in arr.shape]
slicelist_edge_dst[axis] = slice(-1, None)
slicelist_edge_p1 = [slice(None) for _ in arr.shape]
slicelist_edge_p1[axis] = slice(-1, None)
slicelist_edge_m1 = [slice(None) for _ in arr.shape]
slicelist_edge_m1[axis] = slice(-2, -1)
ret[tuple(slicelist_edge_dst)] = (arr[tuple(slicelist_edge_p1)] - arr[tuple(slicelist_edge_m1)]) / dx
return ret
fd2 = FiniteDifferenceUniform(1, 2)
fd4 = FiniteDifferenceUniform(1, 4)
fd8 = FiniteDifferenceUniform(1, 8)
fd22 = FiniteDifferenceUniform(2, 2)
fd24 = FiniteDifferenceUniform(2, 4)
fd28 = FiniteDifferenceUniform(2, 8)
|
