# 1st-order linear ODE system gives inaccurate/biased solutions

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.

• Comments are not for extended discussion; this conversation has been moved to chat.
– Kuba
Jan 30 '19 at 9:54

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];

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 // 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];

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 // MatrixForm
Det[M] // FullSimplify
M // Eigenvalues // MatrixForm

• 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? Jan 24 '19 at 6:21
• @KraZug My answer was for an OP that was edited. After this edition the operator looks with an additional parameter $s$. Jan 24 '19 at 9:26
• 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$. Jan 24 '19 at 9:34
• @KraZug Regarding the constant term now it is included. Thanks. Jan 24 '19 at 9:35
• @KraZug The behavior of $a(y), b(y)$ will be dictated by the series convergence. Jan 24 '19 at 9:38