-iSSKrSSKJr SSKJrJrJr S/r"SS\5r "SS\5r "SS \5r "S S \5r "S S 5r g)N)LinearOperator)kron eye_array dia_array LaplacianNdc^\rSrSrSrS\R S.U4SjjrSrSSjr Sr S r SS jr S r S rS rSrSrSrSrU=r$)r aI 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.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>US;a[SU<S35eXlX l[R"U5n[ TU]X4U4S9 g)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.13/site-packages/scipy/sparse/linalg/_special_sparse_arrays.pyrLaplacianNd.__init__s_ &J J !4 78DD  %#6 GGJ  uF3c URnUc-[R"U5n[R"U5nOX[ U[ [R "U5U-55n[R"U5n[R"U5n[X25HupgURS:Xa>US[R"[RUS--SUS--- 5S--- nMSURS:Xa8US[R"[RU-SU-- 5S--- nMUS[R"[R[R"US-S- 5-U- 5S--- nM UR5n[R"U5n Xn Ub X*Sn X*Sn X4$)zCompute `m` largest eigenvalues in each of the ``N`` directions, i.e., up to ``m * N`` total, order them and return `m` largest. Nrr )rrindiceszerosmintuple ones_likezipr sinpifloorravelargsort) rmrr"Leiggrid_shape_minjn Leig_ravelind eigenvaluess r_eigenvalue_ordering LaplacianNd._eigenvalue_orderings__ 9jj,G88J'D !&r||J'?!'C!DFNjj0G88N+D,DA'';6RVVBEEQUOqAE{$CDIII))Y6RVVBEEAIQ$78A===RVVBEEBHHa!eq[,A$AA$EF!KKK -ZZ\ jj$ o =%bc*Kbc(Crc,URU5up#U$)a>Return 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. )r5)rr-r4_s rr4LaplacianNd.eigenvalues%s2215 rctURS:Xaj[R[R"U5S--US-- n[R"SUS-- 5[R "X1S--5-nGOuURS:Xah[R[R"U5S--U- n[R"US:XaSOSU- 5[R "X1-5-nOUS:Xa1[R"SU- 5[R"U5-nOUS-U:Xa@US-S:Xa7[R"SU- 5[R"SS /US-5-nO}S[R-[R"U5S--U- n[R"SU- 5[R "U[R"US-S- 5-5-nS U[R"U5[R"[R5R:'U$) zYReturn 1 eigenvector in 1d with index `j` and number of grid points `n` where ``j < n``. rr g@?r ?rr!g)r rr)arangesqrtr(cosonestiler*absfinfofloat64eps)rr0r1ievs r_ev1dLaplacianNd._ev1d6s  # #{ 21)*a!e4Aq2v'"&&!e*==B  % % 21+,q0AQ"B!34rvvae}DBAvWWR!V_rwwqz1Q!A WWR!V_rww2w1'==J"))A,"459WWR!V_rvva"((AEQ;2G.G'HH57266":,00 01 rc[XR5VVs/sHup#URX#5PM nnnUSnUSSHn[R"XTSS9nM [R "U5R 5$s snnf)zjReturn 1 eigenvector in Nd with multi-index `j` as a tensor product of the corresponding 1d eigenvectors. rr N)axes)r'rrIr tensordotasarrayr+)rkr0r1phiresults r_one_eveLaplacianNd._one_eveMsu-0??,CD,CDAtzz!,CDQqr7C\\&A6Fzz&!'')) EsBcURU5up#Uc URnO@[UR[[R "UR5U-55n[R "X45n[U6Vs/sHn[U5PM nnUVs/sHopRU5PM nn[R"U5$s snfs snf)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. ) r5rr$r%rr& unravel_indexr'rR column_stack) rr-r8r3r/ N_indicesxrOeigenvectors_lists r eigenvectorsLaplacianNd.eigenvectorsWs**1- 9!__N  %bll4??&Ca&G HJN$$S9 '*I7!U1X 77@Ay!]]1-yA0118As CCc URn[R"U5n[R"X"/[RS9n[R "U5n[R "U5n[ U5GHupgSUSS&S[R"SUSU2SU245SS&S[R"SUSUS- 2SU245SS&S[R"SUSU2SUS- 245SS&URS:XaSUS 'SXGS- US- 4'OGURS :Xa7US:a$USUS- 4==S- ss'XGS- S4==S- ss'O US ==S- ss'UnUS:aT[R"USU5n [SU 5H+n USU2SU24XJU-U S-U-2X-U S-U-24'X- nM- USU2SU24USU2SU24'[[R"XS-S55n SUSU2SU24'[U 5V s/sHoPM n n USU2SU24URXX45SS2U SS2U 4'X4- nGM URUR5$s sn f) 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->ir r r=rrr)rrrr#int8 empty_like enumerateeinsumr rangeintreshapeastyper ) rrr1LL_iLtempr3dimnew_dimtilesr0rXidxs rtoarrayLaplacianNd.toarrayrs__ GGJ  HHaV277 +mmA a !*-HCCF68BIIgs4C4#: / 2;(?E(7(HWH$ % q56 234E&'C(7(" ##El+l1lC+%*(7(HWH*<$= KK! S!S. " HAW.Zxx ##,s I%c [UR5n[R"UR5n[ X"4[R S9n[ U5GHanURUn[R"SU/[R S9nUSSS24==S-ss'URS:Xa SUS'SUS '[ U/S Q4XU4[R S 9nURS :XaF[ XU4[R S9nURS/U*S-S 9 URS/US- S 9 Xx- n[ U5H2n [[URU [R S9U5nM4 [ US-U5H2n [U[URU [R S95nM4 X7- nGMd URUR5$)z 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]r Nr^r r=r r)r r=r=rr )rr r)rO)lenrrrrr`rdrAr setdiagrrrgr ) rrprhrGrkdataritr0s rtosparseLaplacianNd.tosparses   GGDOO $ qfBGG ,qA//!$C77As82773D AJ" J''94T  U T:.sj"$''C'':5sj8 1##a ( 1#Q '1X9T__Q%7rwwGM1q5!_3 $//!*".*26>$@'$R*Brcz^\rSrSrSr\R 4U4SjjrS SjrSr Sr Sr Sr S r S rS rS rU=r$)SakuraiiaT 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 cB>XlX lX4n[TU] X#5 gr)r1r rr)rr1r rrs rrSakurai.__init__as!  &rc TUc URn[R"URS-U- URS-5n[R"S[R"[R "SU-[R -URS-- 5S5-5$)aEReturn 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. r g0@r<)r1rr>flippowerr@r))rr-rOs rr4Sakurai.eigenvaluesgsw 9A IIdffqj!mTVVaZ 0wwsRXXbffS1Wruu_ -K&LaPPQQrc[RSS[R"URS- URS9-S4nS[R"URURS9-n[R"URURS9n[R "X2U/5R UR5$)z1 Construct the Sakurai matrix as a banded array. r!r]r)rr_rAr1r arrayrg)rd0d1d2s rtobandedSakurai.tobandedzsUU1a"''$&&1*DJJ??B C "''$&& 3 3 WWTVV4:: .xx %,,TZZ88rcSSKJn UR5nU"USUSUSUSUS//SQURUR4URS9$)z2 Construct the Sakurai matrix in a sparse format. r diags_arrayr r!)r^r=rr r!offsetsrr ) scipy.sparserrr1r )rrds rrzSakurai.tosparses] - MMOAaD!A$!adAaD9CT"&&&$&&!1B Brc>UR5R5$rrzrors rroSakurai.toarray}}&&((rcLURURS5n[R"URUR5n[R "XS9nSUSSS24-SUSSS24-- USSS24-USSS24'SUSSS24-SUS SS24-- US SS24-USSS24'S USS2SS24-SUSS 2SS24USS2SS24--- [R "USS 2SS24S 5-[R "US S2SS24S5-USS2SS24'U$)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. r=r]rrNrr r!r^r)rsr_rr))rr r_)rfr1r promote_typesr zeros_likepad)rrX result_dtypesxs rrSakurai._matvecs9 IIdffb !''< ]]1 1qAw;Qq!tW,qAw61a4"a%L1qQx</!BE(:2q5 AaeQhK!q"ay1QRU8/C*DDq"ay*:;<qQx)9:;1b5!8  rc$URU5$)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 rrSakurai._matmat ||ArcU$rrrs rrSakurai._adjointrrcU$rrrs rrSakurai._transposerr)r r1r)rrrrrrr`rr4rrzrorrrrrrrs@rrrsGYt!#' R&9 B) rrcv^\rSrSrSr\R 4U4SjjrSrSr Sr Sr Sr S r S rS rS rU=r$) MikotaMia' 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. c<>XlX l[TU] X!5 gr)rr rr)rrr rs rrMikotaM.__init__s   &rcS[R"SURSS-5- RUR5$)Nr;r r)rr>rrgr rs r_diag MikotaM._diags6RYYq$**Q-!"344<A3!%4::? ?rc|[R"UR55RUR5$r)rdiagrrgr rs rroMikotaM.toarrays&wwtzz|$++DJJ77rcURURSS5nUR5SS2[R4U-$)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. rr=N)rfrrrnewaxisrs rrMikotaM._matvecs: IIdjjmR (zz|ArzzM*Q..rc$URU5$)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 rrMikotaM._matmatrrcU$rrrs rrMikotaM._adjointrrcU$rrrs rrMikotaM._transposerrr)rrrrrrrErrrrzrorrrrrrrs@rrrsD.%'JJ' I ? 8/rrcp^\rSrSrSr\R 4U4SjjrSrSr Sr Sr Sr S r S rS rU=r$) MikotaKiaQ 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>XlX l[TU] X!5 USn[R "SU-S- SSURS9Ul[R "US- SSURS9*Ulg)Nrr!r r^r]r=)rr rrrr>_diag0_diag1)rrr r1rs rrMikotaK.__init__sh   & !HiiA 1b C  !a%BdjjAA rc[R"[R"URSS5UR/5$)Nrsconstant)rrrrrrs rrMikotaK.tobandeds+xx VZ@$++NOOrcSSKJn U"URURUR//SQURUR S9$)Nrrrtr)rrrrrr rs rrzMikotaK.tosparses6,DKKdkkBJ!%4::? ?rc>UR5R5$rrrs rroMikotaK.toarray rrcURURSS5n[R"URUR5n[R "XS9nUR nURnUSUSSS24-USUSSS24--USSS24'USUSSS24-USUSSS24--USSS24'USS2S4USS2SS24-USS2S4USS2SS24--USS2S4USS2SS24--USS2SS24'U$)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. rr=r]Nr r^r!)rfrrrr rrr)rrXrkxrrs rrMikotaK._matvec#s7 IIdjjmR (''< ]]1 1 [[ [[a51QT7?RUQq!tW_41a4rFQr1uX%22q5(992q5 3B39 $B$' 2QUD[/AaeQhK78QRX,12q5121b5!8  rc$URU5$)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 rrMikotaK._matmat5rrcU$rrrs rrMikotaK._adjoint=rrcU$rrrs rrMikotaK._transpose@rr)rrr r)rrrrrrint32rrrzrorrrrrrrs@rrrs@0%'HHBP? )$rrcB\rSrSrSr\R 4SjrSSjrSr g) MikotaPairiDat 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]) cXlX lX4Ul[URUR5Ul[ URUR5Ulgr)r1r rrr-rrO)rr1r s rrMikotaPair.__init__sA V TZZ0TZZ0rNcvUc URn[R"SUS-[RS9nX"-$)aDReturn 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. r r])r1rr>uint64)rr- arange_plus1s rr4MikotaPair.eigenvaluess5 9AyyAE; **r)r rOr-r1rr) rrrrrrrErr4rrrrrrDsWp!# 1+rr)numpyrscipy.sparse.linalgrrrrr__all__rrrrrrrrrs`.33 / y.yxgngTBnBJLnL^q+q+r