Eigenvalues of a fourth-order ODE

Question

Consider the following ODE for $y(x)$ over $x\in\left[0,\frac{1}{2}\right]$ with an eigenvalue $\lambda$

$\qquad 2x\,y''''+ 4y'''=\lambda\, y''$

The boundary conditions at $x=\frac{1}{2}$ are $y'\left(\frac{1}{2}\right)=0$ and $y''\left(\frac{1}{2}\right)=\lambda y\left(\frac{1}{2}\right)$.

At $x=0$ we enforce $y(0)=0$ and require $y'(0)$ to be a finite number.

How can I solve numerically for the eigenvalues (values of $\lambda$ corresponding to non-zero solutions) using Mathematica?

Answer 1

This problem can be solved symbolically, although perhaps not with DEigensystem. Instead, begin with DSolve.

s = DSolveValue[{2 x y''''[x] + 4 y'''[x] == lamda y''[x], y[0] == 0, 
    y'[1/2] == 0, y''[1/2] == lamda y[1/2]}, y[x], x] // FullSimplify
(* (16 C[2] (BesselI[0, Sqrt[lamda]] (-1 + 2 x + Sqrt[2] Sqrt[lamda] Sqrt[x]
    BesselK[1, Sqrt[2] Sqrt[lamda] Sqrt[x]]) + x (-2 + lamda BesselK[0, Sqrt[lamda]])
    Hypergeometric0F1Regularized[2, (lamda x)/2]))/(lamda^2 ((8 + lamda) 
    Hypergeometric0F1Regularized[3, lamda/4] + lamda Hypergeometric0F1Regularized[4, lamda/4])) *)

Next, expand the derivative of this expression about x = 0 to determine the eigenvalues.

ser = (Normal@Series[D[s, x], {x, 0, 0}]) // FullSimplify
(* (8 C[2] (-4 + 2 lamda BesselK[0, Sqrt[lamda]] + BesselI[0, Sqrt[lamda]] 
   (4 + 2 EulerGamma lamda + lamda (Log[lamda/2] + Log[x]))))/
   (lamda^2 ((8 + lamda) Hypergeometric0F1Regularized[3, lamda/4] + 
   lamda Hypergeometric0F1Regularized[4, lamda/4])) *)

The Log[x] term will not be finite at x = 0 unless its coefficient vanishes. It is

(* (8 BesselI[0, Sqrt[lamda]] C[2])/(lamda ((8 + lamda) Hypergeometric0F1Regularized[3, lamda/4] + 
   lamda Hypergeometric0F1Regularized[4, lamda/4])) *)

Replacing lamda by - z^2 allows further simplification

FullSimplify[coef /. lamda -> -z^2]
(* (2 BesselJ[0, z] C[2])/((-8 + z^2) BesselJ[2, z] + 2 z BesselJ[3, z]) *)

Clearly, the numerator vanishes at BesselJZero[0, k] for all nonnegative integer k, suggesting that

lamda -> - BesselJZero[0, k]^2

Unfortunately, this substitution also causes the denominator to vanish, as can be shown by evaluating it numerically for various k. (FullSimplify cannot do the required simplification.) The expression for lamda thus is plausible but not necessarily correct.

Instead, solve the original system with this expression for lamda assumed from the outset. (We choose a particular value of k at random, because FullSimplify can make no headway otherwise.)

ss = DSolveValue[{2 x y''''[x] + 4 y'''[x] == lamda y''[x], y[0] == 0, y'[1/2] == 0, 
    y''[1/2] == lamda y[1/2]} /. lamda -> -BesselJZero[0, 17]^2, y[x], x] // FullSimplify
(* (2 x C[1] Hypergeometric0F1Regularized[2, -(1/2) x BesselJZero[0, 17]^2])/ 
   BesselJZero[0, 17]^2 *)

along with the warning message,

Solve::ztest1: Unable to decide whether numeric quantity 2 Hypergeometric0F1Regularized[3,-(1/4) BesselJZero[0,17]^2]-1/4 BesselJZero[0,17]^2 Hypergeometric0F1Regularized[3,-(1/4) BesselJZero[0,17]^2]-1/4 BesselJZero[0,17]^2 Hypergeometric0F1Regularized[4,-(1/4) BesselJZero[0,17]^2] is equal to zero. Assuming it is.

In fact, it is equal to zero, as can be seen by evaluating it numerically. FullSimplify does not know an identity to demonstrate that this expression vanishes.

Finally, take the derivative of ss and evaluate it at x = 0.

FullSimplify[D[ss, x]]
% /. x -> 0
(* (2 C[1] Hypergeometric0F1Regularized[1, -(1/2) x BesselJZero[0, 17]^2])/
   BesselJZero[0, 17]^2 *)
(* -((2 C[1])/BesselJZero[0, 17]^2) *)

demonstrating that the earlier expression for lamda is correct.

Answer 2

First solving

Clear[x, y]
sol = DSolve[2 x y''''[x] + 4 y'''[x] == lambda y''[x], y, x][[1]]
yx = y[x] /. sol

and then the boundary conditions

equ1 = (D[yx, x] /. {x -> 1/2}) // FullSimplify
equ2 = ((D[yx, {x, 2}] - lambda yx) /. {x -> 1/2}) // FullSimplify
equ3 = Limit[yx, x -> 0] // FullSimplify
equ4 = Limit[D[yx, x], x -> 0] - dy0

giving

(* (2 (-1 + BesselI[0, Sqrt[lambda]]) C[1] - 4 BesselK[0, Sqrt[lambda]] C[2])/lambda + C[4]*)
(*C[1] - 1/2 lambda (2 C[3] + C[4])*)
(*(4 C[2])/lambda^2 + C[3]*)
(*-dy0 + DirectedInfinity[-Sign[C[2]]]/lambda*)*)

Analyzing those results we conclude that to be feasible, the conditions should dictate C[2] = C[3] = dy0 = 0 so we follow

equ10 = equ1 /. {C[2] -> 0, C[3] -> 0};
equ20 = equ2 /. {C[2] -> 0, C[3] -> 0};
M = Grad[{equ10, equ20}, {C[1], C[4]}];
equ0 = Det[M] // FullSimplify
(* -BesselI[0, Sqrt[lambda]] *)

and then the eigenvalues are the real roots for BesselI[0, Sqrt[lambda]] or

Plot[-BesselI[0, Sqrt[lambda]], {lambda, -150, 0}]

or

Solve[{BesselI[0, Sqrt[lambda]] == 0, lambda < 0}, lambda][[1]]

(*lambda -> ConditionalExpression[-BesselJZero[0, C[1]]^2, C[1] \[Element] Integers && C[1] >= 1]*)