# NDEigenvalues vs. FindRoot for finding the eigenvalues of Airy equation?

Say I am trying to find the first 5 eigenvalues of the differential equation $$f''(x)=\lambda x f(x)$$, on the interval [-1,0], with boundary conditions $$f(-1)=f(0)=0$$.

I will try to do this 3 ways, and compare them to values for the eigenvalues found via a WKB approximation.

Trying to solve this using NDEigenvalues without Dirichlet condition.

NDEigenvalues[{f''[x]/x}, f[x], {x, -1, 0}, 5]
{0., 25.6383, 95.9537, 210.72, 370.024}


Trying to solve this using NDEigenvalues with Dirichlet condition.

NDEigenvalues[{f''[x]/x, DirichletCondition[f[x] == 0, True]}, f[x], {x, -1, 0}, 5]
{0., 19.6448, 84.2639, 194.087, 349.122}


Trying to solve this using FindRoot.

First solve the differential equation analytically.

DSolve[{y''[x] == x A y[x]}, y[x], x]
{{y[x] -> AiryAi[A^(1/3) x] C + AiryBi[A^(1/3) x] C}}


The solution is in terms of Airy functions.

Now plug in the boundary conditions at $$x=0$$, and $$x=-1$$.

$$f(0)=f(-1)\\ \rightarrow c_{1}\mathrm{Ai}(0)+c_{2}\mathrm{Bi}(0) = c_{1}\mathrm{Ai}(-A^{1/3})-c_{2}\mathrm{Bi}(-A^{1/3})\\ \rightarrow c_{1}\mathrm{Ai}(0)+c_{2}\mathrm{Bi}(0)-c_{1}\mathrm{Ai}(-A^{1/3})+c_{2}\mathrm{Bi}(-A^{1/3})=0$$

Choose the constants to be $$1$$.

If FindRoot is used now.

FindRoot[AiryAi + AiryBi - AiryAi[-A^(1/3)] - AiryBi[-A^(1/3)], {A, {0, 20, 85, 200, 350}}, WorkingPrecision -> 10]
{A -> {0, 0.2867307855, 88.39876753, 798.9989135, 354.8666727}}


And there is an error.

FindRoot:: Encountered a singular Jacobian at the point {A} = {0,0.2867307855,88.39876753,798.9989135,354.8666727}. Try perturbing the initial point(s).


This result isn't too surprising, as the Airy functions are not very nice.

If we try to clean up the equation that is being put in FindRoot.

FindRoot[AiryBi[-A^(1/3)]/AiryAi[-A^(1/3)] - Sqrt, {A, {0, 20, 85, 200, 350}}, WorkingPrecision -> 10]
{A -> {0, 18.95626559, 81.88658338, 189.2209333, 340.9669591}}


WKB Approximation

A quick WKB approximation will tell you $$A=\left(\frac{3 n \pi}{2}\right)^{2}$$.

n = {0, 1, 2, 3, 4};
N[((3 n \[Pi])/2)^2]
{0., 22.2066, 88.8264, 199.859, 355.306}


## Questions

Is one of these approaches more correct than the others? Comparing the results of the NDEigenvalues to FindRoot, the FindRoot solutions are closer to the WKB approximation, especially at larger eigenvalues (e.g. 20 eigenvalues out).

Which version of the FindRoot method is more correct? It is possible I am misusing/misunderstanding some things, so if anyone can point out any mistakes, optimizations, or even other methods, that would be very helpful. Thank You.

• For the first FindRoot, the BCs are $f(-1)=f(0)=0$, not $f(-1)=f(0)$. You cannot choose both constants $c_1,c_2$ freely. – Michael E2 Oct 5 '18 at 10:45
• Yes, thank you. I knew you couldn't really do that, but was interested in the results it gave. – guest84924657 Oct 5 '18 at 19:33

## 5 Answers

From Mariusz's answer, I managed to figure out a slightly less ad hoc route.

First, solve the ODE sans boundary conditions:

DSolveValue[y''[x] == x λ y[x], y[x], x]
AiryAi[x λ^(1/3)] C + AiryBi[x λ^(1/3)] C


We want to find values of C and C so that the boundary condition $$y(-1)=y(0)=0$$ is satisfied. This can be recast as a matrix eigenvalue problem:

cmat = {{AiryAi[x λ^(1/3)], AiryBi[x λ^(1/3)]} /. x -> -1,
{AiryAi[x λ^(1/3)], AiryBi[x λ^(1/3)]} /. x -> 0};


Take the determinant and equate to zero:

Det[cmat] // Simplify
(3^(5/6) AiryAi[-λ^(1/3)] - 3^(1/3) AiryBi[-λ^(1/3)])/(3 Gamma[2/3])


which can be simplified to Sqrt AiryAi[-λ^(1/3)] - AiryBi[-λ^(1/3)] == 0, or

Det[cmat] // FullSimplify
(λ^(1/3) Hypergeometric0F1[4/3, -λ/9])/π


Finding C and C corresponds to getting the eigenvector of cmat corresponding to its zero eigenvalue:

FullSimplify[Eigensystem[cmat /. AiryBi[-λ^(1/3)] -> Sqrt AiryAi[-λ^(1/3)]]]
{{0, AiryAi[-λ^(1/3)] + 1/(3^(1/6) Gamma[2/3])},
{{-Sqrt, 1}, {3^(2/3) AiryAi[-λ^(1/3)] Gamma[2/3], 1}}}


which means the solution wanted is

%[[2, 1]].{AiryAi[x λ^(1/3)], AiryBi[x λ^(1/3)]}
-Sqrt AiryAi[x λ^(1/3)] + AiryBi[x λ^(1/3)]


As I noted in a comment, one can use Plot[] to get roots of transcendental equations (see this thread as well).

Applied to this problem:

With[{λmax = 9000},
eigs = Sort[Cases[Normal[Plot[λ Hypergeometric0F1[4/3, -λ^3/9],
{λ, 0, Surd[λmax, 3]}, Mesh -> {{0}},
MeshFunctions -> {#2 &}, MeshStyle -> {}]],
Point[{x_, _}] :>
(λ /. FindRoot[λ Hypergeometric0F1[4/3, -λ^3/9], {λ, x}]), ∞]^3]]
{18.956265591373203, 81.88658337813662, 189.22093329303374, 340.9669590647526,
537.1257454350938, 777.6975694239662, 1062.6825270308009, 1392.0806585113041,
1765.8919831156077, 2184.116510979998, 2646.754247850686, 3153.8051971811865,
3705.269361147913, 4301.146741177415, 4941.43733823621, 5626.141152997897,
6355.258185943543, 7128.788437424108, 7946.731907700498, 8809.088596969903}


Visualize the first 7:

   Plot[Table[AiryBi[x λ^(1/3)] - Sqrt AiryAi[x λ^(1/3)], {λ, Take[eigs, 7]}],
{x, -1, 0}, Evaluated -> True] As it turns out, the eigenvalues are expressible in terms of BesselJZero[] (via the relationship between Hypergeometric0F1[] and BesselJ[]):

eigs = N[Table[9 BesselJZero[1/3, k]^2/4, {k, 20}]]
{18.9562655913732, 81.88658337813663, 189.22093329303374, 340.96695906475253,
537.1257454350937, 777.6975694239661, 1062.6825270308009, 1392.080658511304,
1765.8919831156077, 2184.1165109799977, 2646.7542478506866, 3153.8051971811865,
3705.269361147912, 4301.146741177416, 4941.437338236211, 5626.141152997898,
6355.258185943542, 7128.788437424109, 7946.731907700495, 8809.088596969903}


Thus, we can define the eigenfunctions as

efuns[k_Integer, x_] := With[{l = (3 BesselJZero[1/3, k]/2)^(2/3)},
AiryBi[l x] - Sqrt AiryAi[l x]]

• You could also just use NSolve: NSolve[λ^(1/3) Hypergeometric0F1[4/3, -λ/9]==0 && 0<λ<10000, λ] – Carl Woll Oct 5 '18 at 14:58

You ask for other methods of finding eigenvalues, so I'll mention my package, which calculates the Evans function, an analytic function whose roots correspond to the eigenvalues. Some details are available at these two questions, or this PDF.

Install the package (also available on my github page):

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


Load the package and setup the system:

Needs["CompoundMatrixMethod"]
sys = ToMatrixSystem[f''[x] == λ x f[x], {f[-1] == 0, f == 0}, f, {x, -1, 0}, λ]


We can plot the Evans function:

Plot[Evans[λ, sys], {λ, 0, 2000}] The findAllRoots function given by Jens here is pretty good at finding all the roots, which are the same as those given by DSolve in another answer.

roots = findAllRoots[Evans[λ, sys], {λ, 0, 10000}]
(* {18.9563, 81.8866, 189.221, 340.967, 537.126, 777.698, 1062.68,
1392.08, 1765.89, 2184.12, 2646.75, 3153.81, 3705.27, 4301.15,
4941.44, 5626.14, 6355.26, 7128.79, 7946.73, 8809.09, 9715.86} *)


Higher precision can be achieved via WorkingPrecision, e.g:

FindRoot[Evans[λ, sys, WorkingPrecision -> 30], {λ, 80}, WorkingPrecision -> 30]
{λ -> 81.8865834393931580592680765082}


You can then plot the first 7 eigenfunctions (different normalisation to those in the other examples):

 Plot[Evaluate[Table[NDSolveValue[{f''[x] == λ x f[x], f[-1] == 0,
f == -1}, f, {x, -1, 0}][x], {λ, roots[[1 ;; 7]]}]], {x, -1, 0}] One advantage of going numerically means that you can go to equations that DSolve can't solve, as well as being more numerically stable than finding roots of a determinant which is essentially what you are trying to do in your code (which can also give you spurious solutions where the analytic solution has repeated roots for instance).

I'm more than happy to talk further about the method and package if you want.

Only long comment.I have no answer to all your questions, but code below work for me.

By DSolve:

sol = y[x] /. DSolve[{y''[x] == x*λ*y[x], y[-1] == 0, y == 0}, y[x],
x, Assumptions -> {0 < λ < 9000}][];
eigvals = {ToRules[sol[[1, 1, 2]]]} // N

(*{{λ -> 18.9563}, {λ -> 81.8866},
{λ -> 189.221}, {λ -> 340.967},
{λ -> 537.126}, {λ -> 777.698},
{λ -> 1062.68}, {λ -> 1392.08},
{λ -> 1765.89}, {λ -> 2184.12},
{λ -> 2646.75}, {λ -> 3153.81},
{λ -> 3705.27}, {λ -> 4301.15},
{λ -> 4941.44}, {λ -> 5626.14},
{λ -> 6355.26}, {λ -> 7128.79},
{λ -> 7946.73}, {λ -> 8809.09}} *)


By NDEigenvalues:

NDEigenvalues[{f''[x]/x, DirichletCondition[f[x] == 0, x == -1],
DirichletCondition[f[x] == 0, x == 0]}, f[x], {x, -1, 0}, 21,
Method -> {"PDEDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 0.00001}}}]

(*Try: "MaxCellMeasure" -> 0.000001}*)

(* {-1.14445,
18.9564, 81.887,
189.222, 340.968,
537.128, 777.7,
1062.69, 1392.09,
1765.9, 2184.12,
2646.76, 3153.81,
3705.28, 4301.16,
4941.45, 5626.16,
6355.27, 7128.81,
7946.75, 8809.11} *)

• Based on this result, here is one way to get a pile of eigenvalues: With[{λmax = 9000}, Sort[Cases[Normal[Plot[Sqrt AiryAi[-λ] - AiryBi[-λ], {λ, 0, Surd[λmax, 3]}, Mesh -> {{0}}, MeshFunctions -> {#2 &}, MeshStyle -> {}]], Point[{x_, _}] :> (λ /. FindRoot[Sqrt AiryAi[-λ] - AiryBi[-λ], {λ, x}]), ∞]^3] – J. M.'s ennui Oct 5 '18 at 10:45
• Thanks a lot :) – Mariusz Iwaniuk Oct 5 '18 at 10:50
• The equation for the eigenvalues can also be expressed as λ^(1/3) Hypergeometric0F1[4/3, -λ/9] == 0 – J. M.'s ennui Oct 5 '18 at 10:51
• You can use: neigvals = NDEigenvalues[{f''[x]/x, DirichletCondition[f[x] == 0, True]}, f[x], {x, -1, 0}, 20, Method -> {"PDEDiscretization" -> {"FiniteElement", "MeshOptions" -> {"MaxCellMeasure" -> 0.00001}}, "Eigensystem" -> {"Arnoldi", "Criteria" -> "RealPart"} }]; to get rid of the spurious first eigenvalue. – user21 May 27 '20 at 6:00

Here's a spectral method based on another Airy equation problem, How to solve ODE with boundary at infinity.

With[{n = 128, prec = MachinePrecision}, (* incr prec above 16 for greater acc *)
tvec = -Rescale@N[Cos[Range[0, n] Pi/n], prec]; (* so t=Cos[theta] incr -1 to 0 *)
(* derivative operators *)
dm2 = NDSolveFiniteDifferenceDerivative[2, tvec,
"DifferenceOrder" -> "Pseudospectral"]["DifferentiationMatrix"];
(* boundary bordering *)
dm2 = dm2[[2 ;; -2, 2 ;; -2]];
tvec = Reverse@tvec[[2 ;; -2]]; (* reverse t so theta is increasing *)
(* construct & solve linear system *)
opL = dm2/tvec;
]

eigs = Eigenvalues[opL, -50]
(*
{55331.9, 53137.2, 50986.8, 48880.9, 46819.4, 44802.3, 42829.6, \
40901.3, 39017.5, 37178., 35383., 33632.4, 31926.2, 30264.4, 28647., \
27074., 25545.5, 24061.3, 22621.6, 21226.3, 19875.4, 18568.9, \
17306.8, 16089.2, 14915.9, 13787.1, 12702.6, 11662.6, 10667., \
9715.86, 8809.09, 7946.73, 7128.79, 6355.26, 5626.14, 4941.44, \
4301.15, 3705.27, 3153.81, 2646.75, 2184.12, 1765.89, 1392.08, \
1062.68, 777.698, 537.126, 340.967, 189.221, 81.8866, 18.9563}
*)


Check eigenvalues:

sols = First@
NDSolve[{f''[x] == #*x f[x], f[-1] == 0, f'[-1] == 1},
f, {x, -1, 0}, WorkingPrecision -> Precision@eigs,
MaxSteps -> 100000] & /@ eigs;
f /. sols // ListPlot  (* check BC at x == 0 *) ListLinePlot[f /. sols[[-5 ;;]]] Comments on the OP's methods.

The NDEigensystem seems to shift the BC at x == 0 to the next mesh point, probably because of the x in the denominator of the differential operator.

In the first FindRoot, the BCs should be $$f(-1) = f(0) = 0$$ and only one of the constants $$c_1, c_2$$ may be freely chosen. The following yields accurate eigenvalues:

A /. FindRoot[
{AiryAi + k*AiryBi == 0, AiryAi[-A^(1/3)] + k*AiryBi[-A^(1/3)] == 0},
{A, {20, 85, 200, 350}}, {k, -{1, 1, 1, 1}/2}]


The second FindRoot works fine.

I don't like Airy fns, so I convert them with the following.

AiryToBesP = {AiryAi[x_] -> (1/3)*Sqrt[x]*(BesselI[-3^(-1), (2/3)*x^(3/2)] -
BesselI[1/3, (2/3)*x^(3/2)]), AiryBi[x_] ->
Sqrt[x/3]*(BesselI[-3^(-1), (2/3)*x^(3/2)] + BesselI[1/3, (2/3)*x^(3/2)])}
AiryToBesM = {AiryAi[x_] -> (1/3)*Sqrt[-x]*(BesselJ[-3^(-1), (2/3)*Sqrt[-x]^3] +
BesselJ[1/3, (2/3)*Sqrt[-x]^3]), AiryBi[x_] ->
Sqrt[-x/3]*(BesselJ[-3^(-1), (2/3)*Sqrt[-x]^3] - BesselJ[1/3, (2/3)*Sqrt[-x]^3])}


where the AiryToBesP is for x > 0 and AiryToBesM is for x < 0. We have the ode

ode = y''[x] == x A y[x];

DSolve[ode, y[x], x] // Flatten
(*{y[x] -> C AiryAi[A^(1/3) x] + C AiryBi[A^(1/3) x]}*)

y[x_] = y[x] /. % /. {C -> c1, C -> c2};


Given that x < 0 in this case.

y[x_] = y[x] /. AiryToBesM // Simplify
(*1/3 Sqrt[-A^(1/3) x] ((c1 + Sqrt c2) BesselJ[-(1/3), 2/3 (-A^(1/3) x)^(3/2)] + (c1 - Sqrt c2) BesselJ[1/3,
2/3 (-A^(1/3) x)^(3/2)])*)


Simplify the constants.

y[x_] = y[x] /. {c1 + Sqrt c2 -> c1, c1 - Sqrt c2 -> c2}
(*1/3 Sqrt[-A^(1/3) x] (c1 BesselJ[-(1/3), 2/3 (-A^(1/3) x)^(3/2)] +
c2 BesselJ[1/3, 2/3 (-A^(1/3) x)^(3/2)])*)


Check that the converted y satisfies the ode.

ode // FullSimplify
(*True*)


BesselJ[-1/3, x] is unbounded at x = 0 so we can set.

c1 = 0;


Now

y[-1] == 0
(*1/3 A^(1/6) c2 BesselJ[1/3, (2 Sqrt[A])/3] == 0*)


A in terms of a Bessel root.

Solve[(2*Sqrt[A])/3 == root, A]
(*{{A -> (9*root^2)/4}}*)


Now find the eigenvalues from BesselJZero

Table[9/4 BesselJZero[1/3, k]^2, {k, 1, 10}] // N
(*{18.9563, 81.8866, 189.221, 340.967, 537.126, 777.698, 1062.68, 1392.08, 1765.89, 2184.12}*)
`