31

Very good observation. Indeed, this issue is really frustrating. To single out the issue: It seems that Arnoldi's method is to blame: Max@Abs@Im@Eigenvalues[M, -Nless] Max@Abs@Im@Eigenvalues[M, -Nless, Method -> "Arnoldi"] Max@Abs@Im@Eigenvalues[M, -Nless, Method -> "Direct"] Max@Abs@Im@Eigenvalues[M, -Nless, Method -> "FEAST"] 0.000610613 ...


30

Not a solution but too big for a comment. There seems to be a catastrophic failure in Eigenvalues happening that is not due to the matrix being crazy. As a diagnostic, let's calculate the smallest (by absolute value) eigenvalue of the upper-left $n\times n$ part of the matrix M = mat[xlist[[3]]]; For odd $n$ the answer is zero, so let's only do this for ...


27

You could use MatrixRank. Here is a function that does this: eigenvectorQ[matrix_, vector_] := MatrixRank[{matrix . vector, vector}] == 1 For your example: eigenvectorQ[h, y] False


21

This is a bug. The problem seems to be related to the fact that some rows of your matrix are packed arrays and some are not. Developer`PackedArrayQ /@ M (* {True, False, False, True, False, False, False, False, False, False, False, False, True, False, False, True} *) Unpacking seems to avoid the problem. Eigenvalues@Developer`FromPackedArray[M] (* {-1, -...


19

Use Eigenvalues[mat, Cubics -> True] Eigenvectors[mat, Cubics -> True] sometimes Quartics -> Truecan be needed. or ToRadicals @ Eigenvalues[ mat] ToRadicals @ Eigenvectors[ mat] In general one cannot find roots (of higher order) polynomials in terms of radicals. The reason that Mathematica allows this option is that in general it is ...


18

Update: This implementation is now a package called CompoundMatrixMethod, hosted on github. It can be installed easily by evaluating: Needs["PacletManager`"] PacletInstall["CompoundMatrixMethod", "Site" -> "http://raw.githubusercontent.com/paclets/Repository/master"] This version also includes a function ToMatrixSystem which converts a system of ODEs ...


18

The following is an attempt to recreate a similar sort of interactive visualization, showing the eigenvectors (when real), and how the various points of the unit circle are transformed by the matrix. The matrix can be chosen by moving its two column vectors using the mouse. I used EventHandler for this, instead of Locators, for greater customizability and a ...


12

Eigenvalues and Eigenvectors return the eigenvalues sorted from highest to lowest. From the documentation: Eigenvectors with numeric eigenvalues are sorted in order of decreasing absolute value of their eigenvalues. This means that the last three eigenvectors returned by Eigenvectors will correspond to the lowest-energy states: ListPlot[{evecs[100, 98],...


11

I have a package for solving eigenvalue boundary value problems using the Compound Matrix Method with the Evans function, which is ideally suited to this. The package is available on my GitHub, a partial explanation is at my previous question, and a good introduction to the method is in this pdf. You can install the package using: PacletInstall["...


11

Edit: Eigensystem fixed in 12.1 in addition to Eigenvalues I attempted a workaround, to see if Eigensystem had any issues also. It does. This is very unfortunate. (Will we have to wait for 12.1 for the fix (?!)) (We waited for 12.1 for the fix (!!)) My code here: e3[n_?EvenQ] := Eigensystem[M[[;; n, ;; n]]][[1]] // Abs // Min Produces the following, ...


10

Chebyshev Series solution A change of variables $\left\{ x \mapsto \frac{L}{2} (t+1)\,,\ w(x) \mapsto \frac{L}{2} u(t)\right\}$ converts the OP's differential equation to $$(t-1) u''(t)+u'(t)=\left(\frac{2}{L}\right)^3 u^{(4)}(t),\quad -1\le t\le 1$$ Integrating four times, we get (with $\int^k$ denoting repeated integration with respect to $t$) the ...


10

NDEigenValues handles the pair of first-order equations in the question much more accurately, when it is converted into a single second-order equation. eqe = Eliminate[Thread[Flatten[{D[(lhs - λ variables) // First, r], lhs - λ variables}] == 0], {α[r], α'[r]}] // First (* β''[r] == -36 β[r] + β[r]/r^2 + 100 r^2 β[r] - λ^2 β[r] - β'[r]/r *) where λ ...


10

The additional problem added to the end of the question can be solved in a similar manner. Begin with variables = {α[r], β[r], γ[r], δ[r]}; lhs = {Fop1[δ[r], l + 1, -1] + b γ[r] + m0[R, r] α[r], Fop1[γ[r], l, 1] - b δ[r] + m0[R, r] β[r], Fop1[β[r], l + 1, -1] + b α[r] - m0[R, r] γ[r], Fop1[α[r], l, 1] - b β[r] - m0[R, r] δ[r]} and ...


10

We want to find the value of $\lambda$ for which there exists a solution of the differential equation $u'' = \lambda ( - u + u^2)$, subject to the boundary conditions $u(0) = u(1) = 0$. We can do this by a shooting algorithm, i.e., solving the related initial value problem $$ u'' = \lambda (- u + u^2), \qquad u(0) = 0, \; u'(0) = v. $$ The value of $u(1) \...


10

Here is a modified version of my answer to How to numerically solve a 1-d time-independent Schrödinger equation? that solves this problem. I use NDEigensystem after manually shifting the potential so its bottom is at zero energy. This way the eigenvalues are automatically sorted in ascending order. In the output, I undo that shift so the eigenvalues are ...


10

There are three conditions when we want to get eigenfunctions in Cartesian coordinates, similar to eigenfunctions in cylindrical coordinates. The first is the correspondence of boundaries. The second is the azimuth number match, eg $l_1=l_2=0$.Third, the radius of the outer circle must meet the boundary condition. All three conditions are violated in the ...


10

I think you hit an Indeterminate form. H = {{0, 0, Exp[I k1] + m1, Exp[I k2] + m2}, {0, 0, Exp[I Phi] (Exp[-I k2] + m2), Exp[-I k1] + m1}, {Exp[-I k1] + m1, Exp[-I Phi] (m2 + Exp[I k2]), 0, 0}, {Exp[-I k2] + m2, Exp[I k1] + m1, 0, 0}}; {eigval, eigvec} = Eigensystem[H]; vvv = {eigvec[[1]], eigvec[[3]]}; Now see what happens here (vvv /. ...


9

Here is how to get the result directly from @Vsevolod analytic solution: seeking a solution to the case w[0]==w'[0]==w''[L]==w'''[L]==0, substituting w'[x]=v[x] : vf[x_, L_, c1_, c2_] = DSolveValue[{v'''[x] + (L - x) v'[x] - v[x] == 0, v[0] == 0}, v[x], x] /. {C[1] -> c1, C[2] -> c2}; now for w'''[L] == v''[L] ==0: Limit[(D[vf[x, L, c1, c2]...


9

Weak method, the solution converges to the exact one at Nn >= 200. It takes time. Figure 1 shows the solution for Nn=100, 200 in comparison with the exact solution. Nn = 200; h = \[Pi]/Sqrt[2 Nn]; h1 = h; exactdelta2 = h1^2 D[S1[\[Phi]1[x], j, h1], {\[Phi]1[x], 2}]; ex5 = (exactdelta2 /. x -> \[Psi]1[x]) /. \[Phi]1[\[Psi]1[x_]] :> x /. x -&...


9

The reason is that your second matrix is a rank-one update of your first matrix: $$ B\equiv uv^t $$ where $u=(1,1,1)$ and $v=(d,e,f)$. Therefore, the new eigenvalues are typically a small perturbation of the old ones, and there are some known formulas for special cases. See e.g. these lectures or the references in this math.OF post.


8

In versions 10.3 and higher, there's DEigensystem and NDEigensystem: DEigensystem[{-Laplacian[y[x], {x}], DirichletCondition[y[x] == 0, True]}, y[x], {x, 0, 1}, 4] (* {{π^2, 4 π^2, 9 π^2, 16 π^2}, <-- eigenvalues {Sin[π x], Sin[2 π x], Sin[3 π x], Sin[4 π x]}} <-- eigenfunctions *) They require a linear operator, but it ...


8

I'm not sure this will be robust enough for a 48 x 48 matrix. The idea is common enough: Build an interpolation of a function by integrating its derivative with NDSolve, using a particular value of the function as an initial condition. If the derivative is computed as continuous, then the functions we get should track an eigenvalue trajectory through a ...


8

Serious Update I decided this was useful enough shit to package up. The package is here. One of these days I'll integrate this into my main chemistry package. In any case, for our purposes here we can load it like this: Get["https://github.com/b3m2a1/mathematica-tools/raw/master/BasisSetSchrodinger.wl"] Then we can call WavefunctionEigensystem to get the ...


8

The ODE itself can be solved symbolically. For generality, assume that k and b are undefined and that ll is the unknown eigenvalue. Then, s = DSolveValue[lhs == ll f[y], f[y], y] // Simplify (* E^(-(((-2 + k) (-2 + 3 b y))/(9 b))) (-6 + 9 b y)^(Sqrt[16 + 32 k + 16 k^2 - 9 ll]/(9 b)) (C[1] HypergeometricU[(-4 + 9 b - 4 k + Sqrt[16 + 32 k + 16 k^2 - 9 ...


8

By default, the eigenvalues are ordered by absolute value. All the eigenvalues of this particular matrix have the same absolute value plus some rounding errors. Thus, it can easily happen, that the fourth eigenvalue is positive or negative, depending on the parameters. You can use Max to plot the largest eigenvalue: Plot3D[Max@Eigenvalues[H[0.1, 0.5, 0.7, ...


8

Fixed in 12.1 ClearAll[x, n]; NN = 374; R = 0.05; t1 = -1 + Cos[x] - I Sin[x] + I R; t1p = -1 + Cos[x] + I Sin[x] + I R; mat[x_] = DiagonalMatrix[Table[If[EvenQ[n], t1, -1], {n, 0, 2 NN - 1 - 1}], 1] + DiagonalMatrix[ Table[If[EvenQ[n], t1p, -1], {n, 0, 2 NN - 1 - 1}], -1] + DiagonalMatrix[Table[If[EvenQ[n], -1, 0], {n, 0, 2 NN - 1 - 3}], ...


8

Use the Arnoldi method with shift-inversion: Eigenvalues[A, 3, Method -> {"Arnoldi", "Criteria" -> "Magnitude", "Shift" -> 0}] gives you the three smallest eigenvalues by absolute value (by magnitude). See here: Efficiently find all values of parameter such that any of the eigenvalues of a matrix is equal to 1 After comments by @HenrikSchumacher ...


8

This question is particularly interesting to me, and I have a package that may be helpful to you here. This particular equation is: Fourth order Linear Inhomogeneous in the independent variable Contains Robin-like boundary conditions Includes the eigenvalue in the boundary conditions I don't think that NDEigensystem can handle fourth-order systems, or ...


8

Classical problem. You want to compute the eigensystem of the second fundamental form with respect to the first fundamental form. Thus you have to solve a generalized eigensystem. This can be done with Eigensystem[{FF2, FF1}], but it does not work very well with symbolic functions. Also, one has to normalize the eigenvectors for some reason I don't get. But ...


7

Here's a range of incremental improvements: Table[Normalize[v[[i]]], {i, 2}] Table[Normalize[v[[i]]], {i, Length[v]}] (* like a for loop in other languages *) Table[Normalize[ev], {ev, v}] (* like a for each loop in other languages *) Map[Normalize, v] Normalize /@ v (* shorthand for Map *) Look up Map and all the syntaxes of Table and read this: ...


Only top voted, non community-wiki answers of a minimum length are eligible