# Using NDEigensystem to solve coupled eigenvalue problem

I want to find the Eigenvalues and Eigenfunctions for the following Eigenvalue problem I tried to solve this numerically

pars=Rationalize[{L-> 6.9,A->0.0133,a0->410,
c1a->0.003,c2a->0.002,c1p->-2.5 10^-6,c2p->4.3 10^-6}];


I defined to Operator like above

op={-(a0^2/A)D[\[Phi]2[xi],xi], -A D[\[Phi]1[xi],xi]}/.pars


I am not shure how to use the

DirichletCondition

but I tired this

bc={DirichletCondition[\[Phi]2[xi]==c1p \[Phi]1[xi],xi==0],
DirichletCondition[\[Phi]2[xi]==c2p \[Phi]1[xi],xi==L]}/.pars


Unfortunately Mathematica could not understand my input

NDEigensystem[Flatten[{op,bc}],{\[Phi]1[xi],\[Phi]2[xi]},{xi,0,L}/.pars,4] Did I make a mistake?

Thanks!

• The missing message text is: NDEigenvalues::fembdcc: Cross-coupling of dependent variables in DirichletCondition[\[Phi]2==-(\[Phi]1/400000),xi==0] is not supported in this version. Aug 22, 2017 at 14:04
• Thanks! So I can't solve it with this version, at all? Aug 22, 2017 at 14:54
• NDEigensystem can not handle inhomogeneous DirichletConditions and thus it can not handle cross coupling of boundary conditions. You could try to manually make a solution. But the coupled ODEs are also convection dominant which makes it harder. Do you have an expected solution? Aug 22, 2017 at 16:04

I know this question is old and has an explicit analytic solution given, but here is a numerical method to confirm the roots, which works for cross-coupling of the boundary conditions as seen here.

I have a package for numerically calculating solutions of eigenvalue problems using the Evans function via the method of compound matrices, which is hosted on github. See my answers to other questions or the github for some more details.

First we install the package (only need to do this the first time):

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


Then we first need to turn the ODEs into a matrix form $$\mathbf{y}'=\mathbf{A} \cdot \mathbf{y}$$, using my function ToMatrixSystem:

Needs["CompoundMatrixMethod"]

pars = Rationalize[{L -> 6.9, A -> 0.0133, a0 -> 410, c1a -> 0.003,
c2a -> 0.002, c1p -> -2.5 10^-6, c2p -> 4.3 10^-6}];
eqns = {-(a0^2/A) D[ϕ2[xi], xi] == λ ϕ1[xi], -A D[ϕ1[xi], xi] == λ ϕ2[xi]}
bcs = {ϕ2 == c1p ϕ1, ϕ2[L] == c2p ϕ1[L]}
sys = ToMatrixSystem[eqns, bcs, {ϕ1, ϕ2}, {xi, 0, L}, λ] /. pars


The object sys contains the matrix $$\mathbf{A}$$, as well as similar matrices for the boundary conditions and the range of integration.

Now the function Evans will calculate the Evans function (also known as the Miss-Distance function) for any given value of $$\lambda$$; this is an analytic function whose roots coincide with eigenvalues of the original equation.

FindRoot will then find solutions for a given start point:

FindRoot[Evans[λ, sys], {λ, 1}]
(* {λ -> -12.5116} *)


And this gives the same values as for the analytic solution given in the other answer.

λ/.FindRoot[Evans[λ, sys], {λ, -12 + # I}] & /@ Range[-1860, 1860, 186] // Quiet
(* {-12.5116 - 1866.74 I, -12.5116 - 1680.07 I, -12.5116 - 1493.39 I, -12.5116 - 1306.72 I,
-12.5116 - 1120.05 I, -12.5116 - 933.372 I, -12.5116 - 746.697 I, -12.5116 - 560.023 I,
-12.5116 - 373.349 I, -12.5116 - 186.674 I, -12.5116, -12.5116 + 186.674 I,
-12.5116 + 373.349 I, -12.5116 + 560.023 I, -12.5116 + 746.697 I, -12.5116 + 933.372 I,
-12.5116 + 1120.05 I, -12.5116 + 1306.72 I, -12.5116 + 1493.39 I, -12.5116 + 1680.07 I,
-12.5116 + 1866.74 I} *)

• for the coupled odes, if an initial guess of the eigval is not available, is it possible to replace the FindRoot with NSolve, for example, \[Lambda] /. NSolve[Evans[\[Lambda], sys] == 0]. Or is there any way to obtain all roots of Evans function without an initial guess? Thank you. Dec 23, 2021 at 6:22
• I don't think that you can find all roots easily I'm afraid, I don't think it will work in NSolve. There is a function to calculate the Winding Number of a contour, which tells you how many roots lie within a circle, although be slightly careful that you have enough points when using that. Dec 24, 2021 at 8:34

I have computed the analytic solution of the Eigenvalue problem above

With the Eigenvalues

\[Beta] -> Log[(Sqrt[A + a0 c1p] Sqrt[A - a0 c2p])/Sqrt[(A - a0 c1p) (A + a0 c2p)]]

ev[k_] := (a0 (I k \[Pi] + \[Beta]))/L


with some constant beta and the Eigenfunctions

ef[xi_,k_] := {{Cosh[(xi ev[k])/a0]/c1p-(a0 Sinh[(xi ev[k])/a0])/A},
{Cosh[(xi ev[k])/a0]-(A Sinh[(xi ev[k])/a0])/(a0 c1p)}}