Eigenvalue problem with NDSolve

Question

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.

Answer 1

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

Answer 2

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"]

Load the package:

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;

And your equations:

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

Answer 3

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.