Eigenvalues of a fourth-order ODE
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.
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]*)