(* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% *) (* 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 %%%%%%%%%%%%%%%%%%%%% *)