Differential equation system, Jacobian matrix, characteristic equation

Question

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[3]), 
Y -> -((Sqrt[β] Sqrt[-1 + 3 ρ])/Sqrt[3]), 
Z -> 1/3 (-1 + 3 ρ)}, {X -> (Sqrt[β] Sqrt[-1 + 3 ρ])/Sqrt[3], 
Y -> (Sqrt[β] Sqrt[-1 + 3 ρ])/Sqrt[3], Z -> 1/3 (-1 + 3 ρ)}}[[2]]

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

a_ 3 = First[{{X -> 0, Y -> 0, 
Z -> 0}, {X -> -((Sqrt[β] Sqrt[-1 + 3 ρ])/Sqrt[3]), 
Y -> -((Sqrt[β] Sqrt[-1 + 3 ρ])/Sqrt[3]), 
Z -> 1/3 (-1 + 3 ρ)}, {X -> (Sqrt[β] Sqrt[-1 + 3 ρ])/Sqrt[3], 
Y -> (Sqrt[β] Sqrt[-1 + 3 ρ])/Sqrt[3], 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?

Answer 1

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

Answer 2

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][[1]]
P2 = Solve[F == 0, V][[2]]
P3 = Solve[F == 0, V][[3]]

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]

Answer 3

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[[1]]
j /. eq[[2]]
j /. eq[[3]]

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

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

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