Eigenvalue problem with NDSolve

I am trying to solve the following system of linear ODEs. It is an eigenvalue problem.

\[Chi] = 1;
m = 1.495;
\[Theta] = -\[Chi]^(-1/m)*(-1 + \[Chi]^(1/m));
n = Log[Abs[\[Chi]]];

l = 2;
k = 0;
a = Sqrt[(l^2 + k^2)];

Pr = 0.1;
Ta = 10^5;
\[Phi] = Pi/4;
Ra = 2*10^(05);

sol1 = NDSolve[{(a^2 Pr + \[Sigma]) Subscript[Z, 0][z] + (
m Pr \[Theta] Subscript[Z, 1][z])/(1 + z \[Theta]) ==
Pr (Sqrt[
Ta] ((I l Cos[\[Phi]] + (m \[Theta] Sin[\[Phi]])/(
1 + z \[Theta])) Subscript[W, 0][z] +
Sin[\[Phi]] Subscript[W, 1][z]) +
Derivative[1][Subscript[Z, 1]][z]),
Pr (((3 (-2 + m) m \[Theta]^4 + (a + a z \[Theta])^4) Subscript[W,
0][z])/(1 + z \[Theta])^4 + (
m \[Theta] ((-4 + m) \[Theta] Subscript[W, 2][z] +
2 (1 + z \[Theta]) Subscript[W, 3][z]))/(1 + z \[Theta])^2 +
Derivative[1][Subscript[W, 3]][z]) ==
a^2 Pr Ra Subscript[S, 0][z] + (
2 a^2 m^2 Pr \[Theta]^2 Subscript[W, 0][z])/(
3 (1 + z \[Theta])^2) + (
3 (-2 + m) m Pr \[Theta]^3 Subscript[W, 1][z])/(1 +
z \[Theta])^3 + (2 a^2 m Pr \[Theta] Subscript[W, 1][z])/(
1 + z \[Theta]) +
2 a^2 Pr Subscript[W, 2][
z] + \[Sigma] ((-a^2 - (
m \[Theta]^2)/(1 + z \[Theta])^2) Subscript[W, 0][z] + (
m \[Theta] Subscript[W, 1][z])/(1 + z \[Theta]) +
Subscript[W, 2][z]) +
Pr Sqrt[Ta] (I Cos[\[Phi]] ((k m \[Theta] Subscript[W, 0][z])/(
1 + z \[Theta]) + l Subscript[Z, 0][z]) +
Sin[\[Phi]] Subscript[Z, 1][z]), (a^2 (1 +
z \[Theta])^-m + \[Sigma]) Subscript[S, 0][
z] == (1 + z \[Theta])^(-1 -
m) (\[Theta] Subscript[S, 1][z] + (1 + z \[Theta])^
m Subscript[W, 0][z] + (1 + z \[Theta]) Derivative[1][
Subscript[S, 1]][z]),
Subscript[W, 1][z] == Derivative[1][Subscript[W, 0]][z],
Subscript[W, 2][z] == Derivative[1][Subscript[W, 1]][z],
Subscript[W, 3][z] == Derivative[1][Subscript[W, 2]][z],
Subscript[S, 1][z] == Derivative[1][Subscript[S, 0]][z],
Subscript[Z, 1][z] == Derivative[1][Subscript[Z, 0]][z],
Subscript[S, 0][0] == 0, Subscript[S, 0][1] == 0,
Subscript[W, 0][0] == 0, Subscript[W, 0][1] == 0,
Subscript[Z, 1][0] == 0, Subscript[Z, 1][1] == 0,
Subscript[W, 2][0] +
m \[Theta] Subscript[W, 1][0]/(1 + 0 \[Theta]) == 0,
Subscript[W, 2][1] +
m \[Theta] Subscript[W, 1][1]/(1 + \[Theta]) == 0}, {Subscript[W,
0], Subscript[Z, 0], Subscript[S, 0], \[Sigma]}, {z, 0, 1}]


Additional condition that can be used is

Subscript[Z, 0][0]==1


I want to solve for the unknown functions and sigma. But right now I am unable to do it. I tried to follow the approach given in https://reference.wolfram.com/language/tutorial/NDSolveBVP.html but it didn't work for me. Please help.

Thank you.

• Have you tried NDEigensystem? Jan 4, 2022 at 11:50
• @user21 I have tried it but since sigma in the above equation is unknown, I am not sure how to place the equations in the NDEigensystem function. Jan 4, 2022 at 12:18
• Your code contains Subscript[Z, 1][0] == 1, Subscript[Z, 1][0] == 0, one of which is incorrect. What are the correct boundary conditions? Also, you need a ninth boundary condition to solve this problem, because the eigenvalue also is an unknown. Jan 4, 2022 at 20:04
• @bbgodfrey thank you for pointing it out. Indeed, those conditions are incorrect. I have made changes to those boundary conditions. Jan 5, 2022 at 3:15
• Thanks for the additional condition. However, the code now contains Subscript[Z, 1][0] == 0, Subscript[Z, 1][0] == 0, which is a duplicate. Another boundary condition is needed to give a total of nine. Eight is sufficient only if all boundary conditions are homogeneous. Jan 5, 2022 at 5:10

Although this question probably can be solved with NDSolve, it is more convenient to solve it with DSolve. To begin, name the ODEs and boundary conditions given in the question as eqs with all quantities Rationalized. Then solve the ODEs with the boundary conditions at z = 0. (Including all boundary conditions at once produces ungainly results.)

v = Variables[Subtract @@@ eqs[[;; 8]]][[2 ;; 9]]
dsol = DSolve[Join[eqs[[;; 8]], eqs[[{9, 11, 13, 15}]]], v, z] // Flatten


With a LeafCount of 9320, dsol is too large to reproduce here. It consists of four equations involving RootSums and four constants of integration.

c = Cases[dsol, C[_], Infinity] // Union
(* {C[2], C[4], C[6], C[7]} *)


Now, apply the four boundary conditions at z = 1 and express the result as a matrix.

eqs[[10 ;; 16 ;; 2]] /. (dsol /. z -> 1);
CoefficientArrays[%, c] // Last;
disp = Det[%];


The roots of disp == 0 are the desired eigenvalues, the first four of which are

FindRoot[disp, {σ, -.1}, WorkingPrecision -> 30] // Chop
(* {σ -> -0.0873219119392251168866661775144} *)


and similarly

(* {σ -> -1.27987460528269767318780120707} *)
(* {σ -> -68.8450688277416266702085885341} *)
(* {σ -> -81.1196686184117486120744862461} *)


A plot of negative σ is given by

Show[
Plot[Chop[disp], {σ, -260, -67}, ImageSize -> Large, AxesLabel -> {σ, "disp"},
LabelStyle -> {15, Bold, Black}, PlotRange -> {{-260, 0}, Automatic},
Plot[Re[disp], {σ, -2, 1/4}, ImageSize -> Large] ]


• it is an interesting approach to solving this problem. Thank you so much. I am going to use it. Could you please also write how you entered the equations? Also, I am posting my solution obtained with NDSolve. Please comment on that. Jan 5, 2022 at 9:18
• @GAURAVMAURYA, In answer to your question about the equations I used, I used eqs = { ... }, where { ... } is the first set of curly brackets in NDSolve in the question and all the code it contains. I also set m = 1495/1000 and Pr = 1/10 in the constants to obtain a more accurate answer. Jan 5, 2022 at 14:19
• @bbgodgrey Thank you. I was able to do that by carefully going through your code. How can we approach the problem when both Ra and Sigma are unknown parameters? It is true that the code I presented fails at many points. Jan 6, 2022 at 5:30
• @GAURAVMAURYA Just as SPPearce noted in response to a similar question, my approach can provide σ as a function of Ra. To determine values for both of then requires an additional relation between the two quantities. Is one available? Jan 6, 2022 at 13:09
• The additional relation would be $a(Ra)=0$, where $\sigma=a+i b$. Can that be used? Jan 6, 2022 at 14:11

Here is an alternative to bbgodfrey's typically excellent analytical solution, that will work for equations that can't be solved by DSolve.

I happen to have written a package to solve eigenvalue equations of this type, using the Evans function.

The Evans function is an analytic function whose roots correspond to the eigenvalues. Some details are available at these two questions, my github page and this PDF, or search for CompoundMatrixMethod to see my previous answers here.

Install the package (also available on my github page):

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


Needs["CompoundMatrixMethod"]


Parameters (without setting Ra, and rationalising):

χ = 1;
m = Rationalize[1.495];
θ = -χ^(-1/m)*(-1 + χ^(1/m));
n = Log[Abs[χ]];

l = 2;
k = 0;
a = Sqrt[(l^2 + k^2)];

Pr = 1/10;
Ta = 10^5;
ϕ = Pi/4;


eqs = {(a^2 Pr + σ) Z0[z] + (m Pr θ Z1[z])/(1 + z θ) ==
Pr (Sqrt[Ta] ((I l Cos[ϕ] + (m θ Sin[ϕ])/(1 + z θ)) W0[z] + Sin[ϕ] W1[z]) + Z1'[z]),
Pr (((3 (-2 + m) m θ^4 + (a + a z θ)^4) W0[z])/(1 + z θ)^4 + (m θ ((-4 + m) θ W2[z] +  2 (1 + z θ) W3[z]))/(1 + z θ)^2 + W3'[z]) == a^2 Pr Ra S0[z] + (2 a^2 m^2 Pr θ^2 W0[z])/(3 (1 + z θ)^2) + (3 (-2 + m) m Pr θ^3 W1[z])/(1 + z θ)^3 + (2 a^2 m Pr θ W1[z])/(1 + z θ) +
2 a^2 Pr W2[z] + σ ((-a^2 - (m θ^2)/(1 + z θ)^2) W0[z] + (m θ W1[z])/(1 + z θ) + W2[z]) +
Pr Sqrt[Ta] (I Cos[ϕ] ((k m θ W0[z])/(1 + z θ) + l Z0[z]) +
Sin[ϕ] Z1[z]), (a^2 (1 + z θ)^-m + σ) S0[z] == (1 + z θ)^(-1 - m) (θ S1[z] + (1 + z θ)^m W0[z] + (1 + z θ) S1'[z]), W1[z] == W0'[z],
W2[z] == W1'[z], W3[z] == W2'[z], S1[z] == S0'[z], Z1[z] == Z0'[z],
S0[0] == 0, S0[1] == 0, W0[0] == 0, W0[1] == 0, Z1[0] == 0,
Z1[1] == 0, W2[0] + m θ W1[0]/(1 + 0 θ) == 0,
W2[1] + m θ W1[1]/(1 + θ) == 0};


We use the function ToMatrixSystem to convert this system into a matrix form

sys[Ra_] = ToMatrixSystem[eqs[[;; 8]],eqs[[9 ;; 16]], {S0, S1, W0, W1, W2, W3, Z0, Z1}, {z, 0, 1}, σ];


Now we can evaluate the Evans function for a given value of $$\sigma$$ , say $$\sigma = 1$$ and $$Ra = 200000$$.

Evans[1, sys[200000]]
(* 0.000545625 + 0. I *)


This is not zero, so $$\sigma=1$$ is not an eigenvalue of the system. So we can use FindRoot:

FindRoot[Evans[σ, sys[200000]], {σ, 1}]
{σ -> -1.27987 - 3.69498*10^-17 I}


Agreeing with bbgodfrey's solution.

If we plot the Evans function, we can see that it has further negative roots, again copying the analytical solution:

Plot[Evans[σ, sys[200000]], {σ, -250, 0}]


To find a value of $$Ra$$ where an eigenvalue crosses the y-axis is a bit tricky, as there are multiple eigenvalues, so you'll need to be very careful. However, one method is to ensure you start the FindRoot close to the previous root. I have some examples of this in the notebook on my github.

Some quick playing around suggests there is an eigenvalue that crosses between $$8000<=Ra=10000$$:

FindRoot[Evans[σ, sys[8000]], {σ, 10 I}]
FindRoot[Evans[σ, sys[10000]], {σ, 10 I}]

(* {σ -> -0.324265 + 14.6518 I} *)

(* {σ -> 0.523407 + 14.1909 I} *)


So we will try and see where this one crosses.

Make a temporary function:

Clear[rootfn]
rootfn[Ra_?NumericQ] := σ /.
FindRoot[Evans[σ, sys[Ra]], {σ, 15 I}]


And then find where the real part of this disappears:

FindRoot[Re[rootfn[Ra]], {Ra, 9000}]
(* {Ra -> 8758.34} *)

• Great approach (+1). Jan 5, 2022 at 14:30
• What is your extra condition/equation to include another parameter? The matrix system can depend on another parameter, which you can then vary to see how the eigenvalues change for instance. Jan 6, 2022 at 8:53
• Ok, you want to find a value of $Ra$ for which $\sigma$ becomes purely imaginary? Jan 7, 2022 at 9:02
• I have found one such value of $Ra$. But as there are a number of eigenvalues, you'll need to be careful to follow particular roots and there could easily be a smaller value that crosses. Please also check the equations, I have rewritten them without using Subscript. Also you don't need to rewrite them into a set of 1st order ODEs like Matlab requires (my package does that in ToMatrixSystem), you could just keep your actual structure. Jan 7, 2022 at 9:48
• @SPPearce This works nicely for me. Thank you so much for these solutions. Also, the Package works perfectly. Jan 7, 2022 at 16:20

I have used NDSolve to obtain the solution. Thanks to @bbgodfrey for pointing out the incorrect boundary conditions.

sol1 = NDSolve[{(a^2 Pr + I \[Sigma][z]) Subscript[Z, 0][z] + (
m Pr \[Theta] Subscript[Z, 1][z])/(1 + z \[Theta]) ==
Pr (Sqrt[
Ta] ((I l Cos[\[Phi]] + (m \[Theta] Sin[\[Phi]])/(
1 + z \[Theta])) Subscript[W, 0][z] +
Sin[\[Phi]] Subscript[W, 1][z]) +
Derivative[1][Subscript[Z, 1]][z]),
Pr (((3 (-2 + m) m \[Theta]^4 + (a + a z \[Theta])^4) Subscript[W,
0][z])/(1 + z \[Theta])^4 + (
m \[Theta] ((-4 + m) \[Theta] Subscript[W, 2][z] +
2 (1 + z \[Theta]) Subscript[W, 3][z]))/(1 + z \[Theta])^2 +
Derivative[1][Subscript[W, 3]][z]) ==
a^2 Pr Ra[z] Subscript[S, 0][z] + (
2 a^2 m^2 Pr \[Theta]^2 Subscript[W, 0][z])/(
3 (1 + z \[Theta])^2) + (
3 (-2 + m) m Pr \[Theta]^3 Subscript[W, 1][z])/(1 +
z \[Theta])^3 + (2 a^2 m Pr \[Theta] Subscript[W, 1][z])/(
1 + z \[Theta]) + 2 a^2 Pr Subscript[W, 2][z] +
I \[Sigma][
z] ((-a^2 - (m \[Theta]^2)/(1 + z \[Theta])^2) Subscript[W, 0][
z] + (m \[Theta] Subscript[W, 1][z])/(1 + z \[Theta]) +
Subscript[W, 2][z]) +
Pr Sqrt[Ta] (I Cos[\[Phi]] ((k m \[Theta] Subscript[W, 0][z])/(
1 + z \[Theta]) + l Subscript[Z, 0][z]) +
Sin[\[Phi]] Subscript[Z, 1][z]), (a^2 (1 + z \[Theta])^-m +
I \[Sigma][z]) Subscript[S, 0][z] == (1 + z \[Theta])^(-1 -
m) (\[Theta] Subscript[S, 1][z] + (1 + z \[Theta])^
m Subscript[W, 0][z] + (1 + z \[Theta]) Derivative[1][
Subscript[S, 1]][z]),
Subscript[W, 1][z] == Derivative[1][Subscript[W, 0]][z],
Subscript[W, 2][z] == Derivative[1][Subscript[W, 1]][z],
Subscript[W, 3][z] == Derivative[1][Subscript[W, 2]][z],
Subscript[S, 1][z] == Derivative[1][Subscript[S, 0]][z],
Subscript[Z, 1][z] ==
Derivative[1][Subscript[Z, 0]][z], \[Sigma]'[z] == 0, Ra'[z] == 0,
Subscript[S, 0][0] == 0, Subscript[S, 0][1] == 0,
Subscript[W, 0][0] == 0, Subscript[W, 0][1] == 0,
Subscript[Z, 1][0] == 0, Subscript[Z, 1][1] == 0,
Subscript[W, 2][0] +
m \[Theta] Subscript[W, 1][0]/(1 + 0 \[Theta]) == 0,
Subscript[W, 2][1] +
m \[Theta] Subscript[W, 1][1]/(1 + \[Theta]) == 0,
Subscript[Z, 0][0] == 1, Subscript[Z, 0][1] == 10}, {Subscript[W,
0], Subscript[Z, 0], Subscript[S, 0], \[Sigma], Ra}, {z, 0, 1},
Method -> {"Shooting",
"ImplicitSolver" -> {"Newton", "StepControl" -> "LineSearch"}}]


Here, you will see that I have kept Ra as an unknown parameter in addition to sigma. The eigenvalue and Ra can be extracted as

\[Sigma][0]/.sol1
Ra[0]/.sol1


The problem that I see here, is the extra conditions required to obtain a solution. These extra conditions might not be known beforehand. Other problem is if the extra conditions are homogeneous then you would not get the true eigenvalues. I have not checked the solution using another method.

• Your code fails with a FindRoot error, and the resulting answers are incorrect. You might try restoring Ra as a constant and deleting one of the extra boundary conditions. Jan 5, 2022 at 14:11
• Your code can be made to work by converting Ra back to a constant, removing Ra'[z] == 0 and Subscript[Z, 0][1] == 10, and returning I σ[z] to σ[z] (for consistency with the question). It also helps to add "StartingInitialConditions" -> σ[0] == -68845/1000 (or whatever initial guess you wish). This method, which is described in the documentation, has the advantage of giving the eigenfunction along with the eigenvalue, but it is slow and requires a reasonable initial guess. Jan 5, 2022 at 16:36
• Having read your discussion with SPPearce, I am confident that requiring Re[σ] = 0 will not determine Ra. (I presume you mean Re[σ] = 0, where σ is defined as in your answer just above, which would be the same as Im[σ] = 0, where σ is defined as in your question.) Jan 7, 2022 at 6:05
• Do you expect Ra` determined in this way to be real or complex? Jan 7, 2022 at 14:08
• There may be multiple values of Ra for which an eigenvalue becomes positive. I have found one such value below. An issue is finding where the eigenvalues that don't lie upon the negative real axis are, as the negative real ones seem to attract the FindRoot solution strongly. At Ra=200000 there is one at $\sigma = 8.049 + 11.36 I$. This moves around the complex plane as Ra is varied. Jan 7, 2022 at 16:44