3
$\begingroup$

Consider an ODE eigensystem $$ t(y+\frac{1}{s})a(y)+[(q+\frac{1}{2}+\frac{s}{2}y)+s(y\partial_y+\frac{1}{2})]b(y)=\lambda a(y)\\ t(y+\frac{1}{s})b(y)+[(q+\frac{1}{2}+\frac{s}{2}y)-s(y\partial_y+\frac{1}{2})]a(y)=\lambda b(y) $$ where $q,t,s$ are (small, if necessary) positive numbers. Assume Dirichlet boundary condition at the infinity. From the background of this problem, it only requires nondivergence at the singularity $y=0$ and vanishing at the infinity. We are mostly interested in a few, say, the first 8 smallest, eigenvalues.

  • I tried to reduce it to a 2nd-order ODE of $b$.

    This often solves the eigensystem more accurately. But unfortunately, it messes up the eigenstructure, i.e., $\lambda$ appears in various terms of the 2nd-order equation. And DSolve does not proceed, either.

  • Then I tried the following to solve the coupled ODE system directly.

For the parameters in the code below (same form, but slightly different parameters than the Latex equations above), I have some exact data from some other approach (not by solving ODE, irrelevant here, I think).
The goal here is actually to use ODE solutions to match these 8 eigenvalues near 0.

{-0.0811827, -0.0660165, -0.0462057, 0.00176381, 0.0051576, 0.0497911, 0.06966, 0.084884}

The code below gives

{-0.0829815, -0.06781, -0.0479909, -5.28964*10^-13, 5.28957*10^-13, 0.0480429, 0.0679139, 0.0831374}

There's some correspondence but certainly not quite accurate, especially the middle two close to zero.

a = 1/Sqrt[3]; s = 2*^-3; t = 3*^-5; eps0 = 1*^-10; Nband = 8;
yR = -200; yL = -480;
p = (1 - (a q)/2); v = a q; r = a (3/2 + (a q)/4); u = a (1 + (a q)/2);
Fop1[F_, pm_] := (v + s p (y + u/(r s))) F + 
   pm (s r (y D[F, y] + 1/2 F));
variables = {\[Alpha], \[Beta]};
lhs[q_] = {Fop1[\[Beta][y], 1] + t (y + u/(r s)) \[Alpha][y], 
   Fop1[\[Alpha][y], -1] + t (y + u/(r s)) \[Beta][y]};
bc = DirichletCondition[
   Table[component@y == 0, {component, variables}], True];
vals = NDEigenvalues[{lhs[0.0], 
     bc}, {\[Alpha][y], \[Beta][y]}, {y, yL, yR}, Nband, 
    Method -> {"SpatialDiscretization" -> {"FiniteElement", \
{"MeshOptions" -> {"MaxCellMeasure" -> 0.01, "MeshOrder" -> 2}}}}];
vals = Sort[vals, Re@#1 < Re@#2 &]

Since it is a seemingly simple ODE system, I was wondering if any other way could solve it better.

$\endgroup$
  • $\begingroup$ Comments are not for extended discussion; this conversation has been moved to chat. $\endgroup$ – Kuba Jan 30 at 9:54
-1
$\begingroup$

The eigenvalues appear as the condition of non-trivial solution for the linear system

$$ M(\lambda,q,t)\left(\begin{array}{c}\alpha\\ \beta\end{array}\right) = 0 $$

obtained after the series expansion proposed solution. The condition is

$$ \det M(\lambda,q,t) = 0 $$

which gives

$$ (\cdots)\left(\lambda ^2-(q-1) (q+2)\right) \left(\lambda ^2-(q-2) (q+3)\right) \left(\lambda ^2-(q-3) (q+4)\right) \left(\lambda ^2-(q-4) (q+5)\right) \left(\lambda ^2-(q-5) (q+6)\right) \left(\lambda ^2-(q-6) (q+7)\right) \left(\lambda ^2-(q-7) (q+8)\right) \left(\lambda ^2-(q-8) (q+9)\right) \left(\lambda ^2-(q-9) (q+10)\right) \left(\lambda ^2-(q-10) (q+11)\right)=0 $$

as can be extracted from the following script

Clear[q, t, y]
n = 10;
Op11 = (q + (y + 1)/2 + 1/2) # + y D[#, y] &;
Op12 = (q + (y + 1)/2 - 1/2) # - y D[#, y] &;
Op2 = t y # &;
ay = Sum[Subscript[alpha, k] y^k, {k, 0, n}];
by = Sum[Subscript[beta, k] y^k, {k, 0, n}];
eq10 = Op2[ay] + Op11[by] - lambda ay;
eq20 = Op2[by] + Op12[ay] - lambda by;
coef1 = CoefficientList[eq10, y];
coef2 = CoefficientList[eq20, y];

equs1 = Thread[coef1 == 0];
equs2 = Thread[coef2 == 0];
Ta = Table[Subscript[alpha, k], {k, 1, n}];
Tb = Table[Subscript[beta, k], {k, 1, n}];
vars = Join[Ta, Tb];
equs = Join[Take[coef1, {2, Length[coef1] - 1}], Take[coef2, {2, Length[coef2] - 1}]];
M = Grad[equs, vars];
M // MatrixForm
Det[M] // FullSimplify
M // Eigenvalues // MatrixForm

Curiously the eigenvalues appear non dependent on $t$

In the case of Laurent series the script is

Clear[q, t, y]
n = 10;
Op11 = (q + (y + 1)/2 + 1/2) # + y D[#, y] &;
Op12 = (q + (y + 1)/2 - 1/2) # - y D[#, y] &;
Op2 = (t y ) # &;
ay = Sum[Subscript[alpha, k] y^k, {k, -n, n}];
by = Sum[Subscript[beta, k] y^k, {k, -n, n}];
eq10 = (Op2[ay] + Op11[by] - lambda ay) y^n;
eq20 = (Op2[by] + Op12[ay] - lambda by) y^n;
coef1 = CoefficientList[eq10, y];
coef2 = CoefficientList[eq20, y];

equs1 = Thread[coef1 == 0];
equs2 = Thread[coef2 == 0];
equs12 = Join[equs1, equs2];
Ta = Table[Subscript[alpha, k], {k, -n, n}];
Tb = Table[Subscript[beta, k], {k, -n, n}];
vars = Join[Ta, Tb];
equs = Join[Take[coef1, {1, Length[coef1] - 1}], Take[coef2, {1, Length[coef2] - 1}]];
MatrixForm[equs12]
M = Grad[equs, vars];
M // MatrixForm
Det[M] // FullSimplify
M // Eigenvalues // MatrixForm
$\endgroup$
  • $\begingroup$ Will your series expansion for $a$ and $b$ not tend to infinity as $y->\infty$? And why do you not include the constant term in the expansion? $\endgroup$ – KraZug Jan 24 at 6:21
  • $\begingroup$ @KraZug My answer was for an OP that was edited. After this edition the operator looks with an additional parameter $s$. $\endgroup$ – Cesareo Jan 24 at 9:26
  • $\begingroup$ Yes, I note the stealth edit that changes the question significantly. But I think my initial questions still apply for the original equation, you are supposing a series expansion in powers of $y^k$, which are not going to meet the condition of $a(y) \to 0 $ and $y \to \infty$. $\endgroup$ – KraZug Jan 24 at 9:34
  • $\begingroup$ @KraZug Regarding the constant term now it is included. Thanks. $\endgroup$ – Cesareo Jan 24 at 9:35
  • $\begingroup$ @KraZug The behavior of $a(y), b(y)$ will be dictated by the series convergence. $\endgroup$ – Cesareo Jan 24 at 9:38

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.