Eigen values of a third order linear homogenous ODE

Question

From a system of PDEs where i used the following ansatz: $$\theta_w(x,y) = e^{-\beta_h x} f(x) e^{-\beta_c y} g(y)$$. $F(x) := \int f(x) \, \mathrm{d}x$ and $G(y) := \int g(y) \, \mathrm{d}y$ So, $$\theta_w(x,y) = e^{-\beta_h x} F'(x) e^{-\beta_c y} G'(y)$$ I have two linear ODEs of which one i am mentioning here $$ \lambda_hF'''-2\lambda_h\beta_hF''+((\lambda_h\beta_h-1)\beta_h-\mu)F'+\beta_h^2F=0 $$ with the b.c. $$ F(0)=0 $$ $$ \frac{F''(0)}{F'(0)}=\beta_h $$

$$ \frac{F''(1)}{F'(1)}=\beta_h$$

Here $\mu$ is the separation constant

A general solution of this ODE will be of the form:

$$ F(x) = \sum_k C_k e^{-\delta_k(\mu)x}$$ $k$ attains values $1,2,3$ because of the cubic nature of the characterstic equation.

I need to find the eigenvalues of this ODE. I know i need to consider $\mu><=$$0$. A second order ODE has a well defined form of its roots which is not the case here.

I will mention here that i am not too well accustomed to Mathematica But is there any way through which these eigen values could be determined using Mathematica ? I am at a dead end here.

(Finally i need to repeat the complete procedure for $G$ and then proceed to find $\theta_w$)

Answer 1

I have a package for numerically calculating solutions of eigenvalue problems using the Evans function via the method of compound matrices, which is hosted on github. See my answers to other questions or the github for some more details.

First we install the package (only need to do this the first time):

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

Then we first need to turn the ODEs into a matrix form $\mathbf{y}'=\mathbf{A} \cdot \mathbf{y}$, using my function ToMatrixSystem:

Needs["CompoundMatrixMethod`"]
sys = ToMatrixSystem[{λ F'''[x] - 2 λ β F''[x] + ((λ β - 1) β - μ) F'[x] + β^2 F[x] == 0}, 
         {F[0] == 0, F''[0] == β F'[0], F''[1] == β F'[1]}, F, {x, 0, 1}, μ]

The object sys contains the matrix $\mathbf{A}$, as well as similar matrices for the boundary conditions and the range of integration.

Now the function Evans will calculate the Evans function (also known as the Miss-Distance function) for any given value of $\mu$ provided that $\lambda$ and $\beta$ are specified; this is an analytic function whose roots coincide with eigenvalues of the original equation.

Evans[1, sys /. {λ -> 1, β -> 2}]
(* -0.023299 *)

We can also plot this:

Plot[Evans[μ, sys /. {λ -> 1, β -> 1}], {μ, -300, 300}, AxesLabel -> {"μ", "E(μ)"}]

For these parameter values you can see there are six eigenvalues within the plot region, the function findAllRoots (by user Jens, available from this post) will give you them all:

findAllRoots[Evans[μ, sys /. {λ -> 1, β -> 1}], {μ, -300, 300}]
(* {-247.736, -158.907, -89.8156, -40.4541, -10.7982, -0.635232} *)

For the values $\lambda=0.02, \beta=10$, it helps to remove the normalisation that is applied by default in Evans, and I've also plotted the values given by bbgodfrey's solution:

Plot[Evans[μ, sys /. {λ -> 1/50, β -> 10}, NormalizationConstants -> 0], 
    {μ, -21, -2}, Epilog -> 
  Point[{#, 0} & /@ {-20.88, -17.48, -14.36, -11.53, -9.03, -6.92, -5.26, -4.08, -3.16}]]

It looks like there is an infinite set of solutions for negative $\mu$ for this set of parameters. For $\beta=0$ the eigenvalues are $-\lambda (n \pi)^2$, the general behaviour as $\mu->-\infty$ looks the same.

As you vary $\beta$ and $\lambda$ you can see how the function changes and eigenvalues can move and coalesce.

Note there is a discrepancy with bbgodfrey's results, which missed the root at $-3.392$ and gave a spurious one at $-3.165$. That isn't to say that you can't manage to get spurious solutions by applying FindRoot to this technique if you aren't careful.

Additionally, as the Evans function is (complex) analytic, you can use the Winding Number to find how many roots are within a contour. I have some functions for this, so you can see that there are only these 9 roots within the circular contour of size 21 located at the origin:

PlotEvansCircle[sys /. {λ -> 1/50, β -> 10}, ContourRadius -> 21, 
   nPoints -> 500, Joined -> True]

Here the left hand side is a circle in the complex plane, and the right hand side is the Evans function applied to the points on that circle. The winding number is how many times the right hand contour goes round the origin. There is also a PlotEvansSemiCircle for seeing if there is a root in the positive half-plane (as that is a common question for instability).

To get the eigenfunctions, replace the right BC with $F'(0)$ to some arbitrary value and use NDSolve. This is at high precision to check the boundary conditions are being met:

Block[{β = 10, λ = 1/50},
 μc = μ /. FindRoot[Evans[μ, sys, NormalizationConstants -> 0, 
      WorkingPrecision -> 50], {μ, -3}, WorkingPrecision -> 50, 
     AccuracyGoal -> 15, PrecisionGoal -> 15] // Quiet;
 sol = First[NDSolve[{λ F'''[x] - 2 λ β F''[x] + ((λ β - 1) β - μ) F'[x] + β^2 F[x] == 0, 
   F[0] == 0, F''[0] == β F'[0], F'[0] == 1/10} /. μ -> μc, F, 
   {x, 0, 1}, WorkingPrecision -> 50, AccuracyGoal -> 30, PrecisionGoal -> 30]];
 Plot[(F[x]/F[1]) /. sol, {x, 0, 1}, 
  PlotLabel -> 
   "Eigenfunction for μ=" <> ToString@Round[μc, 0.0001] <> ", with BC error: " <> 
    ToString@Round[((F''[1] - β F'[1]) /. sol), 0.0001], 
  AxesLabel -> {"x", "F(x) (scaled)"}, PlotRange -> All]
 ]

The eigenfunctions get increasingly oscillatory as the eigenvalue gets more negative:

roots = Reverse@findAllRoots[Evans[μ, sys /. {λ -> 1/50, β -> 10}, 
    NormalizationConstants -> 0], {μ, -200, -2}]
Plot[Evaluate[Table[(F[x]/F[1]) /. 
    NDSolve[{λ F'''[x] - 2 λ β F''[x] + ((λ β - 1) β - μ) F'[x] + β^2 F[x] == 0, 
 F[0] == 0, F''[0] == β F'[0], F'[0] == 1/1000} /. {λ -> 1/50, β -> 10}, 
     F, {x, 0, 1}], {μ, roots}]], {x, 0, 1}, PlotRange -> All, 
 AxesLabel -> {"x", "F(x), scaled"}]

The method will work for higher order systems (up to 10th order generally), and does not require DSolve to be able to solve the underlying ODE.

Answer 2

Error corrected by deleting Numerator from definition of s.

This problem can be transformed symbolically to an algebraic equation by the method described in my earlier comment.

eq = ll f'''[x] - 2 ll b f''[x] + ((ll b - 1) b - m) f'[x] + b^2 f[x] == 0;
s = Simplify[DSolveValue[{eq, f''[0] == b f'[0], f''[1] == b f'[1]}, f[x], x]
    /. x -> 0 /. C[1] -> 1] 

The resulting expression in terms of Root is too lengthy to be reproduced here and must, in any case, be solved numerically to provide values for m. For the parameters given by the OP in a comment,

Table[FindRoot[s /. {ll -> 1/50, b -> 10}, {m, m0}] // Chop // Values, 
    {m0, -21, -1, 1/10}] // Flatten;
Union[%, SameTest -> (Abs[#1 - #2] < 10^-4 &)]

(* {-20.8885, -17.4846, -14.3644, -11.5367, -9.03846, -6.92843, 
    -5.26567, -4.08628, -3.39259} *)

There appear to be no larger values of m for this set of parameters. There are, of course, many more-negative values, perhaps an infinite number.

Here are a few more cases. For {ll -> 1/10, b -> 10},

(* {-13.9638, -7.91839, -3.94417, -1.94864} *)

and for {ll -> 1/1000, b -> 10},

(* {-20.5949, -19.94, -19.3031, -18.684, -18.0821, -17.4972, -16.9288, 
    -16.3763, -15.8388, -15.3158, -14.806, -14.3085, -13.822, -13.3451, 
    -12.8764, -12.4145, -11.9582, -11.5069, -11.0604, -10.6193, -10.1852, 
    -9.76044, -9.34825, -8.95222, -8.57622, -8.22411, -7.89954, -7.6058, 
    -7.34579, -7.12193, -6.93623, -6.79024, -6.68516, -6.62179} *)