o  h@sddlZddlmZddlmZmZmZdgZGdddeZ GdddeZ Gdd d eZ Gd d d eZ Gd d d Z dS)N)LinearOperator)kroneye dia_array LaplacianNdcseZdZdZdejdfdd ZddZdd d Zd d Z d dZ dddZ ddZ ddZ ddZddZddZddZZS)rai" The grid Laplacian in ``N`` dimensions and its eigenvalues/eigenvectors. Construct Laplacian on a uniform rectangular grid in `N` dimensions and output its eigenvalues and eigenvectors. The Laplacian ``L`` is square, negative definite, real symmetric array with signed integer entries and zeros otherwise. Parameters ---------- grid_shape : tuple A tuple of integers of length ``N`` (corresponding to the dimension of the Lapacian), where each entry gives the size of that dimension. The Laplacian matrix is square of the size ``np.prod(grid_shape)``. boundary_conditions : {'neumann', 'dirichlet', 'periodic'}, optional The type of the boundary conditions on the boundaries of the grid. Valid values are ``'dirichlet'`` or ``'neumann'``(default) or ``'periodic'``. dtype : dtype Numerical type of the array. Default is ``np.int8``. Methods ------- toarray() Construct a dense array from Laplacian data tosparse() Construct a sparse array from Laplacian data eigenvalues(m=None) Construct a 1D array of `m` largest (smallest in absolute value) eigenvalues of the Laplacian matrix in ascending order. eigenvectors(m=None): Construct the array with columns made of `m` eigenvectors (``float``) of the ``Nd`` Laplacian corresponding to the `m` ordered eigenvalues. .. versionadded:: 1.12.0 Notes ----- Compared to the MATLAB/Octave implementation [1] of 1-, 2-, and 3-D Laplacian, this code allows the arbitrary N-D case and the matrix-free callable option, but is currently limited to pure Dirichlet, Neumann or Periodic boundary conditions only. The Laplacian matrix of a graph (`scipy.sparse.csgraph.laplacian`) of a rectangular grid corresponds to the negative Laplacian with the Neumann conditions, i.e., ``boundary_conditions = 'neumann'``. All eigenvalues and eigenvectors of the discrete Laplacian operator for an ``N``-dimensional regular grid of shape `grid_shape` with the grid step size ``h=1`` are analytically known [2]. References ---------- .. [1] https://github.com/lobpcg/blopex/blob/master/blopex_tools/matlab/laplacian/laplacian.m .. [2] "Eigenvalues and eigenvectors of the second derivative", Wikipedia https://en.wikipedia.org/wiki/Eigenvalues_and_eigenvectors_of_the_second_derivative Examples -------- >>> import numpy as np >>> from scipy.sparse.linalg import LaplacianNd >>> from scipy.sparse import diags, csgraph >>> from scipy.linalg import eigvalsh The one-dimensional Laplacian demonstrated below for pure Neumann boundary conditions on a regular grid with ``n=6`` grid points is exactly the negative graph Laplacian for the undirected linear graph with ``n`` vertices using the sparse adjacency matrix ``G`` represented by the famous tri-diagonal matrix: >>> n = 6 >>> G = diags(np.ones(n - 1), 1, format='csr') >>> Lf = csgraph.laplacian(G, symmetrized=True, form='function') >>> grid_shape = (n, ) >>> lap = LaplacianNd(grid_shape, boundary_conditions='neumann') >>> np.array_equal(lap.matmat(np.eye(n)), -Lf(np.eye(n))) True Since all matrix entries of the Laplacian are integers, ``'int8'`` is the default dtype for storing matrix representations. >>> lap.tosparse() >>> lap.toarray() array([[-1, 1, 0, 0, 0, 0], [ 1, -2, 1, 0, 0, 0], [ 0, 1, -2, 1, 0, 0], [ 0, 0, 1, -2, 1, 0], [ 0, 0, 0, 1, -2, 1], [ 0, 0, 0, 0, 1, -1]], dtype=int8) >>> np.array_equal(lap.matmat(np.eye(n)), lap.toarray()) True >>> np.array_equal(lap.tosparse().toarray(), lap.toarray()) True Any number of extreme eigenvalues and/or eigenvectors can be computed. >>> lap = LaplacianNd(grid_shape, boundary_conditions='periodic') >>> lap.eigenvalues() array([-4., -3., -3., -1., -1., 0.]) >>> lap.eigenvalues()[-2:] array([-1., 0.]) >>> lap.eigenvalues(2) array([-1., 0.]) >>> lap.eigenvectors(1) array([[0.40824829], [0.40824829], [0.40824829], [0.40824829], [0.40824829], [0.40824829]]) >>> lap.eigenvectors(2) array([[ 0.5 , 0.40824829], [ 0. , 0.40824829], [-0.5 , 0.40824829], [-0.5 , 0.40824829], [ 0. , 0.40824829], [ 0.5 , 0.40824829]]) >>> lap.eigenvectors() array([[ 0.40824829, 0.28867513, 0.28867513, 0.5 , 0.5 , 0.40824829], [-0.40824829, -0.57735027, -0.57735027, 0. , 0. , 0.40824829], [ 0.40824829, 0.28867513, 0.28867513, -0.5 , -0.5 , 0.40824829], [-0.40824829, 0.28867513, 0.28867513, -0.5 , -0.5 , 0.40824829], [ 0.40824829, -0.57735027, -0.57735027, 0. , 0. , 0.40824829], [-0.40824829, 0.28867513, 0.28867513, 0.5 , 0.5 , 0.40824829]]) The two-dimensional Laplacian is illustrated on a regular grid with ``grid_shape = (2, 3)`` points in each dimension. >>> grid_shape = (2, 3) >>> n = np.prod(grid_shape) Numeration of grid points is as follows: >>> np.arange(n).reshape(grid_shape + (-1,)) array([[[0], [1], [2]], [[3], [4], [5]]]) Each of the boundary conditions ``'dirichlet'``, ``'periodic'``, and ``'neumann'`` is illustrated separately; with ``'dirichlet'`` >>> lap = LaplacianNd(grid_shape, boundary_conditions='dirichlet') >>> lap.tosparse() >>> lap.toarray() array([[-4, 1, 0, 1, 0, 0], [ 1, -4, 1, 0, 1, 0], [ 0, 1, -4, 0, 0, 1], [ 1, 0, 0, -4, 1, 0], [ 0, 1, 0, 1, -4, 1], [ 0, 0, 1, 0, 1, -4]], dtype=int8) >>> np.array_equal(lap.matmat(np.eye(n)), lap.toarray()) True >>> np.array_equal(lap.tosparse().toarray(), lap.toarray()) True >>> lap.eigenvalues() array([-6.41421356, -5. , -4.41421356, -3.58578644, -3. , -1.58578644]) >>> eigvals = eigvalsh(lap.toarray().astype(np.float64)) >>> np.allclose(lap.eigenvalues(), eigvals) True >>> np.allclose(lap.toarray() @ lap.eigenvectors(), ... lap.eigenvectors() @ np.diag(lap.eigenvalues())) True with ``'periodic'`` >>> lap = LaplacianNd(grid_shape, boundary_conditions='periodic') >>> lap.tosparse() >>> lap.toarray() array([[-4, 1, 1, 2, 0, 0], [ 1, -4, 1, 0, 2, 0], [ 1, 1, -4, 0, 0, 2], [ 2, 0, 0, -4, 1, 1], [ 0, 2, 0, 1, -4, 1], [ 0, 0, 2, 1, 1, -4]], dtype=int8) >>> np.array_equal(lap.matmat(np.eye(n)), lap.toarray()) True >>> np.array_equal(lap.tosparse().toarray(), lap.toarray()) True >>> lap.eigenvalues() array([-7., -7., -4., -3., -3., 0.]) >>> eigvals = eigvalsh(lap.toarray().astype(np.float64)) >>> np.allclose(lap.eigenvalues(), eigvals) True >>> np.allclose(lap.toarray() @ lap.eigenvectors(), ... lap.eigenvectors() @ np.diag(lap.eigenvalues())) True and with ``'neumann'`` >>> lap = LaplacianNd(grid_shape, boundary_conditions='neumann') >>> lap.tosparse() >>> lap.toarray() array([[-2, 1, 0, 1, 0, 0], [ 1, -3, 1, 0, 1, 0], [ 0, 1, -2, 0, 0, 1], [ 1, 0, 0, -2, 1, 0], [ 0, 1, 0, 1, -3, 1], [ 0, 0, 1, 0, 1, -2]], dtype=int8) >>> np.array_equal(lap.matmat(np.eye(n)), lap.toarray()) True >>> np.array_equal(lap.tosparse().toarray(), lap.toarray()) True >>> lap.eigenvalues() array([-5., -3., -3., -2., -1., 0.]) >>> eigvals = eigvalsh(lap.toarray().astype(np.float64)) >>> np.allclose(lap.eigenvalues(), eigvals) True >>> np.allclose(lap.toarray() @ lap.eigenvectors(), ... lap.eigenvectors() @ np.diag(lap.eigenvalues())) True neumann)boundary_conditionsdtypecsF|dvr td|d||_||_t|}tj|||fddS)N) dirichletrperiodiczUnknown value zv is given for 'boundary_conditions' parameter. The valid options are 'dirichlet', 'periodic', and 'neumann' (default).)r shape) ValueError grid_shapernpprodsuper__init__)selfrrr N __class__^/var/www/vscode/kcb/lib/python3.10/site-packages/scipy/sparse/linalg/_special_sparse_arrays.pyrs  zLaplacianNd.__init__c Cs@|j}|durt|}t|}nt|tt||}t|}t|}t||D]O\}}|jdkrM|dt tj |dd|dd7}q-|jdkre|dt tj |d|d7}q-|dt tj t |dd|d7}q-| }t |} || } |dur| | d} | | d} | | fS)zCompute `m` largest eigenvalues in each of the ``N`` directions, i.e., up to ``m * N`` total, order them and return `m` largest. Nr r)rrindiceszerosmintuple ones_likeziprsinpifloorravelargsort) rmrrLeiggrid_shape_minjn Leig_ravelind eigenvaluesrrr_eigenvalue_orderings,     . &0 z LaplacianNd._eigenvalue_orderingNcCs||\}}|S)aReturn the requested number of eigenvalues. Parameters ---------- m : int, optional The positive number of smallest eigenvalues to return. If not provided, then all eigenvalues will be returned. Returns ------- eigenvalues : float array The requested `m` smallest or all eigenvalues, in ascending order. )r/)rr'r._rrrr.%szLaplacianNd.eigenvaluescCs\|jdkr&tjt|d|d}td|dt||d}nx|jdkrLtjt|d|}t|dkr?dnd|t||}nR|dkr]td|t|}nA|d|kr{|ddkr{td|tdd g|d}n#dtjt|d|}td|t|t |dd}d |t |t tj j k<|S) zjReturn 1 eigenvector in 1d with index `j` and number of grid points `n` where ``j < n``. r rg@?r?rrg)rrr#arangesqrtr"cosonestiler$absfinfofloat64eps)rr*r+ievrrr_ev1d6s & *$*zLaplacianNd._ev1dcsRfddt|jD}|d}|ddD] }tj||dd}qt|S)z{Return 1 eigenvector in Nd with multi-index `j` as a tensor product of the corresponding 1d eigenvectors. csg|] \}}||qSr)r?).0r*r+rrr Qsz(LaplacianNd._one_eve..rrN)axes)r!rr tensordotasarrayr%)rkphiresultrrAr_one_eveMs zLaplacianNd._one_evecst|\}}|durj}ntjttj|}t||}ddt|D}fdd|D}t|S)aReturn the requested number of eigenvectors for ordered eigenvalues. Parameters ---------- m : int, optional The positive number of eigenvectors to return. If not provided, then all eigenvectors will be returned. Returns ------- eigenvectors : float array An array with columns made of the requested `m` or all eigenvectors. The columns are ordered according to the `m` ordered eigenvalues. NcSsg|]}t|qSr)rr@xrrrrBnsz,LaplacianNd.eigenvectors..csg|]}|qSr)rI)r@rFrArrrBos) r/rrrrr unravel_indexr! column_stack)rr'r0r-r) N_indiceseigenvectors_listrrAr eigenvectorsWs  zLaplacianNd.eigenvectorsc Csh|j}t|}tj||gtjd}t|}t|}t|D] \}}d|dd<dtd|d|d|fdd<dtd|d|dd|fdd<dtd|d|d|dfdd<|jdkryd|d <d||d|df<n*|jd kr|dkr|d|dfd7<||ddfd7<n|d d7<|}|dkrt|d|} t d| D]$} |d|d|f|| || d|| || d|f<||7}q|d|d|f|d|d|f<t t||dd} d|d|d|f<d d t | D} |d|d|f| || || fdd| dd| f<||7}q | |j S) z Converts the Laplacian data to a dense array. Returns ------- L : ndarray The shape is ``(N, N)`` where ``N = np.prod(grid_shape)``. r rNzii->irrr3rrr cSsg|]}|qSrrrJrrrrBsz'LaplacianNd.toarray..)rrrrint8 empty_like enumerateeinsumrrangeintreshapeastyper ) rrr+LL_iLtempr-dimnew_dimtilesr*idxrrrtoarrayrsL    $((  < $  zLaplacianNd.toarrayc CsZt|j}t|j}t||ftjd}t|D]}|j|}tjd|gtjd}|dddfd9<|jdkrBd|d<d|d <t|gd f||ftjd }|jd krwt||ftjd}|j dg| dd |j dg|dd ||7}t|D]} t t |j| tjd|}q{t|d|D]} t |t |j| tjd}q||7}q| |j S)a) Constructs a sparse array from the Laplacian data. The returned sparse array format is dependent on the selected boundary conditions. Returns ------- L : scipy.sparse.sparray The shape is ``(N, N)`` where ``N = np.prod(grid_shape)``. rQrNrRrr3rr)rr3r3rrr r r )rF)lenrrrrrTrXr7rsetdiagrrr[r ) rrpr\r=r_datar]tr*rrrtosparses0        zLaplacianNd.tosparsec Cs(|j}t|}||d}d||}t|D]}|tj|d|d7}|tj|d|d7}|jdvr |tdf|dtdf||dtj|d|dtdf|dtdf||d8<|tdf|dtdf||dtj|d|dtdf|dtdf||d8<|jdkr |tdf|dtdf||dtj|d |dtdf|dtdf||d7<|tdf|dtdf||dtj|d |dtdf|dtdf||d7<q|d|jdS) N)r3rRr)axisr3)rr )rrr) rrhrZrXrrollrslicer )rrKrrXYr=rrr_matvecsJ   ,&&& &&&&zLaplacianNd._matveccC ||SNrsrrKrrr_matmats zLaplacianNd._matmatcC|SrurrArrr_adjointzLaplacianNd._adjointcCryrurrArrr _transposer{zLaplacianNd._transposeru)__name__ __module__ __qualname____doc__rrTrr/r.r?rIrPrcrmrsrxrzr| __classcell__rrrrr s"l  @)!csheZdZdZejffdd ZdddZddZd d Z d d Z d dZ ddZ ddZ ddZZS)Sakuraia Construct a Sakurai matrix in various formats and its eigenvalues. Constructs the "Sakurai" matrix motivated by reference [1]_: square real symmetric positive definite and 5-diagonal with the main diagonal ``[5, 6, 6, ..., 6, 6, 5], the ``+1`` and ``-1`` diagonals filled with ``-4``, and the ``+2`` and ``-2`` diagonals made of ``1``. Its eigenvalues are analytically known to be ``16. * np.power(np.cos(0.5 * k * np.pi / (n + 1)), 4)``. The matrix gets ill-conditioned with its size growing. It is useful for testing and benchmarking sparse eigenvalue solvers especially those taking advantage of its banded 5-diagonal structure. See the notes below for details. Parameters ---------- n : int The size of the matrix. dtype : dtype Numerical type of the array. Default is ``np.int8``. Methods ------- toarray() Construct a dense array from Laplacian data tosparse() Construct a sparse array from Laplacian data tobanded() The Sakurai matrix in the format for banded symmetric matrices, i.e., (3, n) ndarray with 3 upper diagonals placing the main diagonal at the bottom. eigenvalues All eigenvalues of the Sakurai matrix ordered ascending. Notes ----- Reference [1]_ introduces a generalized eigenproblem for the matrix pair `A` and `B` where `A` is the identity so we turn it into an eigenproblem just for the matrix `B` that this function outputs in various formats together with its eigenvalues. .. versionadded:: 1.12.0 References ---------- .. [1] T. Sakurai, H. Tadano, Y. Inadomi, and U. Nagashima, "A moment-based method for large-scale generalized eigenvalue problems", Appl. Num. Anal. Comp. Math. Vol. 1 No. 2 (2004). Examples -------- >>> import numpy as np >>> from scipy.sparse.linalg._special_sparse_arrays import Sakurai >>> from scipy.linalg import eig_banded >>> n = 6 >>> sak = Sakurai(n) Since all matrix entries are small integers, ``'int8'`` is the default dtype for storing matrix representations. >>> sak.toarray() array([[ 5, -4, 1, 0, 0, 0], [-4, 6, -4, 1, 0, 0], [ 1, -4, 6, -4, 1, 0], [ 0, 1, -4, 6, -4, 1], [ 0, 0, 1, -4, 6, -4], [ 0, 0, 0, 1, -4, 5]], dtype=int8) >>> sak.tobanded() array([[ 1, 1, 1, 1, 1, 1], [-4, -4, -4, -4, -4, -4], [ 5, 6, 6, 6, 6, 5]], dtype=int8) >>> sak.tosparse() >>> np.array_equal(sak.dot(np.eye(n)), sak.tosparse().toarray()) True >>> sak.eigenvalues() array([0.03922866, 0.56703972, 2.41789479, 5.97822974, 10.54287655, 14.45473055]) >>> sak.eigenvalues(2) array([0.03922866, 0.56703972]) The banded form can be used in scipy functions for banded matrices, e.g., >>> e = eig_banded(sak.tobanded(), eigvals_only=True) >>> np.allclose(sak.eigenvalues, e, atol= n * n * n * np.finfo(float).eps) True cs&||_||_||f}t||dSru)r+r rr)rr+r r rrrraszSakurai.__init__Nc CsZ|dur|j}t|jd||jd}tdttd|tj|jddS)aReturn the requested number of eigenvalues. Parameters ---------- m : int, optional The positive number of smallest eigenvalues to return. If not provided, then all eigenvalues will be returned. Returns ------- eigenvalues : `np.float64` array The requested `m` smallest or all eigenvalues, in ascending order. Nrg0@r2)r+rr4flippowerr6r#)rr'rFrrrr.gs0zSakurai.eigenvaluescCsftjddtj|jd|jddf}dtj|j|jd}tj|j|jd}t|||g|jS)zA Construct the Sakurai matrix as a banded array. rrQr)rr_r7r+r arrayr[)rd0d1d2rrrtobandedzs&zSakurai.tobandedcCsHddlm}|}||d|d|d|d|dggd|j|jS)zB Construct the Sakurai matrix is a sparse format. r)spdiagsrr)rRr3rrr) scipy.sparserrr+)rrdrrrrms (zSakurai.tosparsecC |SrurmrcrArrrrc zSakurai.toarraycCsD||jd}t|j|j}tj||d}d|dddfd|dddf|dddf|dddf<d|dddfd|d ddf|d ddf|dddf<d |ddddfd|dd ddf|ddddft|dd ddfd t|d dddfd|ddddf<|S)z Construct matrix-free callable banded-matrix-vector multiplication by the Sakurai matrix without constructing or storing the matrix itself using the knowledge of its entries and the 5-diagonal format. r3rQrrNrrrrRr)rerSrd))rrrS)rZr+r promote_typesr zeros_likepad)rrK result_dtypesxrrrrssDDBzSakurai._matveccCrt)z Construct matrix-free callable matrix-matrix multiplication by the Sakurai matrix without constructing or storing the matrix itself by reusing the ``_matvec(x)`` that supports both 1D and 2D arrays ``x``. rvrwrrrrx zSakurai._matmatcCryrurrArrrrzr{zSakurai._adjointcCryrurrArrrr|r{zSakurai._transposeru)r}r~rrrrTrr.rrmrcrsrxrzr|rrrrrrsZ   rcsfeZdZdZejffdd ZddZddZdd Z d d Z d d Z ddZ ddZ ddZZS)MikotaMas Construct a mass matrix in various formats of Mikota pair. The mass matrix `M` is square real diagonal positive definite with entries that are reciprocal to integers. Parameters ---------- shape : tuple of int The shape of the matrix. dtype : dtype Numerical type of the array. Default is ``np.float64``. Methods ------- toarray() Construct a dense array from Mikota data tosparse() Construct a sparse array from Mikota data tobanded() The format for banded symmetric matrices, i.e., (1, n) ndarray with the main diagonal. cs||_||_t||dSru)r r rr)rr r rrrrszMikotaM.__init__cCs"dtd|jdd|jS)Nr1rr)rr4r r[r rArrr_diags"z MikotaM._diagcCs|Sru)rrArrrrszMikotaM.tobandedcCs(ddlm}||gdg|j|jdS)Nrdiagsrg)rrrr r rrrrrrms zMikotaM.tosparsecCst||jSru)rdiagrr[r rArrrrcszMikotaM.toarraycCs,||jdd}|ddtjf|S)z Construct matrix-free callable banded-matrix-vector multiplication by the Mikota mass matrix without constructing or storing the matrix itself using the knowledge of its entries and the diagonal format. rr3N)rZr rrnewaxisrwrrrrsszMikotaM._matveccCrt)z Construct matrix-free callable matrix-matrix multiplication by the Mikota mass matrix without constructing or storing the matrix itself by reusing the ``_matvec(x)`` that supports both 1D and 2D arrays ``x``. rvrwrrrrxrzMikotaM._matmatcCryrurrArrrrzr{zMikotaM._adjointcCryrurrArrrr|r{zMikotaM._transpose)r}r~rrrr;rrrrmrcrsrxrzr|rrrrrrs rcs^eZdZdZejffdd ZddZddZdd Z d d Z d d Z ddZ ddZ ZS)MikotaKa Construct a stiffness matrix in various formats of Mikota pair. The stiffness matrix `K` is square real tri-diagonal symmetric positive definite with integer entries. Parameters ---------- shape : tuple of int The shape of the matrix. dtype : dtype Numerical type of the array. Default is ``np.int32``. Methods ------- toarray() Construct a dense array from Mikota data tosparse() Construct a sparse array from Mikota data tobanded() The format for banded symmetric matrices, i.e., (2, n) ndarray with 2 upper diagonals placing the main diagonal at the bottom. cs`||_||_t|||d}tjd|ddd|jd|_tj|ddd|jd |_dS)NrrrrRrQr3)r r rrrr4_diag0_diag1)rr r r+rrrr s  zMikotaK.__init__cCstt|jdd|jgS)Nreconstant)rrrrrrArrrrszMikotaK.tobandedcCs0ddlm}||j|j|jggd|j|jdS)Nrrrfrg)rrrrr r rrrrrms zMikotaK.tosparsecCrrurrArrrrcrzMikotaK.toarraycCs4||jdd}t|j|j}tj||d}|j}|j}|d|dddf|d|dddf|dddf<|d|dddf|d|dddf|dddf<|dddf|ddddf|dddf|ddddf|dddf|ddddf|ddddf<|S)z Construct matrix-free callable banded-matrix-vector multiplication by the Mikota stiffness matrix without constructing or storing the matrix itself using the knowledge of its entries and the 3-diagonal format. rr3rQNrrRr)rZr rrr rrr)rrKrkxrrrrrrs"s<<"""zMikotaK._matveccCrt)z Construct matrix-free callable matrix-matrix multiplication by the Stiffness mass matrix without constructing or storing the matrix itself by reusing the ``_matvec(x)`` that supports both 1D and 2D arrays ``x``. rvrwrrrrx4rzMikotaK._matmatcCryrurrArrrrz<r{zMikotaK._adjointcCryrurrArrrr|?r{zMikotaK._transpose)r}r~rrrint32rrrmrcrsrxrzr|rrrrrrs rc@s(eZdZdZejfddZdddZdS) MikotaPaira Construct the Mikota pair of matrices in various formats and eigenvalues of the generalized eigenproblem with them. The Mikota pair of matrices [1, 2]_ models a vibration problem of a linear mass-spring system with the ends attached where the stiffness of the springs and the masses increase along the system length such that vibration frequencies are subsequent integers 1, 2, ..., `n` where `n` is the number of the masses. Thus, eigenvalues of the generalized eigenvalue problem for the matrix pair `K` and `M` where `K` is the system stiffness matrix and `M` is the system mass matrix are the squares of the integers, i.e., 1, 4, 9, ..., ``n * n``. The stiffness matrix `K` is square real tri-diagonal symmetric positive definite. The mass matrix `M` is diagonal with diagonal entries 1, 1/2, 1/3, ...., ``1/n``. Both matrices get ill-conditioned with `n` growing. Parameters ---------- n : int The size of the matrices of the Mikota pair. dtype : dtype Numerical type of the array. Default is ``np.float64``. Attributes ---------- eigenvalues : 1D ndarray, ``np.uint64`` All eigenvalues of the Mikota pair ordered ascending. Methods ------- MikotaK() A `LinearOperator` custom object for the stiffness matrix. MikotaM() A `LinearOperator` custom object for the mass matrix. .. versionadded:: 1.12.0 References ---------- .. [1] J. Mikota, "Frequency tuning of chain structure multibody oscillators to place the natural frequencies at omega1 and N-1 integer multiples omega2,..., omegaN", Z. Angew. Math. Mech. 81 (2001), S2, S201-S202. Appl. Num. Anal. Comp. Math. Vol. 1 No. 2 (2004). .. [2] Peter C. Muller and Metin Gurgoze, "Natural frequencies of a multi-degree-of-freedom vibration system", Proc. Appl. Math. Mech. 6, 319-320 (2006). http://dx.doi.org/10.1002/pamm.200610141. Examples -------- >>> import numpy as np >>> from scipy.sparse.linalg._special_sparse_arrays import MikotaPair >>> n = 6 >>> mik = MikotaPair(n) >>> mik_k = mik.k >>> mik_m = mik.m >>> mik_k.toarray() array([[11., -5., 0., 0., 0., 0.], [-5., 9., -4., 0., 0., 0.], [ 0., -4., 7., -3., 0., 0.], [ 0., 0., -3., 5., -2., 0.], [ 0., 0., 0., -2., 3., -1.], [ 0., 0., 0., 0., -1., 1.]]) >>> mik_k.tobanded() array([[ 0., -5., -4., -3., -2., -1.], [11., 9., 7., 5., 3., 1.]]) >>> mik_m.tobanded() array([1. , 0.5 , 0.33333333, 0.25 , 0.2 , 0.16666667]) >>> mik_k.tosparse() >>> mik_m.tosparse() >>> np.array_equal(mik_k(np.eye(n)), mik_k.toarray()) True >>> np.array_equal(mik_m(np.eye(n)), mik_m.toarray()) True >>> mik.eigenvalues() array([ 1, 4, 9, 16, 25, 36]) >>> mik.eigenvalues(2) array([ 1, 4]) cCs:||_||_||f|_t|j|j|_t|j|j|_dSru)r+r r rr'rrF)rr+r rrrrs  zMikotaPair.__init__NcCs,|dur|j}tjd|dtjd}||S)aReturn the requested number of eigenvalues. Parameters ---------- m : int, optional The positive number of smallest eigenvalues to return. If not provided, then all eigenvalues will be returned. Returns ------- eigenvalues : `np.uint64` array The requested `m` smallest or all eigenvalues, in ascending order. NrrQ)r+rr4uint64)rr' arange_plus1rrrr.szMikotaPair.eigenvaluesru)r}r~rrrr;rr.rrrrrCsXr)numpyrscipy.sparse.linalgrrrrr__all__rrrrrrrrrs +DO