# Eigen value solution of coupled ODEs

I want an eigen value solution of following coupled ODEs: But the code showing errors.

    gamma = 1.67; ep = 0.01; d = 0.05; H = 1.1; gpara = 0.38;
{svals, sfuns} =
NDEigensystem[{D[V[t, x], t] ==
Rho[t, x]/Exp[-x/H]*gpara -
1/gamma/Exp[-x/H]*(D[Rho[t, x], x] + D[T[t, x], x]) +
1.3*ep*D[V[t, x], {x, 2}],
D[Rho[t, x], t] == - D[Exp[-x/H]*V[t, x], x],
D[T[t, x], t] ==
D[Exp[-x/H]*V[t, x], x] - V[t, x]*D[Exp[-x/H], x] -
gamma*Exp[-x/H]*D[V[t, x], x] + gamma*d*D[T[t, x], {x, 2}]},
{V[t, x], Rho[t, x], T[t, x]}, t, {x, 0, 1}, 4];


Q1: I need to find numerical values of lambda (Eigen Values), Q2: In the above equations gpara is considered constant but it is a variable having 50 rows and 1 column in gpara.txt file, which is in another directory. I import gpara by Import["E\\ ......\\ gpara.txt"] to use in above equations, but its not working too. Please help to get the solution.

• Start with having a look at the reference page of (N)DEigensystem and the syntax used. – user21 Apr 26 at 6:04
• In your first eqn, it looks like you are using braces "{ }" around your derivatives instead of parentheses "( )". Also, you will need braces on the right hand side in your bc= statement. – LouisB Apr 26 at 6:07
• You have a term V[x]*Rho[x] that is going to be an impediment for NDEigensystem. Are you sure the term is correct? – user21 Apr 26 at 6:16
• @LouisB, thanks. Now i corrected as you suggested. – Danial Apr 26 at 6:20
• @user21, you are right, V[x]*Rho[x]*gpara has typing error, correct one is V[x]*P0', now corrected. Thanks – Danial Apr 26 at 6:52

It is necessary to bring the equations to the standard form and add Method->{"Eigensystem" -> "Direct"}

gamma = 1.67; ep = 0.01; d = 0.05; H = 1.1; gpara = 0.38;
{svals, sfuns} =
NDEigensystem[{-D[V[t, x], t] + Rho[t, x]/Exp[-x/H]*gpara -
1/gamma/Exp[-x/H]*D[Rho[t, x], x] -
1/gamma/Exp[-x/H]*D[T[t, x], x] + 1.3*ep*D[V[t, x], {x, 2}] == 0,
D[Rho[t, x], t] + Exp[-x/H]*D[V[t, x], x] - Exp[-x/H]*V[t, x]/H ==
0, -D[T[t, x], t] + Exp[-x/H]*D[V[t, x], x] -
gamma*Exp[-x/H]*D[V[t, x], x] + gamma*d*D[T[t, x], {x, 2}] == 0,
DirichletCondition[{V[t, x] == 0, Rho[t, x] == 0, T[t, x] == 0},
x == 0],
DirichletCondition[{V[t, x] == 0, Rho[t, x] == 0, T[t, x] == 0},
x == 1]}, {V, Rho, T}, t, {x, 0, 1}, 4,
Method -> {"Eigensystem" -> "Direct"}];
svals

(*Out[]= {0.00592203 + 0. I, -0.500739 + 0. I, -0.813679 +
0. I, -1.88774 + 0. I}*)
Table[Plot[ReIm[sfuns[[i, j]][x]], {x, 0, 1}, PlotRange -> All], {j,
1, 3}, {i, 1, 4}]


• Those solutions look very oscillatory except for Rho, do they solve the original PDE when you plug them back in? – KraZug Apr 27 at 19:50
• @KraZug it looks like you don't trust Mathematica 12 :) – Alex Trounev Apr 27 at 20:28
• When it shows these kind of solutions, no I don't. – KraZug Apr 29 at 5:03
• Thanks @AlexTrounev. In these equations ep and d are damping factors, when "ep" ="d" =0, svals ={0. + 0. I, 0. + 0. I, 0. + 0. I, 0. + 0. I} – Danial Apr 29 at 5:49
• @Danial, Do you know what eigenvalues to expect when $\epsilon$ and $d$ are both zero? – KraZug Apr 29 at 5:53

eqn = {D[V[t, x], t] == (gpar[x] ρ[t, x]/ρ0[x]) -
(D[ρ[t, x] + T[t, x], x]/(γ ρ0[x])) + 4/3 ϵ D[V[t, x], x, x],
D[ρ[t, x], t] == -D[ρ0[x] V[t, x], x],
D[T[t, x], t] == D[ρ0[x] V[t, x], x] - V[t, x] p0'[x] - γ p0[x] D[V[t, x], x]
+ d γ D[T[t, x], x, x]}


As I mentioned in a comment, note that there is only 5 spatial derivatives in the unknown variables, and that the highest derivatives for $$\rho'$$ and $$V''$$ both appear together in the same equation.

Assume the ansatz of $$a(t,x) = e^{\lambda t} a(x)$$ for the three functions $$\rho, V, T$$:

eqnSubbed = Simplify[eqn /. {ρ -> Function[{t, x}, Exp[λ t] ρx[x]],
T -> Function[{t, x}, Exp[λ t] Tx[x]], V -> Function[{t, x}, Exp[λ t] Vx[x]]},
{t > 0, x > 0, ρx[x] > 0, γ > 0}] /. {ρx -> ρ, Tx -> T, Vx -> V};


And then we can solve the second equation for $$\rho$$, to reduce to just a pair of second order equations for $$V$$ and $$T$$:

eqnSubbed2 = Simplify[eqnSubbed[[{1, 3}]] /. DSolve[eqnSubbed[[2]], ρ, x][[1]]];


Now 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. Note that this is how I noticed that your system was essentially a DAE, as applying it directly doesn't work (something for me to catch and fix, or at least error gracefully):

Needs["CompoundMatrixMethod"]

subs = {γ -> 1.67, ϵ -> 0.01, d -> 0.05, H -> 1.1, gpar -> Function[{x}, 0.38],
ρ0 -> Function[{x}, Exp[-x/H]], p0 -> Function[{x}, Exp[-x/H]]};
sys = ToMatrixSystem[eqnSubbed2, {T[0] == 0, V[0] == 0, T[1] == 0, V[1] == 0}, {T,
V}, {x, 0, 1}, λ] //. subs


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 $$\lambda$$; this is an analytic function whose roots coincide with eigenvalues of the original equation.

FindRoot will then find solutions for a given start point:

root = FindRoot[Evans[λ, sys], {λ, -1 + 3 I}]
(* {λ -> -0.342689 + 2.65899 I} *)


And we can plot the Evans function to see there are a set of real, negative eigenvalues as well as the complex ones,

Plot[Evans[λ, sys], {λ, -10, 10}]


I can't immediately see any eigenvalues with a positive real part, and suspect there aren't any for this case. For your function gpara, if you manage to read it in and make an interpolation function from the data correctly, the code above should be able to incorporate that.

Unfortunately my code doesn't currently get the eigenfunctions out immediately at the moment. However, you can use NDSolve to get them once you find a root:

  sol = NDSolve[Join[eqnSubbed2 //. subs /. root,
{T[0] == 0, V[0] == 0, T'[0] == 1, V[1] == 0}], {T, V}, {x, 0, 1}];
GraphicsRow[{Plot[Evaluate[ReIm@T[x] /. sol], {x, 0, 1}, PlotLabel -> "T(x)"],
Plot[Evaluate[ReIm@V[x] /. sol], {x, 0, 1}, PlotLabel -> "V(x)"],
Plot[Evaluate[ReIm[ρ[x] /. DSolve[eqnSubbed[[2]], ρ, x][[1]] //.subs /.root /.sol]]
, {x, 0, 1}, PlotLabel -> "ρ(x)"]}, ImageSize -> 1000]


• Thank you @KraZug for complete analysis. I need to go through the method you suggested. What we need to get correct eigenvalues and corresponding eigen vectors, and also how to check the eigen values are correct? – Danial Apr 29 at 7:38
• Eigenvectors? You mean eigenfunctions? I added an example of those. Note that $\rho$ does not meet the boundary conditions you were trying to impose at either end. – KraZug Apr 29 at 8:49
• Thank you very much @KraZug. what should be correct BC for \rho? – Danial Apr 29 at 9:23
• That is up to you to determine with your knowledge of the underlying system. Potentially nothing if you know the other two are correct. I think you can only apply 4 spatial boundary conditions though, due to the highest derivatives appearing in combination. – KraZug Apr 29 at 9:59
• Thanks@KraZug. I m using v11.3 and above code is showing errors related to Evans and FindRoot, even after quiting the kernel. – Danial Apr 29 at 10:08