home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / l / l042 / 1.ddi / CHAP2.ARC / LAGUERRE.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1987-12-30  |  9.5 KB  |  294 lines
  1. program Laguerre_Prog;
  2.  
  3. {---------------------------------------------------------------------------}
  4. {-                                                                         -}
  5. {-     Turbo Pascal Numerical Methods Toolbox                              -}
  6. {-     Copyright (c) 1986, 87 by Borland International, Inc.               -}
  7. {-                                                                         -}
  8. {-               Purpose: This program demonstrates Laguerre's method      -}
  9. {-                                                                         -}
  10. {-               Unit   : RootsEqu    procedure Laguerre                   -}
  11. {-                                                                         -}
  12. {---------------------------------------------------------------------------}
  13.  
  14. {$I-}                { Disable I/O error trapping }
  15. {$R+}                { Enable range checking }
  16.  
  17. uses
  18.   RootsEqu, Dos, Crt, Common;
  19.  
  20. var
  21.   Guess : TNcomplex;                 { Initial approximation }
  22.   InitDegree, Degree : integer;      { Degree of polynomial }
  23.   InitPoly, Poly : TNCompVector;     { Coefficients of polynomial }
  24.   Tol : Float;                       { Tolerance in answer }
  25.   Iter : TNIntVector;                { Number of iterations to find each root }
  26.   MaxIter : integer;                 { Max number of iterations allowed }
  27.   NumRoots : integer;                { Number of roots found }
  28.   Answer : TNCompVector;             { Roots }
  29.   yAnswer : TNCompVector;            { Function evaluated at roots }
  30.   Error : byte;                      { Flags an error }
  31.  
  32. procedure Initial(var Degree  : integer;
  33.                   var Poly    : TNCompVector;
  34.                   var Guess   : TNcomplex;
  35.                   var Tol     : Float;
  36.                   var MaxIter : integer);
  37.  
  38. {-----------------------------------------------------------}
  39. {- Output: Degree, Poly, Guess, Tol, MaxIter               -}
  40. {-                                                         -}
  41. {- This procedure initializes the above variables to zero. -}
  42. {-----------------------------------------------------------}
  43.  
  44. begin
  45.   Degree := 0;
  46.   FillChar(Poly, SizeOf(Poly), 0);
  47.   FillChar(Guess, SizeOf(Guess), 0);
  48.   Tol := 0;
  49.   MaxIter := 0
  50. end; { procedure Initial }
  51.  
  52. procedure UserInput(var Degree  : integer;
  53.                     var Poly    : TNCompVector;
  54.                     var Guess   : TNcomplex;
  55.                     var Tol     : Float;
  56.                     var MaxIter : integer);
  57.  
  58. {---------------------------------------------------------------}
  59. {- Output: Degree, Poly, Guess, Tol, MaxIter                   -}
  60. {-                                                             -}
  61. {- This procedure assigns values to the above variables from   -}
  62. {- keyboard input.  The Degree of the polynomial (Degree), the -}
  63. {- coefficients of the polynomial (Poly), the initial guess    -}
  64. {- (Guess), the tolerance (Tol), and the maximum number of     -}
  65. {- iterations (MaxIter) are all read in here.                  -}
  66. {---------------------------------------------------------------}
  67.  
  68. var
  69.   Ch : char;
  70.  
  71. procedure GetCoefficientsFromKeyboard(var Degree : integer;
  72.                                       var Poly   : TNCompVector);
  73.  
  74. {-------------------------------------------------}
  75. {- Output: Degree, Poly                          -}
  76. {-                                               -}
  77. {- The Degree and coefficients of the polynomial -}
  78. {- are read in from the keyboard.                -}
  79. {-------------------------------------------------}
  80.  
  81. var
  82.   Term : integer;
  83. begin
  84.   Writeln;
  85.   repeat
  86.     Write('Degree of the polynomial (<= ',TNArraySize,')? ');
  87.     Readln(Degree);
  88.     IOCheck;
  89.   until (Degree > 1) and (Degree <= TNArraySize) and not IOerr;
  90.   Writeln;
  91.   Writeln('Input the complex coefficients of the polynomial');
  92.   Writeln('where Poly[n] is the coefficients of x^n');
  93.   Writeln;
  94.   for Term := Degree downto 0 do
  95.   begin
  96.     repeat
  97.       Write('Re(Poly[',Term,']) = ');
  98.       Readln(Poly[Term].Re);
  99.       IOCheck;
  100.     until not IOerr;
  101.     repeat
  102.       Write('Im(Poly[',Term,']) = ');
  103.       Readln(Poly[Term].Im);
  104.       IOCheck;
  105.     until not IOerr;
  106.     Writeln;
  107.   end;
  108. end; { procedure GetCoefficientsFromKeyboard }
  109.  
  110. procedure GetCoefficientsFromFile(var Degree : integer;
  111.                                   var Poly   : TNCompVector);
  112.  
  113. {------------------------------------------------------}
  114. {- Output: Degree, Poly                               -}
  115. {-                                                    -}
  116. {- This procedure reads in values for Degree and Poly -}
  117. {- from a text file.                                  -}
  118. {------------------------------------------------------}
  119.  
  120. var
  121.   Term : integer;
  122.   InFile : text;
  123.   FileName : string[255];
  124.  
  125.  begin
  126.    repeat
  127.      Writeln;
  128.      repeat
  129.        Write('File name? ');
  130.        Readln(FileName);
  131.        Assign(InFile, FileName);
  132.        Reset(InFile);
  133.        IOCheck;
  134.      until not IOerr;
  135.      Degree := 0;
  136.      Read(InFile, Degree);
  137.      IOCheck;
  138.      if not(Degree in [0..TNArraySize]) and not IOerr then
  139.        Writeln('Degree too big');
  140.      Term := Degree;
  141.      while (not IOerr) and (Term >= 0) do
  142.      begin
  143.        Read(InFile, Poly[Term].Re);
  144.        Read(InFile, Poly[Term].Im);
  145.        IOCheck;
  146.        Term := Pred(Term);
  147.      end;
  148.    until not IOerr;
  149.    Close(InFile);
  150.  end; { procedure GetCoefficientsFromFile }
  151.  
  152. procedure GetInitialGuess(var Guess : TNcomplex);
  153. begin
  154.   Writeln;
  155.   Writeln('Initial approximation to the root: ');
  156.   repeat
  157.     Write('Re(Approximation) = ');
  158.     Readln(Guess.Re);
  159.     IOCheck;
  160.   until not IOerr;
  161.   repeat
  162.     Write('Im(Approximation) = ');
  163.     Readln(Guess.Im);
  164.     IOCheck;
  165.   until not IOerr;
  166. end; { procedure GetInitialGuess }
  167.  
  168. procedure GetTolerance(var Tol : Float);
  169. begin
  170.   Writeln;
  171.   repeat
  172.     Tol := 1E-8;
  173.     Write('Tolerance (> 0): ');
  174.     ReadFloat(Tol);
  175.     IOCheck;
  176.     if Tol <= 0 then
  177.     begin
  178.       IOerr := true;
  179.       Tol := 1E-8;
  180.     end;
  181.   until not IOerr;
  182. end; { procedure GetTolerance }
  183.  
  184. procedure GetMaxIter(var MaxIter : integer);
  185. begin
  186.   Writeln;
  187.   repeat
  188.     MaxIter := 100;
  189.     Write('Maximum number of iterations (> 0): ');
  190.     ReadInt(MaxIter);
  191.     IOCheck;
  192.     if MaxIter < 0 then
  193.     begin
  194.       IOerr := true;
  195.       MaxIter := 100;
  196.     end;
  197.   until not IOerr;
  198. end; { procedure GetMaxIter }
  199.  
  200. begin { procedure UserInput }
  201.   case InputChannel('Input Data From') of
  202.     'K' : GetCoefficientsFromKeyboard(Degree, Poly);
  203.     'F' : GetCoefficientsFromFile(Degree, Poly);
  204.   end;
  205.   GetInitialGuess(Guess);
  206.   GetTolerance(Tol);
  207.   GetMaxIter(MaxIter);
  208.   GetOutputFile(OutFile);
  209. end; { procedure UserInput }
  210.  
  211. procedure Results(InitDegree : integer;
  212.                   InitPoly   : TNCompVector;
  213.                   Degree     : integer;
  214.                   Poly       : TNCompVector;
  215.                   Guess      : TNcomplex;
  216.                   Answer     : TNCompVector;
  217.                   yAnswer    : TNCompVector;
  218.                   Tol        : Float;
  219.                   Iter       : TNIntVector;
  220.                   Error      : byte);
  221.  
  222. {------------------------------------------------------------}
  223. {- This procedure outputs the results to the device OutFile -}
  224. {------------------------------------------------------------}
  225.  
  226. var
  227.   Term : integer;
  228.  
  229. begin
  230.   Writeln(OutFile, 'Initial polynomial:');
  231.   Writeln(OutFile);
  232.   for Term := InitDegree downto 0 do
  233.     Writeln(OutFile, 'InitPoly[',Term,']:',
  234.   InitPoly[Term].Re, ' +', InitPoly[Term].Im,'i');
  235.   Writeln(OutFile);
  236.   Writeln(OutFile);
  237.   Write(OutFile,'Initial approximation: ' : 30);
  238.   Writeln(OutFile,Guess.Re, ' +',Guess.Im,'i');
  239.   Writeln(OutFile,'Tolerance: ' : 30, Tol);
  240.   Writeln(OutFile,'Maximum number of iterations: ' : 30, MaxIter);
  241.   Writeln(OutFile);
  242.   if Error = 1 then
  243.     DisplayWarning;
  244.   if Error >= 2 then
  245.     DisplayError;
  246.  
  247.   case Error of
  248.     1 : Writeln(OutFile,'This will take more than ', MaxIter,' iterations.');
  249.  
  250.     2 : Writeln(OutFile,
  251.                 'The degree of the polynomial must be greater than zero.');
  252.  
  253.     3 : Writeln(OutFile,'The tolerance must be greater than zero.');
  254.  
  255.     4 : Writeln(OutFile,
  256.                 'The maximum number of iterations must be greater than zero.');
  257.  
  258.   end; { case }
  259.   Writeln(OutFile);
  260.   if Degree > 0 then
  261.   begin
  262.     Writeln(OutFile,'The deflated polynomial:');
  263.     for Term := Degree downto 0 do
  264.       Writeln(OutFile, 'Poly[',Term,']:',
  265.                         Poly[Term].Re, ' +', Poly[Term].Im,'i');
  266.   end;
  267.  
  268.   if Error <= 1 then
  269.     for Term := 1 to NumRoots do
  270.     begin
  271.       Writeln(OutFile);
  272.       Writeln(OutFile, 'Root ', Term);
  273.       Writeln(OutFile, 'Number of iterations: ':25, Iter[Term] : 3);
  274.       Write(OutFile,'Calculated root: ':25);
  275.       Writeln(OutFile, Answer[Term].Re, ' +', Answer[Term].Im, 'i');
  276.       Writeln(OutFile,'Value of the function at');
  277.       Write(Outfile,'the calculated root: ':25);
  278.       Writeln(OutFile, yAnswer[Term].Re, ' +', yAnswer[Term].Im, 'i');
  279.     end;
  280. end; { procedure Results }
  281.  
  282. begin { program Laguerre }
  283.   ClrScr;
  284.   Initial(InitDegree, InitPoly, Guess, Tol, MaxIter);
  285.   UserInput(InitDegree, InitPoly, Guess, Tol, MaxIter);
  286.   Degree := InitDegree;
  287.   Poly := InitPoly;
  288.   Laguerre(Degree, Poly, Guess, Tol, MaxIter,
  289.            NumRoots, Answer, yAnswer, Iter, Error);
  290.   Results(InitDegree, InitPoly, Degree, Poly, Guess,
  291.           Answer, yAnswer, Tol, Iter, Error);
  292.   Close(OutFile);
  293. end. { program Laguerre }
  294.