home *** CD-ROM | disk | FTP | other *** search
/ Education Sampler 1992 [NeXTSTEP] / Education_1992_Sampler.iso / Mathematics / Notebooks / COSY_PAK / COSYPAK / chap4.m < prev    next >
Text File  |  1992-06-01  |  2KB  |  76 lines
  1.  
  2. (* 
  3. **********************************************************
  4.  
  5.          COSY_PAK package chap4.m        
  6.  
  7. **********************************************************
  8. *)
  9.  
  10.  
  11. If [ TrueQ[ $VersionNumber >= 2.0 ],
  12.      Off[ General::spell ];
  13.      Off[ General::spell1 ] ]
  14.      
  15. (*  B E G I N     P A C K A G E  *)
  16.  
  17.  
  18. BeginPackage["COSYPAK`chap4`"]         
  19.  
  20. Routh::usage = " ";
  21.  
  22.  
  23. Begin["`Private`"]
  24.  
  25. Routh[charpoly_,s_,z_:zero] :=
  26. Module[{Clist,ncoeff,list1,list2,Routhtable,nnn,temp,jj,i,first,k},
  27.     Clist = Reverse[CoefficientList[charpoly,s]];
  28.     ncoeff = Length[Clist];
  29.     nnn = ncoeff;
  30.     
  31.     If[OddQ[ncoeff], ncoeff = ncoeff +1;   (* make Coefflist to be *) 
  32.     AppendTo[Clist,0]];              (* even number    *) 
  33.     list1 = {}; list2= {}; Routhtable = {};
  34.     Do[ AppendTo[list1,Clist[[2 k -1]]],{k,1,ncoeff/2}];    
  35.     Do[ AppendTo[list2,Clist[[2 k]]],{k,1,ncoeff/2}];
  36.     If[list2[[1]] == 0 , list2[[1]] = z ];
  37.     
  38.     Routhtable = Join[{list1},{list2}];
  39.  
  40.     
  41.     For[jj = nnn -2, jj > 0 , jj--,
  42.         temp={};
  43.         Do[elm=Simplify[(list2[[1]] list1[[i+1]]-list1[[1]] list2[[i+1]])/
  44.                         list2[[1]] ];
  45.            If[ (i == 1) && Abs[elm] < 0.00001 , elm= z ];        
  46.            AppendTo[temp,elm],
  47.           {i,1,ncoeff/2 - 1}];
  48.         AppendTo[temp,0];              
  49.         Routhtable = Append[Routhtable,temp];
  50.         list1 = list2; list2 = temp
  51.         ];
  52.         
  53.     first = Table[Routhtable[[i,1]],{i,1,Length[Routhtable]}];    
  54.     Print["The Routh's Table:"];
  55.     Print[" "];
  56.     Print[MatrixForm[Routhtable]];
  57.     Print[" "];
  58.     Print["The first column of the Routh's table is: "];
  59.     Print[" "];
  60.     Print[ first ];
  61.     Print[" "];
  62.     Do[ If[first[[i]] < 0, Print["Sign Change!! Unstable root(s) exist."]],
  63.         {i,1,Length[first]}];
  64.     Return[Routhtable]         
  65.     ]/;PolynomialQ[charpoly,s]
  66.  
  67. End[];
  68. EndPackage[];
  69.  
  70. (*  E N D     P A C K A G E  *)
  71.  
  72. If [ TrueQ[ $VersionNumber >= 2.0 ],
  73.      On[ General::spell ];
  74.      On[ General::spell1 ] ]
  75.  
  76.