Differential equation system, Jacobian matrix, characteristic equation

We assume the differential equation system: $$X'=σ(Y-X),Y'=X(ρ-Z)-Y,Z'=XY-βZ$$. I can calculate the equilibrium points as:

Solve[σ (Y - X) == 0 && 3 X (ρ - Z) - Y == 0 && X Y - β Z == 0, {X, Y,Z}]

a_ 1 = {{X -> 0, Y -> 0,
Z -> 0}, {X -> -((Sqrt[β] Sqrt[-1 + 3 ρ])/Sqrt),
Y -> -((Sqrt[β] Sqrt[-1 + 3 ρ])/Sqrt),
Z -> 1/3 (-1 + 3 ρ)}, {X -> (Sqrt[β] Sqrt[-1 + 3 ρ])/Sqrt,
Y -> (Sqrt[β] Sqrt[-1 + 3 ρ])/Sqrt, Z -> 1/3 (-1 + 3 ρ)}}[]

a_ 2 = Last[{{X -> 0, Y -> 0,
Z -> 0}, {X -> -((Sqrt[β] Sqrt[-1 + 3 ρ])/Sqrt),
Y -> -((Sqrt[β] Sqrt[-1 + 3 ρ])/Sqrt),
Z -> 1/3 (-1 + 3 ρ)}, {X -> (Sqrt[β] Sqrt[-1 + 3 ρ])/Sqrt,
Y -> (Sqrt[β] Sqrt[-1 + 3 ρ])/Sqrt, Z -> 1/3 (-1 + 3 ρ)}}]

a_ 3 = First[{{X -> 0, Y -> 0,
Z -> 0}, {X -> -((Sqrt[β] Sqrt[-1 + 3 ρ])/Sqrt),
Y -> -((Sqrt[β] Sqrt[-1 + 3 ρ])/Sqrt),
Z -> 1/3 (-1 + 3 ρ)}, {X -> (Sqrt[β] Sqrt[-1 + 3 ρ])/Sqrt,
Y -> (Sqrt[β] Sqrt[-1 + 3 ρ])/Sqrt, Z -> 1/3 (-1 + 3 ρ)}}]

I have to calculate the Jacobian matrix for each of the three equilibrium point and then their characteristic polyonymial. My approach is to calculate separate a matrix with the components of the Jacobian matrix and the put $$a_1, a_2, a_3$$. Is there a different method?

• One quick thing, you can't define variables as a_ 1, etc. Easier to save the results of Solve with eq = Solve[... then take parts of eq. Nov 23 '19 at 14:59
• I am a beginner, so it would be helpful (for me) to be more analytic Nov 23 '19 at 15:01

Once we have the equilibrium points as

equs = {sigma (X - Y), X (rho - Z) - Y, X Y - beta Z}
sols = Solve[equs == 0, {X, Y, Z}]

you can calculate the associated jacobian to each equilibrium point as

Js = Grad[equs, {X, Y, Z}] /. sols

and the corresponding characteristic polynomials as

Table[CharacteristicPolynomial[Js[[k]], lambda], {k, 1, Length[Js]}]

or if you know the value of the parameters involved

parms = {sigma -> 10, rho -> 28, beta -> 8/3};
eigen = Table[Eigenvalues[Js[[k]] /. parms], {k, 1, Length[Js]}] // N;
eigen  // MatrixForm
• Huh, didn't know you could just use 0 in Solve like that! Nov 23 '19 at 15:18
• nice, may I ask something? for specific values: $β=\frac{8}{3}, σ=10, ρ=28$ can we calculate $λ$ and then use Eigenvalue command? Nov 23 '19 at 15:35
• The eigenvalues calculation was included. Nov 23 '19 at 15:45
• mathematica.stackexchange.com/questions/253201/…
– dtn
Aug 12 '21 at 4:02

The vector field defined by your system:

f1 = σ (Y - X);
f2 = 3 X (ρ - Z) - Y;
f3 = X Y - β Z;
F = {f1, f2, f3};
V = {X, Y, Z};

The equilibria:

P1 = Solve[F == 0, V][]
P2 = Solve[F == 0, V][]
P3 = Solve[F == 0, V][]

The Jacobian matrix:

J = FullSimplify[D[F, {V}]]
MatrixForm[J]

The three linear approximations:

J1 = FullSimplify[J /. P1];
MatrixForm[J1]
J2 = FullSimplify[J /. P2];
MatrixForm[J2]
J3 = FullSimplify[J /. P3];
MatrixForm[J3]

The three characteristic polynomials:

pol1 = -Collect[CharacteristicPolynomial[J1, λ], λ, Simplify]
pol2 = -Collect[CharacteristicPolynomial[J2, λ], λ, Simplify]
pol3 = -Collect[CharacteristicPolynomial[J3, λ], λ, Simplify]

Let's start by making a list of the equations:

eqns = {σ (Y - X), 3 X (ρ - Z) - Y, X Y - β Z};

Then solve for the equilibria and save the result in eq:

eq = Solve[eqns == {0, 0, 0}, {X, Y, Z}]

Make the generic Jacobian:

j = D[eqns, {{X, Y, Z}}]

and then you can evaluate it at particular equilibria using /.:

j /. eq[]
j /. eq[]
j /. eq[]

Once you have those, it's easy to find the characteristic polynomial or eigenvalues:

CharacteristicPolynomial[j /. eq[], λ]
Eigenvalues[j /. eq[]

(*
(-β - λ) (λ + λ^2 + σ + λ σ - 3ρ σ)
{-β, 1/2 (-1 - σ - Sqrt[1 - 2σ + 12ρ σ + σ^2]),
1/2 (-1 - σ + Sqrt[1 - 2σ + 12ρ σ + σ^2])}
*)