# Numerical eigenvalues of the problem

I have this equations with condition 0 <= x <= 1:

eqns = {((2 - x) y''[x])'' + λ y''[x] == 0, y[0] == 0, y[1] == 0, y'[0] == 0, y''[1] == 0}


So I need to find the first 4 eigenvalues here. I tried to use NDEigenvalues and get this error NDEigenvalues::femcmsd: The spatial derivative order of the PDE may not exceed two. How can I avoid it? I got answer when used z[x] = y''[x]: {1.72307*10^-13, -15.0124, -58.1943, -130.156}, but I'm not sure that I'm right.

• This piece of code ((2 - x) y''[x])'' is not going to work. May 25, 2020 at 14:59
• I can replace it with [D[(2 - x) y''[x], x], x] or -2 y'''[x] + (2 - x) y''''[x] May 25, 2020 at 15:03

Although there are a number of related ways to solve this problem, the cleanest appears be replacing y''[x] by z[x] to obtain

spp = DSolveValue[{D[(2 - x) z[x], {x, 2}] + λ z[x] == 0}, z[x], x]
(* (Sqrt[(-2 + x) λ] BesselI[1, 2 Sqrt[(-2 + x) λ]] C[1])/(2 - x)
- (Sqrt[(-2 + x) λ] BesselK[1, 2 Sqrt[(-2 + x) λ]] C[2])/(2 - x) *)


Then, Integrate to obtain y'[x] and y[x] in turn.

sp = Integrate[spp, x] + C[3]
(* C[1] - BesselI[0, 2 Sqrt[(-2 + x) λ]] C[1]
- BesselK[0, 2 Sqrt[(-2 + x) λ]] C[2] + C[3] *)
s = Simplify@Integrate[sp, x] + C[4]
(* (Sqrt[(-2 + x) λ] BesselK[1, 2 Sqrt[(-2 + x) λ]] C[2])/λ + x (C[1] + C[3]) + C[4]
- (-2 + x) C[1] Hypergeometric0F1Regularized[2, (-2 + x) λ] *)


and compute the determinant of the boundary conditions.

CoefficientArrays[{spp /. x -> 1, sp /. x -> 0, s /. x -> 0, s /. x -> 1},
{C[1], C[2], C[3], C[4]}] // Normal // Last;
disp = Det[%] // FullSimplify
(* 1/2 π Sqrt[λ] (-BesselY[1, 2 Sqrt[λ]] (Hypergeometric0F1Regularized[1, -2 λ]
- 2 Hypergeometric0F1Regularized[2, -2 λ])
+ (Sqrt[λ] BesselY[0, 2 Sqrt[2] Sqrt[λ]] - Sqrt[2] BesselY[1, 2 Sqrt[2] Sqrt[λ]])
Hypergeometric0F1Regularized[2, -λ]) *)


Plotting this function shows its zeroes.

Plot[Chop@disp, {λ, -5, 250}, AxesLabel -> {λ, d}, LabelStyle -> {15, Bold, Black}]


(disp becomes exponentially large at negative λ.) The first four zeroes are 0 and

{FindRoot[disp, {λ, 30}], FindRoot[disp, {λ, 90}], FindRoot[disp, {λ, 170}]}
// Flatten // Chop // Values
(* {29.4788, 87.0186, 173.309} *)


Even more compact is

disp = FullSimplify@Det[Normal@Last@CoefficientArrays[
DSolveValue[{D[(2 - x) y''[x], {x, 2}] + λ y''[x] == 0},
{y[0], y[1], y'[0], y''[1]}, x], {C[1], C[2], C[3], C[4]}]];
Table[FindRoot[disp, {λ, λ0}], {λ0, {30, 90, 170}}] // Flatten // Values

• Thank you very much! I'm really glad to see this full answer. May 25, 2020 at 22:36
• I was impeded by the word "eigenvalues". May 26, 2020 at 8:44

The answer by bbgodfrey is excellent for this particularly problem, because it can be solved directly by DSolve.

Having said that, I have a package which helps solve eigenvalue BVPs by calculating the Evans function, an analytic function whose roots correspond to the eigenvalues. Some details are available at these two questions, or this PDF. Or search for CompoundMatrixMethod to see my previous answers here.

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[D[(2 - x) y''[x], x, x] + λ y''[x] == 0,
{y[0] == 0, y[1] == 0, y'[0] == 0, y''[1] == 0}, y, {x, 0, 1}, λ]


Now we can evaluate the Evans function at a given value of $$\lambda$$ with e.g. Evans[1, sys]. Roots of this function correspond to eigenvalues of the original equation:

Plot[Evans[λ, sys], {λ, 0, 250}]


We can see that this corresponds to the eigenvalues shown in bbgodfrey's answer, although without the $$\lambda = 0$$ root, which isn't a real eigenvalue.

KraZug's entirely numerical solution reminded me of another approach, employing a variation on an example described in the NDSolve documentation, Boundary Value Problems with Parameters. However, the documented approach, which entails treating the parameter as a variable to be computed, has only one solution for the parameter, and the example equation has an adequate number of boundary conditions. Here, the parameter λ can assume an infinite number of values, and treating it as a variable requires an additional boundary condition. Not coincidentally, however, the normalization of the eigenfunctins is arbitrary, and an additional boundary condition for y merely specifies the normalization, which subsequently can be changed (as done below with the variable norm), if desired. And, different values for λ can be sought by providing different initial guesses for it. With this as preface,

s = ParametricNDSolveValue[{D[(2 - x) y''[x], {x, 2}] + λ[x] y''[x] == 0, λ'[x] == 0,
y[0] == 0, y[1] == 0, y'[0] == 0, y''[0] == 1, y''[1] == 0}, {y[x], λ[0]}, {x, 0, 1},
{λ0}, Method -> {"Shooting", "StartingInitialConditions" -> {λ[0] == λ0}}];


The extra boundary condition is y''[0] == 1. I also tried y'[1] == 1, but it works much less well, because the "Shooting" Method then must deal with three unknowns instead of two. The first ten eigenvalues and their associated eigenfunctions then are obtained by

Transpose@Table[s[n], {n, {30, 90, 170, 300, 450, 600, 770, 1050, 1300, 1600}}];
norm = 1/Sqrt@NIntegrate[First[%]^2, {x, 0, 1}];
Plot[Evaluate[norm First[%%]], {x, 0, 1}, ImageSize -> Large, AxesLabel -> {y, x},
LabelStyle -> {15, Bold, Black}, PlotLegends -> Placed[Automatic, {1, .5}]]
Last[%%%]


{29.4788, 87.0186, 173.309, 288.359, 432.171, 604.744, 806.079, 1036.18, 1295.04, 1582.66}

• That is nice! Gets you the eigenfunctions immediately too, which my package doesn't do without some work. Jun 5, 2020 at 14:52