7
$\begingroup$

I am trying to solve the eigensystem of a 1st-order linear ODE system in the region $(-\infty,\infty)$ and with Dirichlet boundary condition at the infinities $$ -i\partial_xu(x)+f^*(x)v(x)=\lambda u(x)\\ f(x)u(x)+i\partial_xv(x)=\lambda v(x) $$ where $f(x)$ is a complex-valued function mainly varying around $x=0$ and $f^*(x)$ is its complex conjugate. I define $f(x)$ by its norm and phase in the code, where I give two forms of $f(x)$ (using Tanh or Lorentzian) of the similar profile.
I tried reducing it to a 2nd-order ODE. But unfortunately, it messes up the eigenstructure, i.e., $\lambda$ appears in several terms. I have no idea if it can be solved analytically.

It is known that the system will have a few (at least one) discrete real eigenvalues in $(-1,1)$ and the eigenfunction is more or less localized around $x=0$. Outside $(-1,1)$, there will be a continuous spectrum with probably not quite confined eigenfunctions. I am interested in the eigenvalue(s) in $(-1,1)$ only.

The code below gives 4 solutions but some eigenfunctions are very saw-tooth like. Therefore, I want to make sure if these solutions are true or spurious and if any other way to solve it reliably.

{0.0856771, -0.169612, -0.881437, 0.94686}

a = 1; b = 0.8 a; c = 1; Nless = 8; cutoff = 20;
xR = cutoff; xL = -cutoff;
f[x_, pm_] = (a + b (Tanh[x/c - 1] - Tanh[x/c + 1])/(2 Tanh[1])) Exp[
    I pm \[Pi] (Tanh[x/c] + 1)/2];
(*f[x_, pm_] = (a + b (-c^2)/(x^2 + c^2)) Exp[
    I pm \[Pi] (Tanh[x/c] + 1)/2];*)
Fop1[F_] := D[F, x];
variables = {u, v};
lhs = {-I Fop1[u[x]] + f[x, 1] v[x], I Fop1[v[x]] + f[x, -1] u[x]};
bc = DirichletCondition[
   Table[component@x == 0, {component, variables}], True];
{vals, funcs} = 
  NDEigensystem[{lhs, bc}, {u[x], v[x]}, {x, xL, xR}, Nless, 
   Method -> {"PDEDiscretization" -> {"FiniteElement", {"MeshOptions" \
-> {"MaxCellMeasure" -> 0.001, "MeshOrder" -> 2}}}}];
vals
Table[Plot[
  If[iband <= Nless, Norm, _]@funcs[[iband, ab]], {x, xL, xR}, 
  PlotRange -> All, ImageSize -> Small], {ab, {1, 2}}, {iband, 1, 
  6}(*,{part,{Norm(*,Re,Im*)}}*)]
Improve this question
$\endgroup$
6
  • $\begingroup$ The short answer is yes, i think the 0.0856771 value is spurious and the other three are fine. I'll hopefully post an answer using my package later. $\endgroup$ – SPPearce Jul 12 '19 at 16:46
  • $\begingroup$ I don't really know why NDEigensystem gives these spurious roots. The Evans function has a double root at zero here, in a similar way to your previous question, that could be part of the reason. $\endgroup$ – SPPearce Jul 12 '19 at 19:25
  • $\begingroup$ @KraZug I notice a weird behaviour. Your package doesn’t give any solution if I change $\pi$ even by a tiny bit in the phase factor of $f$ while NDEigen still produces some solutions. Any remedy or reason? Thanks. $\endgroup$ – xiaohuamao Jul 12 '19 at 23:07
  • $\begingroup$ The Evans function becomes complex for real values of $\lambda$, so Plot won't show anything unless you use ReIm or Abs. Or use ContourPlot or just FindRoot. $\endgroup$ – SPPearce Jul 13 '19 at 7:23
  • $\begingroup$ @KraZug I just use FindRoot and it gives nothing when I set $\pi$ to $1.1\pi$ in the phase factor. Does your answer give any solution for this case? $\endgroup$ – xiaohuamao Jul 13 '19 at 7:33
7
$\begingroup$

I have a package for solving eigenvalue boundary value problems using the Compound Matrix Method with the Evans function, which I'll use here. The package is available on my GitHub (which has a notebook with examples), more details are in my other answers to questions here, and a good introduction to the method is in this pdf.

Install and load the package:

PacletInstall["CompoundMatrixMethod", 
     "Site" -> "http://raw.githubusercontent.com/paclets/Repository/master"];
Needs["CompoundMatrixMethod`"]

First write the system in the form $\mathbf{y}' = \mathbf{A} \cdot \mathbf{y}$ using the function ToMatrixSystem:

sys1 = ToMatrixSystem[Thread[lhs == λ {u[x], v[x]}], {}, {u, v}, {x, -L, L}, λ]

Note by not giving boundary conditions the requirement is that $\mathbf{A}$ tends to the decaying behaviour of $\mathbf{A}_{\pm \infty} =\lim_{x\to\pm\infty} \mathbf{A}(x)$. This will only work when there is a negative (resp. positive) eigenvalue of $\mathbf{A}_{\infty} (\mathbf{A}_{- \infty})$. Here the eigenvalues become imaginary for real (not complex) values outside (-1,1):

Plot[Evaluate[Eigenvalues[Limit[sys1[[1]], x -> ∞]]], {\[FormalLambda], -1, 1}]

enter image description here

and Evans will fail without explicit boundary conditions at those points. This method may be wrong to apply here, I think I need to consider further whether that is only applicable for higher order differential equations (e.g. a second order ODE in one variable).

The Evans function is an analytic function whose zeroes correspond to eigenvalues, so we can plot that and look for roots: 

Plot[Evans[λ, sys1 /. L -> 10], {λ, -0.99, 0.99}, 
 Epilog -> Point[{#, 0} & /@ vals]]

enter image description here

It looks initially like all 4 of the NDEigensystem values are roots, but if we zoom in we can see that there is a double root at zero, so the value at $0.0856771$ is spurious:

enter image description here

The winding number of the circle with radius 0.99 is 5, showing there are no further eigenvalues with $|\lambda|<=0.99$. The double root at zero is contributing two to the count here.

PlotEvansCircle[sys1 /. L -> 10, ContourRadius -> 0.99, nPoints -> 1000, Joined -> True]

enter image description here

Note that providing explicit boundary conditions in gives very similar results, except for the eigenvalue at zero.

Improve this answer
$\endgroup$
8
  • $\begingroup$ Sorry, but it seems that your links point to wrong content. Could you please check them again? $\endgroup$ – Pinti Jul 15 '19 at 6:59
  • $\begingroup$ @Pinti Ah sorry, I copied those from a previous post but forgot to include the links. Duh. Edited them. $\endgroup$ – SPPearce Jul 15 '19 at 8:59
  • $\begingroup$ @Pinti, if you have any questions regarding my package I'm more than happy to discuss them. $\endgroup$ – SPPearce Jul 15 '19 at 9:56
  • $\begingroup$ @KraZug, thank you for the provided link. I read your answers and the Introduction by Fu to the method, which is useful! However, I have difficulty to understand Eq.(16) in that Introduction. As can be seen after integration the argument of the exponential does depend on $d$ since it is the integral limit, so why Evans function $D(\lambda)$ is independent of $d$? Btw, do you have the ref. Chadwick(1976) in that document, as I tried to read it, I found no item in the Reference... $\endgroup$ – Nobody Mar 8 '20 at 5:02
  • $\begingroup$ @Nobody, I'm not entirely sure what that Chadwick reference is. It'll definitely be Peter Chadwick, might be his book on elasticity but with the wrong year. The reference is just being used for the general matrix/tensor inequality, rather than Evans function specific. But if you take the equation (16) and differentiate with respect to 𝑑, using the formula given in (17) and (18), then you can show that $\partial D/\partial d = 0$, and hence $D(\lambda)$ is independent of the matching point $d$. $\endgroup$ – SPPearce Mar 10 '20 at 6:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.