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?

@HenrikSchumacher It is not obvious to me that NDEigensystem can handle this problem, due to the appearance of the eigenvalue in one boundary condition and the finiteness of the solution as another boundary condition. — bbgodfrey, Dec 06 '19 at 02:18
I have a package (CompoundMatrixMethod) for finding eigenvalues numerically that can cope with the eigenvalue in the boundary condition(e.g. https://mathematica.stackexchange.com/questions/197165/eigenvalue-dependent-boundary-conditions-mathematica/197257), but not the finiteness constraint. I can find eigenvalues when the fourth BC is homogeneous (e.g. $y'(0)=0$ or $y''(0)=0$ or $y'''(1/2)=0$, but not when you only wish to impose that it is finite. — SPPearce, Dec 06 '19 at 06:08

score 7 · Accepted Answer · answered Dec 06 '19 at 02:13

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.

Thank you! Indeed, the expression for 'ss' vanishes identically to zero. — user31899, Dec 06 '19 at 02:26

Cesareo · Answer 2 · 2019-12-06T11:05:27.300

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]*)

Eigenvalues of a fourth-order ODE

2 Answers2