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