(* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% *)
(* File: stabilit.m *)
(* Created: October 25, 1991 *)
(* Revised: December 10, 1992 *)
(* *)
(* Tested with the aid of Mathematica Version 2.0 for Microsoft Windows. *)
(* ======================================================================== *)
(* *)
(* Stability Analysis of Matrices and Polynomials *)
(* Based on Routh Criterion and Iteration Method *)
(* *)
(* *)
(* by Alexander Urintsev *)
(* *)
(* ======================================================================== *)
(* PURPOSE: *)
(* *)
(* This package includes functions for stability analysis of square *)
(* numerical matrices and polynomials with numerical coefficients. *)
(* Specifically, the amount of roots z satisfying the unequality *)
(* Re[z] < bound (or Abs[z] < bound) can be computed avoiding the *)
(* explicit finding of roots; number bound is supplied by the user. *)
(* *)
(* *)
(* Remarks: *)
(* *)
(* 1. Function StabilQ uses for stability analysis the criterion *)
(* due Routh. *)
(* *)
(* 2. The method used in function CountSpecialRoots is described in *)
(* the paper: *)
(* *)
(* V.L. Zaguskin, A.V. Kharitonov, Iteration method to solve the *)
(* problem of stability. Journal of Computational Mathematics and *)
(* Mathematical Physics, 1963, Volume 3, Number 2, pp. 361-364 *)
(* (in Russian). *)
(* ======================================================================== *)
(* *)
(* About the Author *)
(* *)
(* Dr. Alexander Urintsev works at the Joint Institute for Nuclear *)
(* Research as a senior research associate. He is 46 years old. He has *)
(* been working in applied research for 23 years. Dr. Urintsev is the *)
(* author/coauthor of about 40 scientific publications in the field of *)
(* applied and computational mathematics, hydrodynamic stability, *)
(* self-organization in nonlinear dynamical systems, application of *)
(* computer algebra systems in the problems of nonlinear magnetic optics *)
(* and accelerator physics. *)
(* *)
(* During 11 years Dr. Urintsev has been combining research work with *)
(* teaching at the Mechanical-Mathematical Department of the University *)
(* (Department of Computational Mathematics, Department of Functions and *)
(* Functional Analysis Theory): lecturing on programming, numerical *)
(* methods, probability theory and mathematical statistics, application *)
(* of computer algebra systems in mechanics and mathematics, supervising *)
(* the students' research work, etc. *)
(* *)
(* Fields of interest: *)
(* development and application of complex computer systems such as *)
(* Mathematica, Maple, Reduce both in applied researches and in *)
(* educational purposes, working out new algorithms of symbolic and *)
(* numerical mathematics, solving differential and integral equations, *)
(* data processing, producing of the special applied software, *)
(* programming in Mathematica, Maple, Reduce, C, Lisp, Fortran, Basic. *)
(* *)
(* Dr. Urintsev wants to apply for a job abroad. He would like to get *)
(* a position in a serious research institution, governmental *)
(* organization, computer firm producing software as a researcher, *)
(* software engineer or a programmer, or as a teacher at university or *)
(* at a college. *)
(* *)
(* Dr. Urintsev would be immensely grateful to colleagues for information *)
(* about possible vacancies, any useful advice, or for kind assistance. *)
(* Necessary documents (detailed CV, list of publications, degrees *)
(* copies, letters of recommendation, application) are available for *)
(* inquiry. *)
(* *)
(* December 1992, JINR, Dubna, RUSSIA. *)
(* *)
(* *)
(* Postal address: *)
(* *)
(* Dr. Alexander Urintsev *)
(* Department of Heavy Quarks *)
(* Laboratory of Particle Physics *)
(* Joint Institute for Nuclear Research *)
(* Head Post Office, P.O. Box 79 *)
(* 101000, Moscow *)
(* RUSSIA *)
(* *)
(* Fax: (7095) 975 23 81 *)
(* Telex: 911621 DUBNA SU *)
(* *)
(* Email: urintsev@lhe24.jinr.dubna.su *)
(* *)
(* ======================================================================== *)
BeginPackage["Stability`"]
StabilQ::usage =
" StabilQ[p, regimen] in case regimen=1 gives True\n
if Re[z]<0 for all roots z of polynomial p and False\n
otherwise. In case regimen=2 StabilQ[p, regimen] yields\n
True if Abs[z]<1 for all roots z and False otherwise.\n
Argument regimen is optional, by default regimen=1.\n
Polynomial p can be only in one variable and must have\n
numerical (rational, real or complex) coefficients. Instead\n
of polynomial you can also put in a square matrix with\n
numerical (rational, real or complex) elements. In this\n
case the part of the above mentioned polynomial is taken\n
actually by characteristic polynomial of this matrix; it is\n
calculated automatically inside StabilQ. The part of the\n
above z-roots are taken now by eigenvalues of the current\n
matrix."
CountSpecialRoots::usage =
" CountSpecialRoots[p, bound, regimen] gives in case\n
regimen=1 an amount of roots z of polynomial p satisfying\n
the unequality Re[z] 0,
If[ a1[[1]] <= 0, stabil = False; Break[]];
l1 = If[EvenQ[l], l/2, (l - 1)/2];
a1k = a1[[1]];
s = a0[[1]]/a1k;
Do[
a0[[k]] = a1k;
a1k = a1[[k + 1]];
a1[[k]] = a0[[k + 1]] - s a1k,
{k, 1, l1}
];
If[ 2 l1 != l, a0[[l1 + 1]] = a1k; a1[[l1 + 1]] = 0 ];
l = l - 1;
];
stabil
] /; (SameQ[Head[p], List] ||
(Not[SameQ[Head[p], List]] && Length[Variables[p]] == 1)) &&
(regimen == 1 || regimen == 2)
CountSpecialRoots[p_, bound_?NumberQ, regimen_Integer:1] :=
Module[{s, x, boundn, a, n, err, b, iter, n1, k, j, se, bj, nextiter,
aj, aj1, ak, testerr, numeric,
maxiter = 17, eps1 = 10^-15, eps2 = 10^-3},
If[ Not[SameQ[Head[p], List]], Goto[LabPoly]];
j = Dimensions[p];
If[Length[j] != 2, Goto[LabMes]];
If[ j[[1]] == j[[2]], Goto[LabMatr]];
Label[LabMes];
Message[CountSpecialRoots::arg1];
Return[Null];
Label[LabMatr];
n = j[[1]];
numeric = True;
Clear[j];
Do[ If[ Not[NumberQ[p[[j, k]] ]], numeric = False], {j, n}, {k, n}];
Clear[j, k];
If[numeric, Goto[LabDet]];
Message[CountSpecialRoots::arg1];
Return[Null];
Label[LabDet];
s = Det[p - x IdentityMatrix[n]];
a = CoefficientList[s, x];
Goto[LabTrans];
Label[LabPoly];
s = Expand[p];
x = Variables[s][[1]];
a = CoefficientList[s, x];
numeric = True;
Do[If[Not[NumberQ[a[[k]]]], numeric = False], {k, Length[a]}];
If[numeric, Goto[LabTrans]];
Message[CountSpecialRoots::arg1];
Return[Null];
Label[LabTrans];
boundn = N[bound];
a = N[a];
If[regimen == 2,
(* then *)
If[boundn != 1., s = s /. x -> boundn x; a = CoefficientList[s, x]];
a = Reverse[a];
n = Length[a] - 1;
Do[
ak = a[[k + 1]];
aj = 0;
Do[
aj1 = a[[j]];
a[[j + 1]] = aj1 + aj;
aj = aj1,
{j, k, 1, -1}
];
a[[1]] = a[[1]] + ak,
{k, n}
];
ak = 2;
Do[
aj = ak a[[k + 1]];
Do[
aj1 = a[[j]];
a[[j + 1]] = aj - aj1;
aj = aj1,
{j, k, 1, -1}
];
ak = ak + ak,
{k, n}
];
a = Reverse[a];
k = n;
While[ a[[k + 1]] == 0 && k != 0, k = k - 1];
a = Take[a, k + 1];
Clear[j, k],
(* else *)
If[boundn != 0.,
s = Expand[s /. x -> x + boundn];
a = CoefficientList[s, x];
];
];
a = N[a];
n = Length[a] - 1;
err = N[Sqrt[eps1 n]];
b = Table[0, {n + 1}];
iter = 0;
Label[LabIter];
n1 = n + 1;
Clear[k, j];
Do[
s = Sum[a[[k - j]] a[[k]], {k, j+1, n1}];
b[[j+1]] = s,
{j, 0, n}
];
b[[1]] = b[[1]]/2;
Do[
bk = b[[k + 1]];
b[[k - 1]] = b[[k - 1]] - bk;
b[[k + 1]] = 2 b[[k + 1]],
{j, 2, n}, {k, n, j, -1}
];
k = n;
While[ b[[k + 1]] == 0 && k != 0, k = k - 1];
b = Take[b, k + 1];
If[k == 0, Return[0]];
s = 1/b[[k + 1]];
Clear[j];
se = 0;
Do[ bj = s b[[j]];
se = se + Abs[a[[j]] - bj];
a[[j]] = bj,
{j, k, 1, -1}
];
a[[k + 1]] = 1;
iter = iter + 1;
nextiter = False;
If[kerr;
If[nextiter && iter0))
End[] (* end Private Context *)
EndPackage[] (* end Stability Context *)
(* %%%%%%%%%%%%%%%%%%%%% The end of stabilit.m file %%%%%%%%%%%%%%%%%%%%% *)