pysdic.assemble_central_finite_difference_matrix#

assemble_central_finite_difference_matrix(order, n_times, n_dim=1, spacing=1.0, accuracy=2, sparse=False, weights=None)[source]#

Assemble the temporal derivation operator matrix for finite difference approximation with central scheme.

The operator matrix is used to approximate the temporal derivative of a time-dependent function using central finite differences.

For a time series with \(N_t\) time steps, the operator matrix \(D\) for the \(n\)-th order derivative is given by the toeplitz matrix:

\[\begin{split}D = \begin{bmatrix} a_{t} & a_{t+1} & 0 & \cdots & 0 & 0 \\ a_{t-1} & a_{t} & a_{t+1} & \cdots & 0 & 0 \\ 0 & a_{t-1} & a_{t} & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & a_{t-1} & a_{t} \\ \end{bmatrix}\end{split}\]

where \(a_k\) are the Stencil coefficients from the central finite difference kernel at the given order, spacing and accuracy.

The dimension n_dim allows to apply the operator to multi-dimensional data at each time step \(f = [f_1(t), f_2(t), \ldots, f_{n_{dim}}(t)]^T\).

Thus the kronecker product with the identity matrix is performed to account for the spatial dimensions.

\[\begin{split}D = T \otimes I_{n_{dim}} = \begin{bmatrix} a_{t} I_{n_{dim}} & a_{t+1} I_{n_{dim}} & 0 & \cdots & 0 & 0 \\ a_{t-1} I_{n_{dim}} & a_{t} I_{n_{dim}} & a_{t+1} I_{n_{dim}} & \cdots & 0 & 0 \\ 0 & a_{t-1} I_{n_{dim}} & a_{t} I_{n_{dim}} & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & a_{t-1} I_{n_{dim}} & a_{t} I_{n_{dim}} \\ \end{bmatrix}\end{split}\]

The operator matrix is then constructed such as \(\frac{\partial^n f}{\partial t^n} \approx D @ F\), where \(F\) is the flattened time series data \(F = [f_1(t_1), f_2(t_1), \ldots, f_{n_{dim}}(t_1), f_1(t_2), f_2(t_2), \ldots, f_{n_{dim}}(t_2), \ldots, f_{n_{dim}}(t_{N_t})]^T\).

Note

If the time series has less time steps than the number of stencil points required for the given order and accuracy, the operator matrix will be truncated accordingly.

Parameters:
  • order (Integral) – The order of the temporal derivative to approximate (e.g., 1 for first derivative, 2 for second derivative).

  • n_times (Integral) – The number of time steps in the time series.

  • n_dim (Optional[Integral], optional) – The number of spatial dimensions of the data at each time step. Default is 1.

  • spacing (Number, optional) – The time step size. Default is 1.0.

  • accuracy (Integral, optional) – The desired accuracy order of the approximation. Must be an even integer. Default is 2.

  • sparse (bool, optional) – Whether to return a sparse matrix. Default is False.

  • weights (Optional[ArrayLike], optional) – An optional array of weights to apply for the identity matrix in the kronecker product. If None, the identity matrix is used. If provided, it must have shape (n_dim,) to weight each spatial dimension accordingly (diagonal coefficients) or (n_dim, n_dim) to provide a full weighting matrix. Default is None.

Returns:

The temporal derivation operator matrix with shape (n_times * n_dim, n_times * n_dim).

Return type:

numpy.ndarray or scipy.sparse.csr_matrix

Raises:

ValueError – If the order or accuracy are not positive integers. If the accuracy is not an even integer. If the spacing is not a strictly positive number.

See also

compute_central_finite_difference_coefficients()

To compute the central finite difference coefficients for given order, spacing and accuracy.

assemble_forward_finite_difference_matrix()

To assemble the forward finite difference operator matrix.

assemble_backward_finite_difference_matrix()

To assemble the backward finite difference operator matrix.

Examples

>>> assemble_central_finite_difference_matrix(2, 5, spacing=1.0, accuracy=2)
array([[ -2.,  1.,  0.,  0.,  0.],
       [  1., -2.,  1.,  0.,  0.],
       [  0.,  1., -2.,  1.,  0.],
       [  0.,  0.,  1., -2.,  1.],
       [  0.,  0.,  0.,  1., -2.]])
>>> assemble_central_finite_difference_matrix(1, 4, n_dim=2, spacing=0.1, accuracy=4)
array([[ -500/2,      0,   400/3,      0, -100/12,       0,       0,       0],
       [      0, -500/2,       0,  400/3,       0, -100/12,       0,       0],
       [  400/3,      0,  -500/2,      0,   400/3,       0, -100/12,       0],
       [      0,  400/3,       0, -500/2,       0,   400/3,       0, -100/12],
       [-100/12,      0,   400/3,      0,  -500/2,       0,   400/3,       0],
       [      0,-100/12,       0,  400/3,       0,  -500/2,       0,   400/3],
       [      0,      0, -100/12,      0,   400/3,       0,  -500/2,       0],
       [      0,      0,       0,-100/12,      0,   400/3,       0,  -500/2]])