Eigenvalue problem: locating complex eigenvalues via frequency scan

[member 428835]

Hi PF!

Here's an ODE (for now let's not worry about the solutions, as A LOT of preceding work went into reducing the PDEs and BCs to this BVP):

$$\lambda^2\phi-0.1 i\lambda\phi''-\phi'''=0$$ which admits analytic eigenvalues

$$\lambda =-2.47433 + 0.17337 I, 2.47433 + 0.17337 I, -10.5087 + 1.16033 I, 10.5087 + 1.16033 I, -22.0562 + 3.13425 I, 22.0562 + 3.13425 I...$$

There's an analytic technique to solve which is not terribly difficult (again, you'd need the entire system of equations to do this). However, another way to solve this ODE is to change the nature of the problem from an eigenvalue problem to a frequency response problem, where now the eigenvalues $\lambda$ are instead input parameters of the forcing frequency:

$$ \lambda^2\phi-0.1 i\lambda\phi''-\phi'''=f \lambda$$

By carefully selecting basis functions $\phi_j = \sin(\pi(j-1/2)x)\cosh( \sqrt{ (\pi(j-1/2))^2 + 1} )$ we can take inner products of this equation, resulting in a square matrix linear system of algebraic equations:

$$\lambda^2 M-0.1 i \lambda V - K = F : M_{ij} = \int_0^1 \int_0^1 \phi_i(x) \phi_j(y)\,dy\,dx,\\
M_{ij} = \int_0^1 \int_0^1 \phi_i(x) \phi_j(y)\,dy\,dx,\\
V_{ij} = \int_0^1 \int_0^1 \phi_i''(x) \phi_j(y)\,dy\,dx,\\
K_{ij} = \int_0^1 \int_0^1 \phi_i'''(x) \phi_j(y)\,dy\,dx,\\
F_{ij} = f\int_0^1 \phi_i(x) \,dx.$$

The vector solution for each $\lambda$ input is denoted $C$. If we plot the magnitude of $C$ against several values of $\lambda$ we would arrive at a frequency response (see attached picture). The $\lambda$ values associated with the peaks are in fact the eigenvalues of the unforced problem.

Okay, now for my question: the $\lambda$ associated with each peak is real. However, I know the analytic solution (written above in decimal form) is complex. How would I separate the forced solutions of $\lambda$ to ascertain the real and imaginary components of the analytic eigenvalues?

EDIT: the preview does not display LaTeX form for me (google chrome, ubuntu 18.04). Also, the "\\" command no longer line breaks within latex.

 

[Paul Colby]

Neat. If there are reasons to assume the peaks are first order poles located just off the real axis, then you could peak fit the curves given,

$f(\lambda) = \left|\sum \frac{A_k}{\lambda-p_k}\right|^2$

extracting the magnitude of the imaginary pole offset, $\text{Im}(p_k)$. The sign of the imaginary components would likely have to come from some principle like causality in the case of particle physics for example.

 

[member 428835]

Neat. If there are reasons to assume the peaks are first order poles located just off the real axis, then you could peak fit the curves given,

$f(\lambda) = \left|\sum \frac{A_k}{\lambda-p_k}\right|^2$

extracting the magnitude of the imaginary pole offset, $\text{Im}(p_k)$. The sign of the imaginary components would likely have to come from some principle like causality in the case of particle physics for example.

Thanks for your reply! Can you elaborate what you mean? I'm lost. What is $A_k,p_k$? And what do you mean by "extracting the magnitude of the imaginary pole".

Also, I should add that the $\lambda$ value I get from the peak is precisely equal to the magnitude of the analytic eigenvalue, if this helps.

 

[Paul Colby]

Well, are you familiar with the method of least squares? Basically one takes data, in your case a function $f(\lambda)$, and fits it with a model. In your case this model is the absolute square of a sum of 3 first order pole terms. Each pole term has an amplitude, $A_k$, and a complex pole location, $p_k$. In rough terms, when the poles are close to the real axis the height of each peak yields the amplitude, $A_k$, while the peak width the distance from the real-$\lambda$ axis. Poles further from the real axis will be broader, closer will be narrower. Try a simple case by generating data with one pole.

 

[member 428835]

Well, are you familiar with the method of least squares? Basically one takes data, in your case a function $f(\lambda)$, and fits it with a model.

You're saying curve-fit the plot $|C|$ by $\lambda$ I've included, right?

In your case this model is the absolute square of a sum of 3 first order pole terms. Each pole term has an amplitude, $A_k$, and a complex pole location, $p_k$.

In the plot I've shown, everything is real, so I'm confused what $p_k$ is? Do you have a reference I could look at?

 

[Paul Colby]

You're saying curve-fit the plot |C| by λ I've included, right?

yes

In the plot I've shown, everything is real, so I'm confused what pk is? Do you have a reference I could look at?

I guess, not knowing anything about the problem beyond what you've written, that you are looking at what's called a resolvent, $R(\lambda)$. The resolvent for an operator, $A$, form this would look like,

$R(\lambda) = \|(A-\lambda I)^{-1}\|$

which becomes singular at an eigenvalue of $A$. In your case this is at a complex locations I've blessed with the label, $p_k$.

[edit] If you're looking for an analytic solution rather than numeric, my suggestion is of little value.

 

[member 428835]

yes
I guess, not knowing anything about the problem beyond what you've written, that you are looking at what's called a resolvent, $R(\lambda)$. The resolvent for an operator, $A$, form this would look like,

$R(\lambda) = \|(A-\lambda I)^{-1}\|$

which becomes singular at an eigenvalue of $A$. In your case this is at a complex locations I've blessed with the label, $p_k$.

[edit] If you're looking for an analytic solution rather than numeric, my suggestion is of little value.

Definitely not looking for an analytic solution, but I confess I'm very confused by what you're suggesting. I have $C(\lambda)$, where $\lambda$ is real and $C$ is complex valued. I'm confused what $A$ refers to here since this is not an eigenvalue problem anymore (and previously was a quadratic eigenvalue problem). And is $p_k$ the decimal complex eigenvalues I listed in post 1?

If you'd like, I can send you the MATLAB script with data so you can see exactly what is happening (pretty short, and takes about 1 min to run). I want to understand what you're getting at, but I'm very lost. Again, if you have any reference I'd be happy to read it.

 

[Paul Colby]

Do you have a reference I could look at?

I could reference my thesis

"An Application of Pole Extrapolation to the Mass Three Helium Isotope"
University of Wisconsin - Madison 1984

but it's an awful read. In it I was extracting poles terms from cross section measurement. In this case the data is actual data and is required to be real. The locations of the poles were determined by particle kinematics and I was after the pole term amplitude, $A_k$.

I want to understand what you're getting at, but I'm very lost.

$A$ is the operator which leads to the eigenvalue problem you posted. The driven problem, much like the resolvent, will have resonances located in the complex $\lambda$ plane. Your plot is a cut taken along the real axis.

So if you have the analytic solution, why bother? I assume because you are interested in extracting numbers when an analytic solution isn't available.

If you'd like, I can send you the MATLAB script with data so you can see exactly what is happening

I'd love to but I don't have Matlab.

 

[Paul Colby]

The vector solution for each λ input is denoted C.

When you generate the complex vector $C$, why must you use real $\lambda$? Any complex value should work, well except at exactly an eigenvalue that is.

 

[member 428835]

When you generate the complex vector $C$, why must you use real $\lambda$? Any complex value should work, well except at exactly an eigenvalue that is.

Okay, so you're suggesting taking a 3D image instead of the 2D slice? Interesting idea!

Also, attached is the problem broken down into better wording, for you viewing convenience (though I think you already understand what I'm talking about).

 

[member 428835]

I could reference my thesis

"An Application of Pole Extrapolation to the Mass Three Helium Isotope"
University of Wisconsin - Madison 1984

but it's an awful read. In it I was extracting poles terms from cross section measurement. In this case the data is actual data and is required to be real. The locations of the poles were determined by particle kinematics and I was after the pole term amplitude, $A_k$.

$A$ is the operator which leads to the eigenvalue problem you posted. The driven problem, much like the resolvent, will have resonances located in the complex $\lambda$ plane. Your plot is a cut taken along the real axis.

So if you have the analytic solution, why bother? I assume because you are interested in extracting numbers when an analytic solution isn't available.
I'd love to but I don't have Matlab.

I think I will check out your thesis. Can't hurt, and could possibly help. So thanks! And yep, this is a toy problem. The real one doesn't admit analytic solutions, sadly.

 

[S.G. Janssens]

I'd love to but I don't have Matlab.

I think that for this problem you could try Octave.
It is useful as a drop-in MATLAB replacement.

 

[Paul Colby]

Nice writeup thank. So, if one replaces the drive vector $F$ with a unit matrix, then you have upon inverting the LHS and taking the norm a generalized resolvent. I assume $F$ as the other matrices which appear are dictated by the real problem.

 

[member 428835]

Nice writeup thank. So, if one replaces the drive vector $F$ with a unit matrix, then you have upon inverting the LHS and taking the norm a generalized resolvent. I assume $F$ as the other matrices which appear are dictated by the real problem.

Thanks! Okay, so turns out your suggestion in post 9 did the trick! Thanks so much! I've included a plot, and you can clearly see the poles lie off the real axis. When I compare these results they match the analytic prediction! Just wanted to post to close the matter. Thanks a ton Paul!