2
$\begingroup$

Consider this code

import numpy as np
from numpy import pi
from matplotlib import pyplot as plt
import scipy.linalg as linalg

## Domain
L = 5
M = 50
dx = L/(M-1)
xx = np.linspace(0, L, M)


## Diff matrix
row0 = np.zeros(M)
row0[[-1, 0, 1]] = [1, -2, 1]
row0 /= dx**2
D = linalg.circulant(row0)

# Dirichlet BCs -- implicitly set boundaries to 0
D = D[1:-1, 1:-1]

# Use negative laplacian
D *= -1


## Eigenvalue decomposition
ews, evs = linalg.eigh(D)

# ews, evs = linalg.eig(D)
# ews = ews.real


## Plot
fig, (ax1, ax2) = plt.subplots(nrows=2)
modulation = np.cos(pi/dx * xx[1:-1])  # [-1, +1, -1, +1, ...]
for ev in evs.T[:3]:
    ax1.plot(xx[1:-1], ev)
    ax2.plot(xx[1:-1], ev*modulation)
ax1.set_title("Unmodulated")
ax2.set_title("Modulated by [-1, +1, -1, +1, ...]")
plt.show(block=False)
fig.tight_layout()

with the following output

enter image description here

Now, if I don't do D *= -1 then the contents of the two panels are swapped (i.e. the modulation becomes necessary to recover the true eigenvectors). But theoretically, the eigenvectors should not have changed. So why do I get different answers?

Note1: both the modulated and unmodulated eigenvectors pass the validation of comparing ev to D @ ev. But the eigenvectors are supposed to be unique (up to a scaling), so what gives?

Note2: The issue seems related to aliasing, since the modulation is by a high-frequency signal. But I can't fathom what that has to do with eigenvalue problems.

Note3: If I uncomment the use of the (general) eig function, then it does not matter whether D *= -1 was done or not (the eigenvectors always require the modulation). Isn't the wide discrepancy of the output disconcerting for the "stability" of these routines? I mean, sure, as iterative methods, they depend on the initial guesses, but the problem seems well-posed and simple to me, but the output varies greatly.

Note4: this is not the same issue as in here, as I take care to transpose the eigenvector matrix.

Note5: by modulation I mean that the original signal/vector/function gets (pointwise) multiplied by the other, as seen in the code (ev*modulation)

$\endgroup$
7
  • 1
    $\begingroup$ Please explain what you call "modulation" by [-1,+1,-1...] in mathematical terms : which action on your initial matrix ? $\endgroup$
    – Jean Marie
    Commented Mar 23, 2021 at 17:24
  • $\begingroup$ @JeanMarie I've added Note5 to explain it. $\endgroup$
    – Patrick
    Commented Mar 23, 2021 at 17:59
  • $\begingroup$ Can you please post a self-contained script that reproduces exactly the issue you describe? This will make it easier to tell whether it is simply a coding mistake. $\endgroup$
    – Simon Segert
    Commented Mar 23, 2021 at 18:08
  • 3
    $\begingroup$ Another possible issue is that the eigenvectors are listed in the opposite order when you multiply the matrix by -1. This would happen if the eigh function by default returns the eigenvalues in order from largest to smallest $\endgroup$
    – Simon Segert
    Commented Mar 23, 2021 at 18:13
  • 2
    $\begingroup$ In fact Mike Hawk has located the issue: you are displaying for some reason the final eigenvectors instead of the first ones... $\endgroup$
    – Jean Marie
    Commented Mar 23, 2021 at 18:33

1 Answer 1

Reset to default
0
$\begingroup$

As pointed out in the comments, I had simply forgotten that I needed to sort the eigenvalues and vectors, whose order otherwise depends on the sign of the matrix operator and the choice of numerical routine.

$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .