gddlZddlmZddlmZmZmZdgZGddeZ GddeZ Gdd eZ Gd d eZ Gd d Z y)N)LinearOperator)kroneye dia_array LaplacianNdceZdZdZdej dfd ZdZddZdZ dZ dd Z d Z d Z d Zd ZdZdZxZS)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_conditionsdtypec|dvrtd|d||_||_tj|}t ||||fy)N) dirichletr periodiczUnknown value zv is given for 'boundary_conditions' parameter. The valid options are 'dirichlet', 'periodic', and 'neumann' (default).)r shape) ValueError grid_shaper npprodsuper__init__)selfrr r N __class__s ]/var/www/html/venv/lib/python3.12/site-packages/scipy/sparse/linalg/_special_sparse_arrays.pyrzLaplacianNd.__init__se &J J !4 78DD  %#6 GGJ  uQF3c t|j}|+tj|}tj|}nUt |t tj ||z}tj|}tj|}t||D]\}}|jdk(r<|dtjtj|dzzd|dzzz dzzz }Q|jdk(r6|dtjtj|zd|zz dzzz }|dtjtjtj|dzdz z|z dzzz }|j}tj|} || } | | | 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_likezipr sinpifloorravelargsort) rmrrLeiggrid_shape_minjn Leig_ravelind eigenvaluess r_eigenvalue_orderingz LaplacianNd._eigenvalue_orderings__ 9jj,G88J'D !&r||J'?!'C!DFNjj0G88N+D, LDAq'';6RVVBEEQUOqAE{$CDIII))Y6RVVBEEAIQ$78A===RVVBEEBHHa!eq[,A$AA$EF!KKK  LZZ\ jj$ o =%qbc*Kqbc(CCrc.|j|\}}|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. )r2)rr*r1_s rr1zLaplacianNd.eigenvalues%s2215 QrcL|jdk(rhtjtj|dzz|dzz }tjd|dzz tj ||dzzz}nf|jdk(retjtj|dzz|z }tj|dk(rdnd|z tj ||zz}n|dk(r/tjd|z tj|z}n|dz|k(r=|dzdk(r5tjd|z tjdd g|dzz}nydtjztj|dzz|z }tjd|z tj |tj|dzdz zz}d |tj|tjtjjk<|S) zjReturn 1 eigenvector in 1d with index `j` and number of grid points `n` where ``j < n``. r rg@?r ?rrg)r rr&arangesqrtr%cosonestiler'absfinfofloat64eps)rr-r.ievs r_ev1dzLaplacianNd._ev1d6s  # #{ 21)*a!e4Aq2v'"&&a!e*==B  % % 21+,q0AQ"B!34rvva!e}DBAvWWR!V_rwwqz1Q!A WWR!V_rww2w1'==J"))A,"459WWR!V_rvva"((AEQ;2G.G'HH57266":,00 01 rct||jDcgc]\}}|j||}}}|d}|ddD]}tj||d}tj |j Scc}}w)z{Return 1 eigenvector in Nd with multi-index `j` as a tensor product of the corresponding 1d eigenvectors. rrN)axes)r$rrDr tensordotasarrayr()rkr-r.phiresults r_one_evezLaplacianNd._one_eveMs-04??,CDDAqtzz!QDDQqr7 7C\\&#A6F 7zz&!'')) EsBc|j|\}}| |j}n?t|jtt j |j|z}t j ||}t|Dcgc] }t|}}|Dcgc]}|j|}}t j|Scc}wcc}w)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. ) r2rr!r"rr# unravel_indexr$rL column_stack) rr*r4r0r, N_indicesxrIeigenvectors_lists r eigenvectorszLaplacianNd.eigenvectorsWs**1-3 9!__N  %bll4??&Ca&G HJN$$S.9 '*I7!U1X7 77@A!T]]1-AA0118As C Cc |j}tj|}tj||gtj}tj |}tj |}t |D]\}}d|dddtjd|d|d|fdddtjd|d|dz d|fdddtjd|d|d|dz fdd|jdk(rd|d <d||dz |dz f<nF|jd k(r7|dkDr%|d|dz fxxdz cc<||dz dfxxdz cc<n |d xxdz cc<|}|dkDrTtj|d|} td| D]-} |d|d|f|| |z| dz|z| |z| dz|zf<||z }/|d|d|f|d|d|f<ttj||dzd} d|d|d|f<t| D cgc]} | } } |d|d|f|j|| || fdd| dd| f<||z }|j|jScc} w) 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->irr r8rrr)rrrr int8 empty_like enumerateeinsumr rangeintreshapeastyper ) rrr.LL_iLtempr0dimnew_dimtilesr-rQidxs rtoarrayzLaplacianNd.toarrayrs__ GGJ  HHaV277 +mmA a !*-+ HCCF68BIIgs4C4#: / 2;(?E(7(HWH$ % 3q56 234E&'C(7(" ##El+1+C+%*(7(HWH*<$= KK%%! S!S. " HAW+ Zxx ##,s I c t|j}tj|j}t ||ftj }t |D]r}|j|}tjd|gtj }|dddfxxdzcc<|jdk(r d|d<d|d <t |gd f||ftj }|jd k(rQt ||ftj }|jdg| dz |jdg|dz ||z }t |D]4} tt|j| tj |}6t |dz|D]4} t|t|j| tj }6||z }u|j|jS)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)``. rUrNrVr r8rr)rr8r8rrrr r)rI)lenrrrrrXr\r<r setdiagrrr_r ) rrpr`rBrcdataratr-s rtosparsezLaplacianNd.tosparses   GGDOO $ q!fBGG ,q A//!$C77As82773D AJ" J''94T  U T:.sCj"$''C'':5sCj8 1##a ( 1#Q 'q1X H3tq1A3G H1q5!_ H3DOOA$6bgg FG H HA/ 0xx ##rc |j}t|}|j|dz}d|z|z}t|D]!}|t j |d|z }|t j |d|z }|j dvsI|tdf|zdztdf||z dz zzxxt j |d|tdf|zdztdf||z dz zzzcc<|tdf|zdztdf||z dz zzxxt j |d|tdf|zdztdf||z dz zzzcc<|j dk(s>|tdf|zdztdf||z dz zzxxt j |d |tdf|zdztdf||z dz zzz cc<|tdf|zdztdf||z dz zzxxt j |d |tdf|zdztdf||z dz zzz cc<$|jd|jdS) N)r8rVr)axisr8)r r )rr r) rrmr^r\rrollr slicer)rrQrrXYrBs r_matveczLaplacianNd._matvecs__  O IIj5( ) FQJq A AA& &A BQ' 'A''+CC5;."T)U4[NAaCE,BBwwq!!,4[NQ&-t!A#a%0HH4[NQ&.%+1Q3q51IIWWQ+4[NQ&.%+1Q3q51II ++y8t*T1U4[Nac!e4LLAA.t*T1U4[Nac!e4LL t*U2eDk^qs1u5MMAA.t*U2eDk^qs1u5MM) 4yyQWWR[))rc$|j|SNryrrQs r_matmatzLaplacianNd._matmats||Arc|Sr{rs r_adjointzLaplacianNd._adjoint rc|Sr{rrs r _transposezLaplacianNd._transposerrr{)__name__ __module__ __qualname____doc__rrXrr2r1rDrLrSrgrrryr~rr __classcell__rs@rrr sUhV&/ww4" >".*26>$@'$R*BrcleZdZdZej ffd Zd dZdZdZ dZ dZ dZ d Z d ZxZS) 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 cJ||_||_||f}t| ||yr{)r.r rr)rr.r rrs rrzSakurai.__init__as) A &rc J| |j}tj|jdz|z |jdz}tjdtjtj d|ztj z|jdzz dzS)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. rg0@r7)r.rr9flippowerr;r&)rr*rIs rr1zSakurai.eigenvaluesgsw 9A IIdffqj!mTVVaZ 0wwsRXXbffS1Wruu_ -K&LaPPQQrctjddtj|jdz |jzdf}dtj|j|jz}tj|j|j}tj |||gj |jS)zA Construct the Sakurai matrix as a banded array. rrUr)rr_r<r.r arrayr_)rd0d1d2s rtobandedzSakurai.tobandedzsUU1a"''$&&1*DJJ??B C "''$&& 3 3 WWTVV4:: .xxR %,,TZZ88rcddlm}|j}||d|d|d|d|dggd|j|jS)zB Construct the Sakurai matrix is a sparse format. r)spdiagsrr)rVr8rrr) scipy.sparserrr.)rrds rrrzSakurai.tosparsesR ) MMO!adAaD!A$!57Hvvtvv' 'rc>|jjSr{rrrgrs rrgzSakurai.toarray}}&&((rcL|j|jd}tj|j|j}tj ||}d|dddfzd|dddfzz |dddfz|dddf<d|dddfzd|d ddfzz |d ddfz|dddf<d |ddddfzd|dd ddf|ddddfzzz tj |dd ddfd ztj |d dddfdz|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. r8rUrrNrrrrVr)rjrWri))rrrW)r^r.r promote_typesr zeros_likepad)rrQ result_dtypesxs rryzSakurai._matvecs; IIdffb !''< ]]1L 1qAw;Qq!tW,qAw61a4"a%L1qQx</!BE(:2q5 AaeQhK!q"ay1QRU8/C*DDq"ay*:;<qQx)9:;1b5!8  rc$|j|S)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``. r|r}s rr~zSakurai._matmat ||Arc|Sr{rrs rrzSakurai._adjointrrc|Sr{rrs rrzSakurai._transposerrr{)rrrrrrXrr1rrrrgryr~rrrrs@rrrsAYt!#' R&9 ') rrcjeZdZdZej ffd ZdZdZdZ dZ dZ dZ d Z d ZxZS) 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. cB||_||_t| ||yr{)rr rr)rrr rs rrzMikotaM.__init__s    &rcdtjd|jddzz j|jS)Nr6rr)rr9rr_r rs r_diagz MikotaM._diags6RYYq$**Q-!"344<|jjSr{rrs rrgzMikotaK.toarrayrrc"|j|jdd}tj|j|j}tj ||}|j }|j}|d|dddfz|d|dddfzz|dddf<|d|dddfz|d|dddfzz|dddf<|dddf|ddddfz|dddf|ddddfzz|dddf|ddddfzz|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. rr8rUNrrVr)r^rrrr rrr)rrQrkxrrs rryzMikotaK._matvec"s9 IIdjjmR (''< ]]1L 1 [[ [[a51QT7?RUQq!tW_41a4rFQr1uX%22q5(992q5 3B39 $B$' 2QUD[/AaeQhK78QRX,12q5121b5!8  rc$|j|S)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``. r|r}s rr~zMikotaK._matmat4rrc|Sr{rrs rrzMikotaK._adjoint<rrc|Sr{rrs rrzMikotaK._transpose?rr)rrrrrint32rrrrrgryr~rrrrs@rrrs;0%'HHBP9 )$rrc6eZdZdZej fdZddZy) 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]) c||_||_||f|_t|j|j|_t |j|j|_yr{)r.r rrr*rrI)rr.r s rrzMikotaPair.__init__sG V TZZ0TZZ0rNcz| |j}tjd|dztj}||zS)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. rrU)r.rr9uint64)rr* arange_plus1s rr1zMikotaPair.eigenvaluess7 9AyyAE; l**rr{)rrrrrr@rr1rrrrrCsWp!# 1+rr)numpyrscipy.sparse.linalgrrrrr__all__rrrrrrrrrs`.-- / y.yxgngTAnAHLnL^q+q+r