Source code for pylops.basicoperators.SecondDerivative

import numpy as np
from numpy.core.multiarray import normalize_axis_index

from pylops import LinearOperator
from pylops.utils._internal import _value_or_list_like_to_tuple
from pylops.utils.backend import get_array_module
from pylops.utils.decorators import reshaped


[docs]class SecondDerivative(LinearOperator): r"""Second derivative. Apply a second derivative using a three-point stencil finite-difference approximation along ``axis``. Parameters ---------- dims : :obj:`list` or :obj:`int` Number of samples for each dimension (``None`` if only one dimension is available) axis : :obj:`int`, optional .. versionadded:: 2.0.0 Axis along which derivative is applied. sampling : :obj:`float`, optional Sampling step :math:`\Delta x`. kind : :obj:`str`, optional Derivative kind (``forward``, ``centered``, or ``backward``). edge : :obj:`bool`, optional Use shifted derivatives at edges (``True``) or ignore them (``False``). This is currently only available for centered derivative dtype : :obj:`str`, optional Type of elements in input array. name : :obj:`str`, optional .. versionadded:: 2.0.0 Name of operator (to be used by :func:`pylops.utils.describe.describe`) Attributes ---------- shape : :obj:`tuple` Operator shape explicit : :obj:`bool` Operator contains a matrix that can be solved explicitly (``True``) or not (``False``) Notes ----- The SecondDerivative operator applies a second derivative to any chosen direction of a multi-dimensional array. For simplicity, given a one dimensional array, the second-order centered first derivative is: .. math:: y[i] = (x[i+1] - 2x[i] + x[i-1]) / \Delta x^2 while the second-order forward stencil is: .. math:: y[i] = (x[i+2] - 2x[i+1] + x[i]) / \Delta x^2 and the second-order backward stencil is: .. math:: y[i] = (x[i] - 2x[i-1] + x[i-2]) / \Delta x^2 """ def __init__( self, dims, axis=-1, sampling=1, kind="centered", edge=False, dtype="float64", name="S", ): dims = _value_or_list_like_to_tuple(dims) super().__init__(dtype=np.dtype(dtype), dims=dims, dimsd=dims, name=name) self.axis = normalize_axis_index(axis, len(self.dims)) self.sampling = sampling self.kind = kind self.edge = edge self._register_multiplications(self.kind) def _register_multiplications(self, kind): # choose _matvec and _rmatvec kind if kind == "forward": self._matvec = self._matvec_forward self._rmatvec = self._rmatvec_forward elif kind == "centered": self._matvec = self._matvec_centered self._rmatvec = self._rmatvec_centered elif kind == "backward": self._matvec = self._matvec_backward self._rmatvec = self._rmatvec_backward else: raise NotImplementedError( "'kind' must be 'forward', 'centered' or 'backward'" ) @reshaped(swapaxis=True) def _matvec_forward(self, x): ncp = get_array_module(x) y = ncp.zeros(x.shape, self.dtype) y[..., :-2] = x[..., 2:] - 2 * x[..., 1:-1] + x[..., :-2] y /= self.sampling**2 return y @reshaped(swapaxis=True) def _rmatvec_forward(self, x): ncp = get_array_module(x) y = ncp.zeros(x.shape, self.dtype) y[..., :-2] += x[..., :-2] y[..., 1:-1] -= 2 * x[..., :-2] y[..., 2:] += x[..., :-2] y /= self.sampling**2 return y @reshaped(swapaxis=True) def _matvec_centered(self, x): ncp = get_array_module(x) y = ncp.zeros(x.shape, self.dtype) y[..., 1:-1] = x[..., 2:] - 2 * x[..., 1:-1] + x[..., :-2] if self.edge: y[..., 0] = x[..., 0] - 2 * x[..., 1] + x[..., 2] y[..., -1] = x[..., -3] - 2 * x[..., -2] + x[..., -1] y /= self.sampling**2 return y @reshaped(swapaxis=True) def _rmatvec_centered(self, x): ncp = get_array_module(x) y = ncp.zeros(x.shape, self.dtype) y[..., :-2] += x[..., 1:-1] y[..., 1:-1] -= 2 * x[..., 1:-1] y[..., 2:] += x[..., 1:-1] if self.edge: y[..., 0] += x[..., 0] y[..., 1] -= 2 * x[..., 0] y[..., 2] += x[..., 0] y[..., -3] += x[..., -1] y[..., -2] -= 2 * x[..., -1] y[..., -1] += x[..., -1] y /= self.sampling**2 return y @reshaped(swapaxis=True) def _matvec_backward(self, x): ncp = get_array_module(x) y = ncp.zeros(x.shape, self.dtype) y[..., 2:] = x[..., 2:] - 2 * x[..., 1:-1] + x[..., :-2] y /= self.sampling**2 return y @reshaped(swapaxis=True) def _rmatvec_backward(self, x): ncp = get_array_module(x) y = ncp.zeros(x.shape, self.dtype) y[..., :-2] += x[..., 2:] y[..., 1:-1] -= 2 * x[..., 2:] y[..., 2:] += x[..., 2:] y /= self.sampling**2 return y