EVP returns correct eigenvalues but clearly wrong eigenvectors

[Christopher Marcotte]

I'm trying to reproduce a simple chebfun eigenvalue example (http://www.chebfun.org/examples/ode-eig/Drum.html) using Dedalus' EVP methods. Fortunately, the eigenvalues are correct. However, the eigenvectors are obviously wrong and exhibit the Nyquist mode when plotted. Have I made a fundamental mistake with the output? 

[Daniel Lecoanet]

Hi,

You only sorted your eigenvalues, not your eigenvectors.  To do this, you could do something like

order = np.argsort(solver.eigenvalues)

np.print(solver.eigenvalues[order[:5]])

for i in range(5):

  plt.plot(x,solver.eigenvectors[:,order[i]])

Daniel

On Tue, Feb 14, 2017 at 5:46 AM, Christopher Marcotte <christophe...@gmail.com> wrote:

I'm trying to reproduce a simple chebfun eigenvalue example (http://www.chebfun.org/examples/ode-eig/Drum.html) using Dedalus' EVP methods. Fortunately, the eigenvalues are correct. However, the eigenvectors are obviously wrong and exhibit the Nyquist mode when plotted. Have I made a fundamental mistake with the output? 

--
You received this message because you are subscribed to the Google Groups "Dedalus Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dedalus-users+unsubscribe@googlegroups.com.
To post to this group, send email to dedalu...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/CAESDi8TXk3PWSpU2z1_RiUOD6X1GXVecvy5CwdxYuyUNwHG4zw%40mail.gmail.com.
For more options, visit https://groups.google.com/d/optout.

[Jeffrey S. Oishi]

Hi Chris,

Dedalus returns eigenvalues and eigenvectors in the order they were received from the numerical eigenvalue solver. This can lead to some unexpected behavior, as you've found. I've attached a corrected script which fixes the following issues:

* when you followed the waves on a string script and sorted the eigenvalues, that had no effect on the order of the eigenvectors. Instead of sorting the eigenvalues, I used argsort, which returns a numpy array of indices in order. Then, you can use that array of indices to sort both the eigenvalues and vectors. Now, order here is defined by the default sort algorithm of argsort; if you wanted to sort in some other way, that too is possible, but requires a bit more doing.

* We have a nice routine on the solver, solver.set_state(mode) which sets the data in the state variable to the mode numbered by "mode". I've used that here in the loop over the first 6 modes.

* In order to use set_state(), you have to make the domain complex. This is a known limitation, since the solutions on the grid for the drum problem are real valued.

* Just for fun, I used scipy.special's jn_zeros(0,6) to compute the first six zeros of J0 so you don't have to copy them in as text in the script.

With these modifications, the eigenvectors now looks pretty reasonable to me.

As a follow on, we have a package outside of Dedalus called eigentools, https://bitbucket.org/jsoishi/eigentools, written by Susan Clark and myself, which does things like select only the "good" eigenvalues by solving the EVP at resolution n and 3*n/2 (or any other fraction you choose) and only selecting modes that are the same between the two runs (for some value of "the same"; for details see Boyd Ch 7). It also does nice things like plot spectra and find the fastest growing mode. eigentools is designed to be most helpful when you have an eigenvalue problem with multiple (possibly overlapping) wave families, in which case determining the "order" for good eigenvalues can be tricky. It's still very much a work in progress, but you might find it useful.

Please let us know if you have any more questions!

Jeff

On Tue, Feb 14, 2017 at 5:46 AM Christopher Marcotte <christophe...@gmail.com> wrote:

I'm trying to reproduce a simple chebfun eigenvalue example (http://www.chebfun.org/examples/ode-eig/Drum.html) using Dedalus' EVP methods. Fortunately, the eigenvalues are correct. However, the eigenvectors are obviously wrong and exhibit the Nyquist mode when plotted. Have I made a fundamental mistake with the output? 

--

You received this message because you are subscribed to the Google Groups "Dedalus Users" group.

To unsubscribe from this group and stop receiving emails from it, send an email to dedalus-user...@googlegroups.com.

[Keaton Burns]

Hi Chris,

To clarify a bit here, besides the sorting, the real issue was that your script was originally plotting the coefficients, not the grid-space values, of the eigenvectors.  This is because Dedalus uses a tau method to construct and solve problems as sparse and (typically) banded matrices acting on the Chebyshev coefficients, as opposed to dense collocation matrices like Chebfun.  The eigenvectors returned by the eigenvalue solver are therefore the interleaved coefficients of the problem variables (i.e. the coefficient state vector).  What the “set_state” method does is properly extract each variables coefficient from those vectors, so you can then use the Dedalus transforms to view the eigenmodes in grid space, as Jeff showed.

Best,

-Keaton

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/CAAOZTN_%3DOaYMB3CTyRjEu3tQTU64pXp3TLpRyiXT5%3D%3DH77kcvA%40mail.gmail.com.

[Christopher Marcotte]

Okay, so not only did I not understand the output but I obfuscated it further by rearranging the eigenvalues and leaving the eigenvectors as is.

Jeff, I will have a look at eigentools.

Thank you Daniel, Jeff, and Keaton!

[Christopher Marcotte]

Perhaps there's a very good programmatic reason for this but why is the output of EVP not an array of fields, as in the BVP and IVP solutions?

CM

[Jeffrey S. Oishi]

Hi Chris,

It is, with the caveat that here, your data is actually n eigenmodes, each of length n*n_fields. The purpose of set_state is to properly copy, in coefficient space, the data for a given mode to the usual state vector. Thus, having done set_state, you have

solver.state['u']['g']

solver.state['u_x']['g']

or

solver.state['u']['c']

solver.state['u_x']['c']

for grid and coefficient space, respectively. We provide the state vector as the way to interface with data from IVP, BVP, and EVPs. EVPs are just a bit different because of the way the data is returned.

j

--

You received this message because you are subscribed to the Google Groups "Dedalus Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dedalus-user...@googlegroups.com.
To post to this group, send email to dedalu...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/7c84d8a2-a7db-4d59-8c80-5d01b7c62619%40googlegroups.com.

[Christopher Marcotte]

I see, thanks Jeff!

CM

[Daniel Lecoanet]

Hi,

The real reason is that I was the one who originally coded the eigenvalue solver, and I never do anything for good programatic reasons.

Daniel

On Wed, Feb 15, 2017 at 1:11 PM, Christopher Marcotte <christophe...@gmail.com> wrote:

I see, thanks Jeff!

CM

On Wednesday, February 15, 2017 at 6:00:00 PM UTC, J. S. Oishi wrote:

To unsubscribe from this group and stop receiving emails from it, send an email to dedalus-users+unsubscribe@googlegroups.com.

To post to this group, send email to dedalu...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/83fe826f-cc00-4b78-9e5b-e0ff5e4b7559%40googlegroups.com.