home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / l / l042 / 1.ddi / CHAP2.ARC / MULLER.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1987-12-30  |  6.4 KB  |  187 lines
  1. program Muller_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 Muller's method.              -}
  9. {-                                                                         -}
  10. {-        Unit   : RootsEqu    procedure Muller                            -}
  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. type
  21.   TNCompVector = array[0..TNArraySize] of TNcomplex;
  22.  
  23. var
  24.   Guess : TNcomplex;               { Initial approximations }
  25.   Tol : Float;                     { Tolerance in answer }
  26.   Iter : integer;                  { Number of iterations }
  27.   MaxIter : integer;               { Max number of iterations allowed }
  28.   Answer, yAnswer : TNcomplex;     { Root and function evaluated at root }
  29.   Error : byte;                    { Flags an error }
  30.  
  31. {$F+}
  32. {-------------- HERE IS THE FUNCTION ------------------}
  33.  
  34. procedure TNTargetF(X : TNcomplex; var Y : TNcomplex);
  35.  
  36. {----------------------------------------------------------}
  37. {-                                                        -}
  38. {-  Here are some standard complex functions:             -}
  39. {-                                                        -}
  40. {-    Cos(X) = Cos(X.Re)-(Exp(-X.Im) + Exp(X.Im))/2 +     -}
  41. {-             (Sin(X.Re)-(Exp(-X.Im) - Exp(X.Im))/2)i    -}
  42. {-                                                        -}
  43. {-    Sin(X) = Sin(X.Re)-(Exp(-X.Im) + Exp(X.Im))/2 -     -}
  44. {-             (Cos(X.Re)-(Exp(-X.Im) - Exp(X.Im))/2)i    -}
  45. {-                                                        -}
  46. {-    Exp(X) = Exp(X.Re)-Cos(X.Im) + Exp(X.Re)-Sin(X.Im)i -}
  47. {-                                                        -}
  48. {-    Sqr(X) = Sqr(X.Re) + Sqr(X.Im) + i(2 - X.Re - X.Im) -}
  49. {-                                                        -}
  50. {----------------------------------------------------------}
  51. begin
  52.   Y.Re := Cos(X.Re) * (Exp(-X.Im) + Exp(X.Im)) / 2 - X.Re;
  53.   Y.Im := Sin(X.Re) * (Exp(-X.Im) - Exp(X.Im)) / 2 - X.Im;
  54.   { This is the complex function Y := Cos(X) - X }
  55. end; { procedure TNTargetF }
  56. {$F-}
  57.  
  58. procedure UserInput(var Guess   : TNcomplex;
  59.                     var Tol     : Float;
  60.                     var MaxIter : integer);
  61.  
  62. {----------------------------------------------------------}
  63. {- Output: Guess, Tol, MaxIter                            -}
  64. {-                                                        -}
  65. {- This procedure assigns values to the above variables   -}
  66. {- from keyboard input.  The initial approximation of the -}
  67. {- guess (Guess), tolerance (Tol), and maximum number of  -}
  68. {- iterations (MaxIter) are all input here.               -}
  69. {----------------------------------------------------------}
  70.  
  71. procedure GetInitialGuess(var Guess : TNcomplex);
  72. begin
  73.   Writeln;
  74.   Writeln('Initial approximation to the root:  ');
  75.   repeat
  76.     Write('Re(Approximation) = ');
  77.     Readln(Guess.Re);
  78.     IOCheck;
  79.   until not IOerr;
  80.   repeat
  81.     Write('Im(Approximation) = ');
  82.     Readln(Guess.Im);
  83.     IOCheck;
  84.   until not IOerr;
  85. end; { procedure GetInitialGuess }
  86.  
  87. procedure GetTolerance(var Tol : Float);
  88. begin
  89.   Tol := 1E-8;
  90.   Writeln;
  91.   repeat
  92.     Write('Tolerance (> 0): ');
  93.     ReadFloat(Tol);
  94.     IOCheck;
  95.     if Tol <= 0 then
  96.     begin
  97.       IOerr := true;
  98.       Tol := 1E-8;
  99.     end;
  100.   until not IOerr;
  101. end; { procedure GetTolerance }
  102.  
  103. procedure GetMaxIter(var MaxIter : integer);
  104. begin
  105.   MaxIter := 100;
  106.   Writeln;
  107.   repeat
  108.     Write('Maximum number of iterations (> 0): ');
  109.     ReadInt(MaxIter);
  110.     IOCheck;
  111.     if MaxIter < 0 then
  112.     begin
  113.       IOerr := true;
  114.       MaxIter := 100;
  115.     end;
  116.   until not IOerr;
  117. end; { procedure GetMaxIter }
  118.  
  119. begin { procedure UserInput }
  120.   GetInitialGuess(Guess);
  121.   GetTolerance(Tol);
  122.   GetMaxIter(MaxIter);
  123.   GetOutputFile(OutFile);
  124. end; { procedure UserInput }
  125.  
  126. procedure Results(Guess   : TNcomplex;
  127.                   Answer  : TNcomplex;
  128.                   yAnswer : TNcomplex;
  129.                   Tol     : Float;
  130.                   MaxIter : integer;
  131.                   Iter    : integer;
  132.                   Error   : byte);
  133.  
  134. {------------------------------------------------------------}
  135. {- This procedure outputs the results to the device OutFile -}
  136. {------------------------------------------------------------}
  137.  
  138. begin
  139.   Writeln(OutFile);
  140.   Writeln(OutFile);
  141.   Write(OutFile, 'Initial approximation: ' : 30);
  142.   Writeln(OutFile, Guess.Re, ' +', Guess.Im, 'i');
  143.   Writeln(OutFile, 'Tolerance: ' : 30, Tol);
  144.   Writeln(OutFile, 'Maximum number of iterations: ' : 30, MaxIter);
  145.   Writeln(OutFile);
  146.   if Error in [1, 2] then
  147.     DisplayWarning;
  148.   if Error >= 3 then
  149.     DisplayError;
  150.  
  151.   case Error of
  152.     1 : Writeln(OutFile, 'This will take more than ', MaxIter, ' iterations.');
  153.  
  154.     2 : begin
  155.           Writeln(OutFile, 'A parabola which intersects the x-axis can not');
  156.           Writeln(OutFile, 'be constructed through these three points.');
  157.         end;
  158.  
  159.     3 : Writeln(OutFile, 'The tolerance must be greater than zero.');
  160.  
  161.     4 : Writeln(OutFile,
  162.                 'The maximum number of iterations must be greater than zero.');
  163.   end; { case }
  164.  
  165.   if Error <= 2 then
  166.   begin
  167.     Writeln(OutFile);
  168.     Writeln(OutFile, 'Number of iterations: ' : 26, Iter : 3);
  169.     Write(OutFile, 'Calculated root: ' : 26);
  170.     Writeln(OutFile, Answer.Re, ' + ', Answer.Im, 'i');
  171.     Writeln(OutFile, 'Value of the function  ' : 26);
  172.     Write(OutFile, 'at the calculated root: ' : 26);
  173.     Writeln(OutFile, yAnswer.Re, ' + ', yAnswer.Im, 'i');
  174.     Writeln(OutFile);
  175.   end;
  176. end; { procedure Results }
  177.  
  178. begin { program Muller }
  179.   ClrScr;
  180.   UserInput(Guess, Tol, MaxIter);
  181.   { Use Muller's method with deflation to find a root }
  182.   Muller(Guess, Tol, MaxIter, Answer, yAnswer, Iter, Error, @TNTargetF);
  183.   Results(Guess, Answer, yAnswer, Tol, MaxIter, Iter, Error);
  184.   Close(OutFile);
  185. end. { program Muller }
  186.  
  187.