Skip to content

Commit

Permalink
refactor(diff): rename indices to offsets
Browse files Browse the repository at this point in the history
  • Loading branch information
alexfikl committed Sep 19, 2023
1 parent 851cbac commit 1fbe547
Show file tree
Hide file tree
Showing 2 changed files with 20 additions and 20 deletions.
34 changes: 17 additions & 17 deletions pycaputo/differentiation/finite_difference.py
Original file line number Diff line number Diff line change
Expand Up @@ -28,10 +28,10 @@ class Stencil:
.. math::
\frac{\mathrm{d}^n\, f}{\mathrm{d}\, x^n}
\approx \frac{1}{h^n} \sum_{i \in \text{indices}} a_i f_i
\frac{\mathrm{d}^n\, f}{\mathrm{d}\, x^n}(x_m)
\approx \frac{1}{h^n} \sum_{k \in \text{offsets}} a_k f_{m + k}
where :math:`a_i` are the given coefficients :attr:`coeffs` and :math:`f_i`
where :math:`a_k` are the given coefficients :attr:`coeffs` and :math:`f_m`
are the point function evaluations. The approximation is to an order and
truncation error can be retrieved from :attr:`trunc`.
"""
Expand All @@ -40,17 +40,17 @@ class Stencil:
derivative: int
#: Coefficients used in the stencil.
coeffs: Array
#: Indices around the centered :math:`0` used in the stencil.
indices: Array
#: Offsets around the centered :math:`0` used in the stencil.
offsets: Array

@cached_property
def padded_coeffs(self) -> Array:
"""Padded coefficients that are symmetric around the :math:`0`
index and can be easily applied as a convolution.
"""
n = np.max(np.abs(self.indices))
n = np.max(np.abs(self.offsets))
coeffs = np.zeros(2 * n + 1, dtype=self.coeffs.dtype)
coeffs[n + self.indices] = self.coeffs
coeffs[n + self.offsets] = self.coeffs

return coeffs

Expand Down Expand Up @@ -81,9 +81,9 @@ def determine_stencil_truncation_error(
.. math::
\frac{\mathrm{d}^n\, f}{\mathrm{d}\, x^n}
- \sum_{i \in \text{indices}} \frac{a_i}{h^n} f_i
= c \frac{\mathrm{d}^p\, f}{\mathrm{d}\, x^p},
\frac{\mathrm{d}^n\, f}{\mathrm{d}\, x^n}(x_m)
- \frac{1}{h^n} \sum_{k \in \text{offsets}} a_k f_{m + k}
= c \frac{\mathrm{d}^p\, f}{\mathrm{d}\, x^p}(x_m),
where :math:`c` is the expected truncation error coefficient and :math:`p`
is the order of the approximation.
Expand All @@ -96,10 +96,10 @@ def determine_stencil_truncation_error(

c = 0.0
i = s.derivative
indices = s.indices.astype(s.coeffs.dtype)
offsets = s.offsets.astype(s.coeffs.dtype)
while i < 64 and np.allclose(c, 0.0, atol=atol, rtol=0.0):
i += 1
c = s.coeffs @ indices**i / math.factorial(i)
c = s.coeffs @ offsets**i / math.factorial(i)

return Truncation(i - s.derivative, c)

Expand All @@ -112,7 +112,7 @@ def modified_wavenumber(s: Stencil, k: Array) -> Array:
km = np.empty(k.shape, dtype=np.complex128)

for n in range(k.shape[0]):
exp = np.exp(1.0j * s.indices * k[n])
exp = np.exp(1.0j * s.offsets * k[n])
km[n] = np.sum(s.coeffs * exp)

return km
Expand Down Expand Up @@ -151,11 +151,11 @@ def make_taylor_approximation(
dtype = np.dtype(dtype)

# construct the system
indices = np.arange(bounds[0], bounds[1] + 1)
offsets = np.arange(bounds[0], bounds[1] + 1)
A = np.array(
[indices**i / math.factorial(i) for i in range(indices.size)], dtype=dtype
[offsets**i / math.factorial(i) for i in range(offsets.size)], dtype=dtype
)
b = np.zeros(indices.shape, dtype=dtype)
b = np.zeros(offsets.shape, dtype=dtype)
b[derivative] = 1.0

# determine coefficients
Expand All @@ -165,5 +165,5 @@ def make_taylor_approximation(
return Stencil(
derivative=derivative,
coeffs=x,
indices=indices,
offsets=offsets,
)
6 changes: 3 additions & 3 deletions tests/test_diff.py
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@
def finite_difference_convergence(d: Stencil) -> EOCRecorder:
eoc = EOCRecorder()

s = np.s_[abs(d.indices[0]) + 1 : -abs(d.indices[-1]) - 1]
s = np.s_[abs(d.offsets[0]) + 1 : -abs(d.offsets[-1]) - 1]
for n in [32, 64, 128, 256, 512]:
theta = np.linspace(0.0, 2.0 * np.pi, n, dtype=d.coeffs.dtype)
h = theta[1] - theta[0]
Expand Down Expand Up @@ -132,8 +132,8 @@ def test_finite_difference_taylor_stencil(*, visualize: bool = False) -> None:
a = np.array(
[-0.02651995, 0.18941314, -0.79926643, 0.0, 0.79926643, -0.18941314, 0.02651995]
)
indices = np.arange(-3, 4)
s = Stencil(derivative=1, coeffs=a, indices=indices)
offsets = np.arange(-3, 4)
s = Stencil(derivative=1, coeffs=a, offsets=offsets)

order, c = determine_stencil_truncation_error(s, atol=1.0e-6)
assert order == 4
Expand Down

0 comments on commit 1fbe547

Please sign in to comment.