home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / l / l042 / 2.ddi / CHAP8.ARC / SHOOT2.PAS < prev   
Encoding:
Pascal/Delphi Source File  |  1987-12-30  |  13.8 KB  |  376 lines
  1. program Shooting_Prog;
  2.  
  3. {----------------------------------------------------------------------------}
  4. {-                                                                          -}
  5. {-     Turbo Pascal Numerical Methods Toolbox                               -}
  6. {-     Copyright (c) 1986, 87 by Borland International, Inc.                -}
  7. {-                                                                          -}
  8. {-       Purpose:  This unit demonstrates the procedure Shooting            -}
  9. {-                 which solves a boundary value problem - a 2nd order,     -}
  10. {-                 non-linear ordinary differential equation with boundary  -}
  11. {-                 conditions specified - using the non-linear shooting     -}
  12. {-                 method.                                                  -}
  13. {-                                                                          -}
  14. {-       Unit   : InitVal    procedure Shooting                             -}
  15. {-                                                                          -}
  16. {----------------------------------------------------------------------------}
  17.  
  18. {$I-}      { Disable I/O error trapping  }
  19. {$R+}      { Enable range checking  }
  20.  
  21. uses
  22.   InitVal, Dos, Crt, Common;
  23.  
  24. var
  25.   LowerLimit, UpperLimit : Float;     { Limits over which to approximate Y  }
  26.   LowerInitial, UpperInitial : Float; { Boundary values of Y at limits  }
  27.   InitialSlope : Float;               { A guess of the slope at LowerLimit  }
  28.   NumReturn : integer;                { Number of points to be returned  }
  29.   NumIntervals : integer;             { Number of sub-intervals  }
  30.   Tolerance : Float;                  { Tolerance in answer  }
  31.   MaxIter : integer;                  { Maximum number of iterations  }
  32.   Iter : integer;                     { Number of iterations  }
  33.   XValues : TNvector;                 { Values of X between limits  }
  34.   YValues : TNvector;                 { Value of Y at XValues  }
  35.   YDerivValues : TNvector;            { Derivative of Y at XValues  }
  36.   Error : byte;                       { Flags if something went wrong  }
  37.  
  38. {$F+}
  39. function TNTargetF(X      : Float;
  40.                    Y      : Float;
  41.                    yPrime : Float) : Float;
  42.  
  43. {--------------------------------------------------------------------}
  44. {-    This is the 2nd order non-linear differential equation        -}
  45. {--------------------------------------------------------------------}
  46.  
  47. begin
  48.   TNTargetF := 192 * Sqr(Y / yPrime);
  49. end; { function TNTargetF }
  50. {$F-}
  51.  
  52. procedure Initialize(var LowerLimit   : Float;
  53.                      var UpperLimit   : Float;
  54.                      var LowerInitial : Float;
  55.                      var UpperInitial : Float;
  56.                      var InitialSlope : Float;
  57.                      var NumIntervals : integer;
  58.                      var NumReturn    : integer;
  59.                      var Error        : byte);
  60.  
  61. {------------------------------------------------------------------}
  62. {- Output: LowerLimit, UpperLimit, LowerInitial, UpperInitial,    -}
  63. {-         InitialSlope, NumIntervals, NumReturn, Error           -}
  64. {-                                                                -}
  65. {- This procedure initializes the above variables to zero         -}
  66. {------------------------------------------------------------------}
  67.  
  68. begin
  69.   LowerLimit := 0;
  70.   UpperLimit := 0;
  71.   LowerInitial := 0;
  72.   UpperInitial := 0;
  73.   InitialSlope := 0;
  74.   NumIntervals := 0;
  75.   NumReturn := 0;
  76.   Error := 0;
  77. end; { procedure Initialize }
  78.  
  79. procedure GetData(var LowerLimit   : Float;
  80.                   var UpperLimit   : Float;
  81.                   var LowerInitial : Float;
  82.                   var UpperInitial : Float;
  83.                   var InitialSlope : Float;
  84.                   var NumReturn    : integer;
  85.                   var NumIntervals : integer;
  86.                   var Tolerance    : Float;
  87.                   var MaxIter      : integer);
  88.  
  89. {------------------------------------------------------------}
  90. {- Output: LowerLimit, UpperLimit, LowerInitial,            -}
  91. {-         UpperInitial, NumReturn, NumIntervals,           -}
  92. {-         Tolerance, MaxIter                               -}
  93. {-                                                          -}
  94. {- This procedure assigns values to the above variables     -}
  95. {- from keyboard input                                      -}
  96. {------------------------------------------------------------}
  97.  
  98. procedure GetLimits(var LowerLimit : Float;
  99.                     var UpperLimit : Float);
  100.  
  101. {------------------------------------------------------------}
  102. {- Output: LowerLimit, UpperLimit                           -}
  103. {-                                                          -}
  104. {- This procedure assigns values to the limits of           -}
  105. {- integration from keyboard input                          -}
  106. {------------------------------------------------------------}
  107.  
  108. begin
  109.   repeat
  110.     repeat
  111.       Write('Lower limit of interval? ');
  112.       Readln(LowerLimit);
  113.       IOCheck;
  114.     until not IOerr;
  115.     Writeln;
  116.     repeat
  117.       Write('Upper limit of interval? ');
  118.       Readln(UpperLimit);
  119.       IOCheck;
  120.     until not IOerr;
  121.     if LowerLimit = UpperLimit then
  122.     begin
  123.       Writeln;
  124.       Writeln('       The limits of integration must be different.');
  125.       Writeln;
  126.     end;
  127.   until LowerLimit <> UpperLimit;
  128. end; { procedure GetLimits }
  129.  
  130. procedure GetBoundaryValues(LowerLimit   : Float;
  131.                             UpperLimit   : Float;
  132.                         var LowerInitial : Float;
  133.                         var UpperInitial : Float);
  134.  
  135. {--------------------------------------------------}
  136. {- Input: LowerLimit, UpperLimit                  -}
  137. {- Output: LowerInitial, UpperInitial             -}
  138. {-                                                -}
  139. {- This procedure assigns a value to LowerInitial -}
  140. {- and UpperInitial from keyboard input.          -}
  141. {--------------------------------------------------}
  142.  
  143. begin
  144.   Writeln;
  145.   repeat
  146.     Write('Enter Y value at X =', LowerLimit : 14, ': ');
  147.     Readln(LowerInitial);
  148.     IOCheck;
  149.   until not IOerr;
  150.   repeat
  151.     Write('Enter Y value at X =', UpperLimit : 14, ': ');
  152.     Readln(UpperInitial);
  153.     IOCheck;
  154.   until not IOerr;
  155. end; { procedure GetBoundaryValues }
  156.  
  157. procedure GetInitialSlope(LowerLimit   : Float;
  158.                           UpperLimit   : Float;
  159.                           LowerInitial : Float;
  160.                           UpperInitial : Float;
  161.                       var InitialSlope : Float);
  162.  
  163. {--------------------------------------------------}
  164. {- Input: LowerLimit, UpperLimit, LowerInitial,   -}
  165. {-        UpperInitial                            -}
  166. {- Output: InitialSlope                           -}
  167. {-                                                -}
  168. {- This procedure assigns a value to InitialSlope -}
  169. {--------------------------------------------------}
  170.  
  171. begin
  172.   Writeln;
  173.   if LowerLimit <> UpperLimit then
  174.     InitialSlope := (UpperInitial - LowerInitial) / (UpperLimit - LowerLimit);
  175.     repeat
  176.       Write('Enter a guess for the slope at X =', LowerLimit : 14, '): ');
  177.       ReadFloat(InitialSlope);
  178.       IOCheck;
  179.     until not IOerr;
  180. end; { procedure GetInitialSlope }
  181.  
  182. procedure GetNumReturn(var NumReturn : integer);
  183.  
  184. {------------------------------------------------------------}
  185. {- Output: NumReturn                                        -}
  186. {-                                                          -}
  187. {- This procedure reads in the number of points which will  -}
  188. {- be returned from the procedure.                          -}
  189. {------------------------------------------------------------}
  190.  
  191. begin
  192.   Writeln;
  193.   repeat
  194.     Write('Number of points returned (1-', TNArraySize, ')? ');
  195.     Readln(NumReturn);
  196.     IOCheck;
  197.   until not IOerr and (NumReturn <= TNArraySize) and (NumReturn >= 1);
  198. end; { procedure GetNumReturn }
  199.  
  200. procedure GetNumIntervals(NumReturn    : integer;
  201.                       var NumIntervals : integer);
  202.  
  203. {---------------------------------------------}
  204. {- Input: LowerInitial, UpperInitial,        -}
  205. {-        NumReturn                          -}
  206. {- Output: NumIntervals                      -}
  207. {-                                           -}
  208. {- This procedure reads in the NumIntervals. -}
  209. {---------------------------------------------}
  210.  
  211. begin
  212.   Writeln;
  213.   NumIntervals := NumReturn;
  214.   repeat
  215.     Write('Number of intervals (>= ', NumReturn, ')? ');
  216.     ReadInt(NumIntervals);
  217.     IOCheck;
  218.     if NumIntervals < NumReturn then
  219.     begin
  220.       IOerr := true;
  221.       NumIntervals := NumReturn;
  222.     end;
  223.   until not IOerr;
  224. end; { procedure GetNumIntervals }
  225.  
  226. procedure GetTolerance(var Tolerance : Float);
  227.  
  228. {------------------------------------------------------------}
  229. {- Output: Tolerance                                        -}
  230. {-                                                          -}
  231. {- This procedure reads in the tolerance in the answer.     -}
  232. {------------------------------------------------------------}
  233.  
  234. begin
  235.   Writeln;
  236.   Tolerance := 1E-6;
  237.   repeat
  238.     Write('Tolerance (> 0)? ');
  239.     ReadFloat(Tolerance);
  240.     IOCheck;
  241.     if Tolerance <= 0 then
  242.     begin
  243.       IOerr := true;
  244.       Tolerance := 1E-6;
  245.     end;
  246.   until not IOerr;
  247. end; { procedure GetTolerance }
  248.  
  249. procedure GetMaxIter(var MaxIter : integer);
  250.  
  251. {-------------------------------------------------------------}
  252. {- Output: MaxIter                                           -}
  253. {-                                                           -}
  254. {- This procedure reads in the maximum number of iterations. -}
  255. {-------------------------------------------------------------}
  256.  
  257. begin
  258.   Writeln;
  259.   MaxIter := 100;
  260.   repeat
  261.     Write('Maximum number of iterations (>0)? ');
  262.     ReadInt(MaxIter);
  263.     IOCheck;
  264.     if MaxIter <= 0 then
  265.     begin
  266.       IOerr := true;
  267.       MaxIter := 100;
  268.     end;
  269.   until not IOerr;
  270. end; { procedure GetMaxIter }
  271.  
  272. begin { procedure GetData }
  273.   GetLimits(LowerLimit, UpperLimit);
  274.   GetBoundaryValues(LowerLimit, UpperLimit, LowerInitial, UpperInitial);
  275.   GetInitialSlope(LowerLimit, UpperLimit, LowerInitial, UpperInitial,
  276.                   InitialSlope);
  277.   GetNumReturn(NumReturn);
  278.   GetNumIntervals(NumReturn, NumIntervals);
  279.   GetTolerance(Tolerance);
  280.   GetMaxIter(MaxIter);
  281.   GetOutputFile(OutFile);
  282. end; { procedure GetData }
  283.  
  284. procedure Results(LowerLimit   : Float;
  285.                   UpperLimit   : Float;
  286.                   LowerInitial : Float;
  287.                   UpperInitial : Float;
  288.                   InitialSlope : Float;
  289.                   NumIntervals : integer;
  290.                   NumReturn    : integer;
  291.                   Tolerance    : Float;
  292.                   MaxIter      : integer;
  293.                   Iter         : integer;
  294.               var XValues      : TNvector;
  295.               var YValues      : TNvector;
  296.               var YDerivValues : TNvector;
  297.                   Error        : byte);
  298.  
  299. {------------------------------------------------------------}
  300. {- This procedure outputs the results to the device OutFile -}
  301. {------------------------------------------------------------}
  302.  
  303. var
  304.   Index : integer;
  305.  
  306. begin
  307.   Writeln(OutFile);
  308.   Writeln(OutFile);
  309.   Writeln(OutFile, 'Lower Limit:' : 29, LowerLimit);
  310.   Writeln(OutFile, 'Upper Limit:' : 29, UpperLimit);
  311.   Writeln(OutFile, 'Value of Y at ' : 19, LowerLimit:8:4, ' :' , LowerInitial);
  312.   Writeln(OutFile, 'Value of Y at ' : 19, UpperLimit:8:4, ' :' , UpperInitial);
  313.   Writeln(OutFile, 'Initial slope at ' : 19, LowerLimit:8:4, ' :', InitialSlope);
  314.   Writeln(Outfile, 'NumIntervals:' : 29, NumIntervals : 3);
  315.   Writeln(OutFile, 'Tolerance:  ' : 31, Tolerance : 10);
  316.   Writeln(OutFile, 'Maximum number of iterations:' : 29, MaxIter);
  317.   Writeln(OutFile);
  318.   Writeln(OutFile, 'Number of iterations:' : 29, Iter:3);
  319.   Writeln(OutFile);
  320.   if Error = 6 then
  321.     DisplayWarning;
  322.   if (Error >= 1) and (Error <> 6) then
  323.     DisplayError;
  324.   case Error of
  325.     1 : Writeln(OutFile,
  326.                 'The number of points returned must be greater than zero.');
  327.  
  328.     2 : begin
  329.           Writeln(OutFile, 'The number of intervals must be greater than');
  330.           Writeln(OutFile, 'or equal to the number of values to return.');
  331.         end;
  332.  
  333.     3 : Writeln(OutFile,
  334.                 'The lower limit must not equal the upper limit.');
  335.  
  336.     4 : Writeln(OutFile, 'The tolerance must be greater than zero.');
  337.  
  338.     5 : Writeln(OutFile,
  339.                 'The maximum number of iterations must be greater than zero.');
  340.  
  341.     6 : begin
  342.           Writeln(OutFile,
  343.                   'Convergence was not reached in ', iter, ' iterations.');
  344.           Writeln(OutFile, 'Results of the last iteration:');
  345.         end;
  346.  
  347.     7 : begin
  348.           Writeln(OutFile, 'Convergence not possible.');
  349.           Writeln(OutFile,
  350.                   'Try another approximation to the initial slope.');
  351.         end;
  352.   end; { case }
  353.  
  354.   if Error in [0, 6] then
  355.   begin
  356.     Writeln(OutFile, 'X':10, 'Y Value' : 30, 'Derivative of Y' : 28);
  357.     for Index := 0 to NumReturn do
  358.       Writeln(OutFile, XValues[Index], YValues[Index] : 26, YDerivValues[Index]:24);
  359.   end;
  360. end; { procedure Results }
  361.  
  362. begin { program Shooting }
  363.   ClrScr;
  364.   Initialize(LowerLimit, UpperLimit, LowerInitial, UpperInitial,
  365.              InitialSlope, NumIntervals, NumReturn, Error);
  366.   GetData(LowerLimit, UpperLimit, LowerInitial, UpperInitial, InitialSlope,
  367.           NumReturn, NumIntervals, Tolerance, MaxIter);
  368.   Shooting(LowerLimit, UpperLimit, LowerInitial, UpperInitial, InitialSlope,
  369.            NumReturn, Tolerance, MaxIter, NumIntervals, Iter, XValues,
  370.            YValues, YDerivValues, Error, @TNTargetF);
  371.   Results(LowerLimit, UpperLimit, LowerInitial, UpperInitial, InitialSlope,
  372.           NumIntervals, NumReturn, Tolerance, MaxIter, Iter, XValues,
  373.           YValues, YDerivValues, Error);
  374.   Close(OutFile);
  375. end. { program Shooting }
  376.