Hello everybody in Mathematica SE. Although my question is related to flow stability analysis, this should be a general application of MMA to solve a system of ODEs. Thank you for your suggestion!

For the special case of zero streamwise wave number ($\alpha=0$), the dispersion relation (frequency $\omega$ as a function of spanwise wave number $\beta$ and other control parameter) and the eigenfunctions could be determined analytically because the Orr-Sommerfeld and Squire equations (OSS eqn, see Fig.1) reduce to constant coefficient ODEs.

Fig.1 snapshot of OSS eqns (from a highly cited textbook in the flow stability field): $\omega$ is a frequency, $\alpha$ and $\beta$ are streamwise and spanwise wave number, respectively, and $k^2=\alpha^2+\beta^2$, $U(y)$ is a mean flow, for plan Pliseuille flow $U=1-y^2$, $D=\frac{d}{d y}$, and a prime stands for differential w.r.t $y$, $\tilde{v}(y)$ and $\tilde{\eta}(y)$ are normal velocity and normal vorticity, and $Re$ is Reynolds number (just a constant parameter). enter image description here

My question arises from the following snapshot from this book when I have been trying to find the solution to the OSS equation and the dispersion relation in the special case with Mathematica.

Fig.2 enter image description here

My traces:

  1. For plan Poiseuille flow: $U(y)=1-y^2$; with $\alpha=0$, it follows that $k=\beta$, then the OSS eqn reduces to

$$\left[-\mathrm{i} \omega (D^2-\beta^2)-\frac{1}{Re}(D^2-\beta^2)^2 \right] \tilde{v} =0, \tag{1}$$

$$\left[-\mathrm{i} \omega -\frac{1}{Re}(D^2-\beta^2) \right] \tilde{\eta} = -\mathrm{i} \beta U^{\prime} \tilde{v}. \tag{2}$$

  1. With the following code, DSolve returns the equations unsolved.

    DSolve[{-I \[Omega]*(v''[y] - \[Beta]^2 v[y]) - 
    1/Rey (v''''[y] - 2 \[Beta]^2 v''[y] + \[Beta]^4 v[y]) == 0,
    -I*\[Omega] \[Eta][y] - 
    1/Rey*(\[Eta]''[y] - \[Beta]^2*\[Eta][y]) == -I \[Beta] (-2 y) v[y],
    v[1] == 0, v[-1] == 0, v'[1] == 0, 
    v'[-1] == 0, \[Eta][-1] == 0, \[Eta][1] == 0}, {v[y], \[Eta][y]}, y]
  2. If I solve the two equations separately

    DSolve[{-I \[Omega]*(v''[y] - \[Beta]^2 v[y]) - 
    1/Rey (v''''[y] - 2 \[Beta]^2 v''[y] + \[Beta]^4 v[y]) == 0,
    v[1] == 0, v[-1] == 0, v'[1] == 0, v'[-1] == 0}, v[y], y]
    (*{{v[y] -> 0}}*)

With $v(y)=0$, the second equation becomes homogeneous. The code gives {{\[Eta][y] -> 0}}:

    DSolve[{-I*\[Omega] \[Eta][y] - 
    1/Rey*(\[Eta]''[y] - \[Beta]^2*\[Eta][y]) == 0, \[Eta][-1] == 
    0, \[Eta][1] == 0}, \[Eta][y], y]

My questions are:

(1) How to modify the code to solve the OSS equation in this special case analytically?

(2) How to introduce an intermediate variable $\mu$ in the solution of $\omega$ (eigenvalue), see eqs.(3.61)--(3.63) in Fig.2, which is related to $\beta$ through a transcendental equation. Similarly, how to introduce another intermediate variabel $n$ in (3.64), which is the solution of Eq.(2). But, I found that $n$ was not defined there in that book.

(3) How could I account for odd and even modes in MMA's DSolve?

Thank you in advance!

  • $\begingroup$ DSolve should be able to give the answer, which is v = η = 0, but it cannot without some human assistance. The appropriate function, I would expect, is DEigensystem, but it returns unevaluated, which is bizarre, because the problem is so simple. So, I suggest you use DSolve to solve the first equation without boundary conditions, after which you can do the same for the second equation, Then, apply the boundary conditions, yielding six linear algebraic equations for the six constants of integration. Construct the determinant of the equations to obtain a dispersion relation. $\endgroup$
    – bbgodfrey
    Sep 8, 2018 at 14:19
  • $\begingroup$ With respect to your second question, the eigenfunctions never appear in the equations for the eigenvalues. $\endgroup$
    – bbgodfrey
    Sep 8, 2018 at 14:22
  • $\begingroup$ Thanks, @bbgodfrey I updated my post according to your first comment. I guess my question 2 was stupid, but could u pls give some simple explanation? Thanks a lot! $\endgroup$
    – lxy
    Sep 8, 2018 at 14:40
  • $\begingroup$ It would be nice to solve these equations numerically for the Poiseuille flow. Investigate the dependence of solutions on parameters. Then construct a non-linear extension of the model. Mathematics can solve all these problems. $\endgroup$ Sep 8, 2018 at 15:49
  • $\begingroup$ @bbgodfrey, although your DEigensystem was unevaluated, could you pls post it that all of us can debug together? However, it depends on you of course :) $\endgroup$
    – lxy
    Sep 9, 2018 at 5:20

2 Answers 2


Dispersion relations for this question can be derived in a straightforward manner, as I described in an earlier comment above. Solve the two ODEs in turn but without boundary conditions.

s = DSolveValue[{-I ω*(v''[y] - β^2 v[y]) - 1/Rey (v''''[y] - 2 β^2 v''[y] 
    + β^4 v[y]) == 0} /. ω -> -I (β^2 - μ^2)/Rey, v[y], y]
(* E^(-y μ) C[1] + E^(y μ) C[2] + E^(-y β) C[3] + E^(y β) C[4] *)
s1 = DSolveValue[(-I*ω η[y] - 1/Rey*(η''[y] - β^2*η[y]) == 
    -I β (-2 y) s) /. ω -> -I (β^2 - μ^2)/Rey, η[y], y] // Simplify
(* (1/(4 (β - μ)^2 μ^3 (β + μ)^2))E^(-y (β + 3 μ)) (8 I E^(3 y μ)
   Rey β μ^3 (-2 β - y β^2 + y μ^2) C[3] - 8 I E^(2 y β + 3 y μ)
   Rey β μ^3 (-2 β + y β^2 - y μ^2) C[4] + E^(y (β + 4 μ)) (β^2 - μ^2)^2 
   (-I Rey β (1 - 2 y μ + 2 y^2 μ^2) C[2] + 4 μ^3 C[5]) + E^(y (β + 2 μ)) 
   (β^2 - μ^2)^2 (I Rey β (1 + 2 y μ + 2 y^2 μ^2) C[1] + 4 μ^3 C[6])) *)

where ω has been replaced by -I (β^2 - μ^2)/Rey, consistent with the substitution shown in the question. Now, applying the six boundary conditions yields six linear, homogeneous equations for the six C[_], which must all equal zero, unless the determinant of the coefficient matrix vanishes. The determinant, therefore, is the dispersion relation.

{s /. y -> 1, s /. y -> -1, D[s, y] /. y -> 1, D[s, y] /. y -> -1, 
    s1 /. y -> -1, s1 /. y -> 1};
Det[CoefficientArrays[%, Array[C, 6]] // Last] // Factor;
Factor@Simplify[-% // ExpToTrig]/64
(* Cosh[μ] Sinh[μ] (μ Cosh[μ] Sinh[β] - β Cosh[β] Sinh[μ]) 
   (β Cosh[μ] Sinh[β] - μ Cosh[β] Sinh[μ]) *)

The third and fourth factors, which also can be obtained by solving the first ODE only, can be rewritten as

Thread[Simplify[-%[[3]]/(Cosh[μ] Cosh[β]), %[[4]]/(Sinh[μ] Sinh[β])}] == 0]
(* -μ Tanh[β] + β Tanh[μ] == 0, -μ Coth[β] + β Coth[μ] == 0} *)

which are similar to, but not identical with, the odd and even mode dispersion results cited in the question. Note that the only solutions of interest to these two equations are pure imaginary.

The first two factors of the determinant arise from the second ODE.

Thread[List @@ %%[[1 ;; 2]] == 0]
(* {Cosh[μ] == 0, Sinh[μ] == 0} *)

Here, too, roots are pure imaginary, given by I (n-1/2) Pi and I n Pi, respectively. The first differs from the corresponding expression in the question by a factor of I; the second is not given in the question.

Even though these solutions differ from those cited in the question, I am inclined to believe them (except, perhaps, for Sinh[μ] == 0, which may be spurious), because the computation is so straightforward.

  • $\begingroup$ Thanks,@bbgodfrey I am inclined to believe there are typos in the cited book. $\endgroup$
    – lxy
    Sep 10, 2018 at 3:17

(1)&&(3) It is possible to construct analytical solutions for even and odd modes, using the boundary conditions at the center of the channel and on one of the walls. These solutions depend on an arbitrary amplitude (A, B). Using these solutions on the second wall, we find the dispersion equation

solodd = DSolve[{-I \[Omega]*(v''[y] - \[Beta]^2 v[y]) - 
     1/Rey (v''''[y] - 2 \[Beta]^2 v''[y] + \[Beta]^4 v[y]) == 0, 
   v[1] == 0, v[0] == 0, v'[1] == 0, v'[0] == A}, v[y], y]
soleven = 
 DSolve[{-I \[Omega]*(v''[y] - \[Beta]^2 v[y]) - 
     1/Rey (v''''[y] - 2 \[Beta]^2 v''[y] + \[Beta]^4 v[y]) == 0, 
   v[1] == 0, v[0] == B, v'[1] == 0, v'[0] == 0}, v[y], y]

(2) The dependence of the dispersion equation on the amplitude appears in the second approximation, whereas the Orr-Sommerfeld equation is the first (linear) approximation. The complete solution is a linear combination of even and odd solutions, we have

V=First[v[y] /. soleven]+First[v[y] /. solodd]; V1=D[V,y];

We calculate the values of these functions on the wall

eq = {V, V1} /. y -> -1

By the conditions of the problem eq==0, from this we find a linear system of equations for A,B, the determinant of which must be equal to zero. This is the dispersion equation. Using V, we can try to solve equation

sol = DSolve[{-I*\[Omega] \[Eta][y] - 
         y] - \[Beta]^2*\[Eta][
          y]) == -I \[Beta] (-2 y) V, \[Eta][-1] == 0, \[Eta][1] == 
    0}, \[Eta][y], y]

But the solution turns out to be so cumbersome that it can not be surveyed.

I derived the dispersion equation in two ways and got the same result.The first method is shown above. The second way is shown below.It is possible to simplify the derivation of the dispersion equation. To do this, we find a solution that satisfies the boundary conditions on one wall:

sol1 = DSolve[{-I \[Omega]*(v''[y] - \[Beta]^2 v[y]) - 
     1/Rey (v''''[y] - 2 \[Beta]^2 v''[y] + \[Beta]^4 v[y]) == 0, 
   v[1] == 0, v'[1] == 0}, v[y], y];

V = ExpToTrig[First[v[y] /. sol1]];

Then on the other wall we have a system of two equations:

    eq = {V, D[V, y]} /. y -> -1
m = CoefficientArrays[eq, Array[C, 2]] // Last

FullSimplify[Det[m] /. {\[Omega] -> -I (\[Beta]^2 - \[Mu]^2)/Rey}]

The answer looks like neither in the book, nor at the fact that bbgodfrey

2 \[Beta] \[Mu] - 
  2 \[Beta] \[Mu] Cosh[2 \[Beta]] Cosh[
    2 \[Mu]] + (\[Beta]^2 + \[Mu]^2) Sinh[2 \[Beta]] Sinh[2 \[Mu]] = 0

However, if we compare the solutions of this equation and the system of equations {-\[Mu] Tanh[\[Beta]] + \[Beta] Tanh[\[Mu]] == 0, -\[Mu] Coth[\[Beta]] + \[Beta] Coth[\[Mu]] == 0} using ContourPlot[], then we see that their solutions are identical $\beta =\pm \mu$. This is a trivial solution $\omega =0$. We can indicate other non-trivial solutions, for example

sol = DSolve[{-I \[Omega]*(v''[y] - \[Beta]^2 v[y]) - 
     1/Rey (v''''[y] - 2 \[Beta]^2 v''[y] + \[Beta]^4 v[y]) == 0, 
   v[1] == 0, v[-1] == 0}, v[y], y];

eq = {(D[V, y] /. y -> -1), (D[V, y] /. y -> 1)};
m = CoefficientArrays[eq, Array[C, 2]] // Last;

V = ExpToTrig[First[v[y] /. sol]]

Here the dispersion equation has the form (with Rey \[Omega] -> \[CapitalOmega])

 deq= -2 \[Beta] Sqrt[\[Beta]^2 - 
       I \[CapitalOmega]] + \[Beta] Sqrt[\[Beta]^2 - I \[CapitalOmega]]
       Cosh[2 \[Beta] - 
        2 Sqrt[\[Beta]^2 - I \[CapitalOmega]]] + \[Beta] Sqrt[\[Beta]^2 - 
       I \[CapitalOmega]]
       Cosh[2 (\[Beta] + 
          Sqrt[\[Beta]^2 - I \[CapitalOmega]])] + (-2 \[Beta]^2 + 
        I \[CapitalOmega]) Sinh[2 \[Beta]] Sinh[2 Sqrt[\[Beta]^2 - I \[CapitalOmega]]]
 Re[deq] == 0, { \[Beta], -10, 10}, {\[CapitalOmega], -10, 10}]

Using ContourPlot, we find the roots of equation


  • $\begingroup$ Thanks a lot, @Alex Trounev, the most difficult part, when v[y] being obtained, could be the Squire equation since it is nonhomogeneous. Actually, \[Eta][y] is what I want. $\endgroup$
    – lxy
    Sep 9, 2018 at 14:11
  • $\begingroup$ Could you expand your answer to show the derivation of the dispersion relation? Thanks. $\endgroup$
    – bbgodfrey
    Sep 9, 2018 at 14:49
  • $\begingroup$ @ bbgodfrey I added an explanation of how to get the dispersion equation. $\endgroup$ Sep 9, 2018 at 18:08
  • $\begingroup$ I was hoping the see the dispersion relation you derived. I have done a similar calculation but obtained a different dispersion relation from that in the question. $\endgroup$
    – bbgodfrey
    Sep 9, 2018 at 18:09
  • $\begingroup$ @bbgodfrey I published the dispersion equation, it is quite identical to yours. $\endgroup$ Sep 10, 2018 at 6:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.