The Routh-Hurwitz criterion in Mathematica implementation

Question

The Routh-Hurwitz criterion is a mathematical test that provides conditions for the stability of a system by analyzing the signs of the determinants of specific matrices derived from the system's characteristic polynomial.

I successfully ported the Python code to Mathematica code. I would love any feedback you might have.

(* Function to calculate the determinant of a matrix *)
(* We can use Mathematica's built-in Det function instead of implementing it ourselves *)
determinant[matrix_] := Det[matrix]
(* Alternatively, if we want to keep the recursive implementation: *)
determinantRecursive[matrix_] := Module[{n = Length[matrix]},
 If[n == 1,
 matrix[[1, 1]],
 Sum[
 (-1)^(i-1) * matrix[[1, i]] * determinantRecursive[Drop[#, {1}, {i}]& /@ Drop[matrix, {1}]],
 {i, 1, n}
 ]
 ]
]
(* Function to calculate determinants for each power from 1 to max_power *)
calculateDeterminants[matrix_, maxPower_] := Table[
 Det[matrix[[1 ;; power, 1 ;; power]]],
 {power, 1, maxPower}
]
(* Function to check if all coefficients are positive *)
checkPositiveCoefficients[coefficients_] := AllTrue[coefficients, # > 0 &]
(* Function to apply the Routh-Hurwitz criterion *)
criterionRouthHurwitz[matrix_] := Module[{maxPower, determinants},
 maxPower = Length[matrix];
 determinants = calculateDeterminants[matrix, maxPower];
 
 If[AllTrue[determinants, # > 0 &],
 Print["The system is stable."],
 If[AnyTrue[determinants, # == 0 &],
 Print["The stability can't be recognized."],
 Print["The system is unstable."]
 ]
 ]
]
(* Main function *)
main[] := Module[{s, M, coefficients, sizePower, matrix},
 s = Symbol["s"];
 M = s^4 + 6*s^3 + 13*s^2 + 12*s + 4;
 
 (* Extract coefficients from the polynomial M *)
 coefficients = CoefficientList[M, s];
 (* Reverse to match Python's all_coeffs order *)
 coefficients = Reverse[coefficients];
 
 (* Determine the size of the matrix *)
 sizePower = Length[coefficients] - 1;
 
 (* Initialize the matrix - fixed implementation *)
 matrix = Table[
 Module[{index},
 index = sizePower - 3 - (i-1) + 2*(j-1);
 If[index >= 0 && index < Length[coefficients],
 coefficients[[index + 1]],
 0
 ]
 ],
 {i, 1, sizePower},
 {j, 1, sizePower}
 ];
 
 (* Print the matrix *)
 Print["Size power: ", sizePower];
 Print["Matrix:"];
 Do[
 Print[matrix[[i]]],
 {i, 1, sizePower}
 ];
 
 If[checkPositiveCoefficients[coefficients],
 criterionRouthHurwitz[matrix],
 Print["The system is unstable due to negative coefficients."]
 ]
]
(* Call the main function *)
main[]

Answer 1

automated tests

We're doing a "refactor" here -- changing the code while not changing the behavior. Sometimes that happens within a language, but here we're making mathematica code behave the same way that python code behaves.

Having a solid test suite is pretty much a prerequisite to doing a refactor. It instills confidence in the newly written code.

We might even have to port some python unit tests to mathematica, or ask them to fork off a child process to run the mathematica target code and evaluate whether the result is correct.

The main driver is very nice and appreciated. But it's not self evaluating. It doesn't know the right answer -- it demands that we eyeball the results each time, and manually decide whether a correct or incorrect answer was displayed.

type annotations

Throwing in some annotations wouldn't hurt.
Seeing mypy --strict *.py lint cleanly is another way to instill greater confidence in a codebase.

redundant functions

determinant[matrix_] := Det[matrix]
determinantRecursive[matrix_] := ...

Both of these are very nice. Please pick one.

(Or exercise both of them from unit tests, to demonstrate that they behave identically.)

similar name, similar concept

This original python code does not appear in the OP:

 determinants = []
 ...
 det = determinant(submatrix)

The code works fine. But there is an unfortunate shifting of gears when we move between singular and plural. Consider the idiom for x in xs:, where x has some type and xs is a container of that type. But here the singular is a float-returning function, and the plural is a container whose elements have its return type.

Each name works fine, and was well chosen in isolation. But the combination comes out slightly jarring. Maybe we should have dets: list[float]?

identically zero

This line

 If[AnyTrue[determinants, # == 0 &],

exactly matches the original source line of

 has_zero = any(det == 0 for det in determinants)

The determinant is an IEEE-754 FP quantity, the sum of several products. Ask any numerical analyst: round off error "is a thing!". I am not familiar with the business problem you hope to solve, but I'm skeptical that a bit pattern of 64 zeros is exactly what your use case calls for. Consider choosing some small epsilon, and testing whether abs(det) is less than that threshold. Essentially we would be asking "is det small?", rather than "is det exactly zero?".