# Eigen values of a third order linear homogenous ODE

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$$)

• Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! – Michael E2 Jan 14 '19 at 17:07
• People here generally like users to post code as Mathematica code instead of just images or TeX, so they can copy-paste it. It makes it convenient for them and more likely you will get someone to help you. You may find this this meta Q&A helpful – Michael E2 Jan 14 '19 at 17:07
• Also the boundary condition at 1 is identical on the top and the bottom of the fraction – KraZug Jan 14 '19 at 18:38
• I can solve it numerically for a given value of $\beta$ and $\lambda$ (assuming it should be $F''(1)/F'(1)=\beta$). Depending on their values, both positive and negative eigenvalues can exist. – KraZug Jan 14 '19 at 19:46
• Try solving the ODE and two of the three boundary conditions with DSolve. Then apply the third boundary condition to the resulting output to obtain an algebraic equation for the eigenvalues. – bbgodfrey Jan 14 '19 at 20:08

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.

• I should probably note that I don't currently have support for getting the eigenfunctions out of my package, which it looks like you'll probably need to proceed. – KraZug Jan 15 '19 at 10:13
• This was great and very useful. My only doubt was regarding the range of $\mu$ selected, which i guess i need to choose judiciously. Further since the general solution of $F$ is of the form $F(x) = \sum_k C_k e^{-\delta_k(\mu)x}$, the eigen functions should be of the form $e^{-\delta_k(\mu_n)t}$ The $\delta_k$ is the root of the characteristic equation – Indrasis Mitra Jan 15 '19 at 11:56
• @IndrasisMitra, added how to get the eigenfunctions and some examples. I would say not to accept an answer immediately, and wait to see if any other answers are forthcoming. – KraZug Jan 15 '19 at 13:28
• can the increasing oscillatory nature of eigenfunctions as the eigenvalues get more negative be inferred in a way that $F$ does not converge ? – Indrasis Mitra Jan 15 '19 at 18:20
• also i see that you have numerically solved the ODE with an arbitrary value of $F'(0)$ to find the eigenfuctions. I fail to grasp the logic behind this step – Indrasis Mitra Jan 15 '19 at 18:43

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

• Thanks a lot for this help. Now i have some idea. Some points i need to clarify. For the case $\lambda_h=0.1$ and $\beta_h=10$, you have posted a set of eigenvalues that come out. You also mention that there could be infinitely many more negative values. Any advice on how many of these eigen values should one consider to generate a solution ? I know the query is not related to mathematica here. Thanks again. – Indrasis Mitra Jan 15 '19 at 3:48
• @IndrasisMitra There is no good way to know without seeing the actual sum you are constructing. Is it a sum of eigenfunctions for the ODE in your question, or something else? In general, if the function to be approximated by a sum is reasonable uniform, then probably only a half-dozen or so eigenvalues are needed. – bbgodfrey Jan 15 '19 at 5:14
• I have edited my original post to show the type of series i expect. Yes, it is a sum of eigenfunctions from the ODE i posted. – Indrasis Mitra Jan 15 '19 at 5:29
• I'm getting a catastrophic precision loss for FindRoot at your largest root for {ll -> 1/50, b -> 10}: s /. {ll -> 1/50, b -> 10} /. FindRoot[s /. {ll -> 1/50, b -> 10}, {m, -3}, WorkingPrecision -> 20] gives me no significant digits to display (I don't believe it to be a root). – KraZug Jan 15 '19 at 9:33
• @KraZug I now see that it is not a root. My error was to use only the Numerator, which I plan to correct soon. Thanks. – bbgodfrey Jan 15 '19 at 18:05