&jddlZddlmZddlmZmZmZdgZGddeZ GddeZ Gdd eZ Gd d eZ Gd d Z dS)N)LinearOperator)kron eye_array dia_array LaplacianNdczeZdZdZdejdfd ZdZddZdZ d Z dd Z d Z d Z d ZdZdZdZxZS)ray" 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_array, 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_array(np.ones(n - 1), offsets=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.28867513, 0.40824829], [ 0.57735027, 0.40824829], [ 0.28867513, 0.40824829], [-0.28867513, 0.40824829], [-0.57735027, 0.40824829], [-0.28867513, 0.40824829]]) >>> lap.eigenvectors() array([[ 0.40824829, 0.5 , 0.28867513, 0.28867513, 0.5 , 0.40824829], [-0.40824829, 0. , -0.57735027, 0.57735027, 0. , 0.40824829], [ 0.40824829, -0.5 , 0.28867513, 0.28867513, -0.5 , 0.40824829], [-0.40824829, 0.5 , 0.28867513, -0.28867513, -0.5 , 0.40824829], [ 0.40824829, 0. , -0.57735027, -0.57735027, 0. , 0.40824829], [-0.40824829, -0.5 , 0.28867513, -0.28867513, 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 |||fdS)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 XC:\Projects\DSCons\.venv\Lib\site-packages\scipy\sparse\linalg\_special_sparse_arrays.pyrzLaplacianNd.__init__s &J J JD!4DDD  %#6 GJ   uQF33333c d|j}|)tj|}tj|}nZt |t tj||z}tj|}tj|}t||D]\}}|jdkr7|dtj tj |dzzd|dzzz dzzz }G|jdkr1|dtj tj |zd|zz dzzz }|dtj tj tj |dzdz z|z dzzz }| }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_ 9j,,G8J''DD !&r|J'?'?!'C!D!DFFNj00G8N++D,, L LDAq';66RVBEQUOqAE{$CDDIII)Y66RVBEAIQ$788A===RVBEBHa!eq[,A,A$AA$EFF!KKKZZ\\ j$$ o =%qbcc*Kqbcc(CCrNc6||\}}|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%s!22155 Qrc|jdkratjtj|dzz|dzz }tjd|dzz tj||dzzz}n|jdkr`tjtj|dzz|z }tj|dkrdnd|z tj||zz}n5|dkr.tjd|z tj|z}n|dz|kr<|dzdkr3tjd|z tjdd g|dzz}n|dzdzdkrZtjtj|dzz|z }tjd|z tj||dzzz}nVtjtj|dzz|z }tjd|z tj||zz}d |tj |tj tj j k<|S) ziReturn 1 eigenvector in 1d with index `j` and number of grid points `n` where ``j < n``. r rg@?r ?rrg) r rr&arangesqrtr%cosonestileabsfinfofloat64eps)rr-r.ievs r_ev1dzLaplacianNd._ev1d6s  #{ 2 21)*a!e4Aq2v''"&a!e*=*==BB  % 2 21+,q0AQ""B!344rva!e}}DBBAvvWR!V__rwqzz1Q!A WR!V__rw2w1'='==a%1!!ERYq\\C/014WR!V__rva1q5k':'::ERYq\\C/014WR!V__rva!e}}45726"::,,0 01 rcfdt|jD}|d}|ddD]}tj||d}tj|S)zzReturn 1 eigenvector in Nd with multi-index `j` as a tensor product of the corresponding 1d eigenvectors. cBg|]\}}||S)rD).0r-r.rs r z(LaplacianNd._one_eve..Ts+DDDDAqtzz!QDDDrrrN)axes)r$rr tensordotasarrayr()rkphiresults` r_one_evezLaplacianNd._one_evePsEDDDC4?,C,CDDDQqrr7 7 7C\&#A666FFz&!!'')))rcT|\}}|j}n.qs777!U1XX777rc:g|]}|SrG)rP)rHrMrs rrIz,LaplacianNd.eigenvectors..rs%AAA!T]]1--AAAr) r2rr!r"rr# unravel_indexr$ column_stack)rr*r4r0r, N_indiceseigenvectors_lists` r eigenvectorszLaplacianNd.eigenvectorsZs**1--3 9!_NN  %bl4?&C&Ca&G H HJJN$S.99 77sI777 AAAAyAAA0111rc |j}tj|}tj||gtj}tj|}tj|}t |D]\}}d|dd<dtjd|d|d|fdd<dtjd|d|dz d|fdd<dtjd|d|d|dz fdd<|jdkrd|d <d||dz |dz f<nL|jd krA|dkr+|d|dz fxxdz cc<||dz dfxxdz cc<n|d xxdz cc<|}|dkr]tj|d|} td| D]0} |d|d|f|| |z| dz|z| |z| dz|zf<||z }1|d|d|f|d|d|f<ttj||dzd} d|d|d|f<d t| D} |d|d|f| || || fdd| dd| f<||z }| |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->irr r8rrrcg|]}|SrGrGrSs rrIz'LaplacianNd.toarray..s+++1+++r)rrrr int8 empty_like enumerateeinsumr rangeintreshapeastyper ) rrr.LL_iLtempr0dimnew_dimtilesr-idxs rtoarrayzLaplacianNd.toarrayus_ GJ   HaV27 + + +mA a  !*--+ + HCCF68BIgs4C4#: / / 2;(?E(7(HWH$ % 3q566 23344E&'C(7(" #++eEll+++C%*(7(HWH*<$= KK%%!  !!S!!!S. " HAAxx ###rc t|j}tj|j}t ||ftj}t |D]v}|j|}tjd|gtj}|dddfxxdzcc<|jdkr d|d<d|d <t |gd f||ftj }|jd krYt ||ftj}| dg| dz | dg|dz ||z }t |D]6} tt|j| tj|}7t |dz|D]6} t|t|j| tj}7||z }x| |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)``. r\rNr]r r8rr)rr8r8rr)rr r)rM)lenrrrrr`rdr<r setdiagrrrgr ) rrprhrBrkdataritr-s rtosparsezLaplacianNd.tosparses   GDO $ $ q!fBG , , ,q  A/!$C7As827333D AAAJJJ" JJJ'944T  U T:::.sCj"$'C':55sCj888 1##a ((( 1#Q '''q1XX N N9T_Q%7rwGGGMM1q5!__ N N3 $/!*"4***22226>$>$>$@'$'$'$R***BrcbeZdZdZejffd Zd dZdZdZ dZ dZ d Z 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 cp||_||_||f}t||dSr)r.r rr)rr.r rrs rrzSakurai.__init__ds: A &&&&&rNc ||j}tj|jdz|z |jdz}tjdtjtjd|ztjz|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. Nrg0@r7)r.rr9flippowerr;r&)rr*rMs rr1zSakurai.eigenvaluesjsq 9A Idfqj!mTVaZ 0 0wsRXbfS1Wru_ -K&L&LaPPPQQQrcXtjddtj|jdz |jzdf}dtj|j|jz}tj|j|j}tj|||g|jS)zA Construct the Sakurai matrix as a banded array. rr\r)rr_r<r.r arrayrg)rd0d1d2s rtobandedzSakurai.tobanded}sU1a"'$&1*DJ????B C "'$& 333 3 WTV4: . . .xR %%,,TZ888rcddlm}|}||d|d|d|d|dggd|j|jf|jS)zB Construct the Sakurai matrix in a sparse format. r diags_arrayrr)r]r8rrroffsetsrr ) scipy.sparserrr.r )rrds rryzSakurai.tosparses| -,,,,, MMOO{AaD!A$!adAaD9CTCTCT"&&$&!1BBB BrcN|Srryrors rrozSakurai.toarray}}&&(((rcp||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. r8r\rrNrrrr]r)rrr^rq))rrr^)rfr.r promote_typesr zeros_likepad)rrT result_dtypesxs rrzSakurai._matvecsz IIdfb ! !'<< ]1L 1 1 1qAAAw;Qq!!!tW,qAAAw61aaa4"aaa%L1qQQQx</!BE(:2qqq5 AaeQQQhK!q"aaay1QRRU8/C*DDq"aaay*:;;<qQQQx)9::;1b5!!!8  rc,||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``. rrs rrzSakurai._matmat ||Arc|SrrGrs rrzSakurai._adjointrrc|SrrGrs rrzSakurai._transposerrr)rrrrrr`rr1rryrorrrrrrs@rrr sYYt!#'''''' RRRR&999 B B B))) rrc`eZdZdZejffd 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. ch||_||_t||dSr)rr rr)rrr rs rrzMikotaM.__init__s1   &&&&&rc~dtjd|jddzz |jS)Nr6rr)rr9rrgr rs r_diagz MikotaM._diags6RYq$*Q-!"3444<A3!%4:??? ?rc~tj||jSr)rdiagrrgr rs rrozMikotaM.toarrays*wtzz||$$++DJ777rc||jdd}|ddtjf|zS)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. rr8N)rfrrrnewaxisrs rrzMikotaM._matvecs? IIdjmR ( (zz||AAArzM*Q..rc,||S)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``. rrs rrzMikotaM._matmatrrc|SrrGrs rrzMikotaM._adjointrrc|SrrGrs rrzMikotaM._transposerr)rrrrrr@rrrryrorrrrrrs@rrrs.%'J'''''' III ??? 888///rrcZeZdZdZejffd ZdZdZdZ dZ dZ dZ d Z xZS) 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. c||_||_t|||d}t jd|zdz dd|j|_t j|dz dd|j |_dS)Nrrrr]r\r8)rr rrrr9_diag0_diag1)rrr r.rs rrzMikotaK.__init__s   &&& !HiA 1b CCC  !a%BdjAAAA rcjtjtj|jdd|jgS)Nrrconstant)rrrrrrs rrzMikotaK.tobandeds)x VZ@@$+NOOOrchddlm}||j|j|jggd|j|jS)Nrrrsr)rrrrrr rs rryzMikotaK.tosparsesK,,,,,,{DKdkBJJJ!%4:??? ?rcN|Srrrs rrozMikotaK.toarray#rrcN||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. rr8r\Nrr]r)rfrrrr rrr)rrTrkxrrs rrzMikotaK._matvec&sb IIdjmR ( ('<< ]1L 1 1 1 [ [a51QT7?RUQq!!!tW_41aaa4rFQr111uX%22qqq5(992qqq5 3B39 $B$' 2QUD[/AaeQQQhK78QRRX,122qqq5121b5!!!8  rc,||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``. rrs rrzMikotaK._matmat8rrc|SrrGrs rrzMikotaK._adjoint@rrc|SrrGrs rrzMikotaK._transposeCrr)rrrrrint32rrryrorrrrrrs@rrrs0%'HBBBBBBPPP??? )))$rrc.eZdZdZejfdZddZdS) 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|_dSr)r.r rrr*rrM)rr.r s rrzMikotaPair.__init__sI V TZ00TZ00rNcf||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. Nrr\)r.rr9uint64)rr* arange_plus1s rr1zMikotaPair.eigenvaluess7 9AyAE;;; l**rr)rrrrrr@rr1rGrrrrGsPWWp!# 1111++++++rr)numpyrscipy.sparse.linalgrrrrr__all__rrrrrrGrrrsP......3333333333 / |||||.|||~gggggngggTBBBBBnBBBJLLLLLnLLL^q+q+q+q+q+q+q+q+q+q+r