home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / l / l042 / 2.ddi / CHAP8.ARC / RKF_1.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1987-12-30  |  8.8 KB  |  246 lines
  1. program RungeKuttaFehlberg_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 procedure RungeKuttaFehlberg    -}
  9. {-                   which solves an initial value problem - a first order  -}
  10. {-                   ordinary differential equation with initial conditions -}
  11. {-                   specified - to a specified degree of accuracy using    -}
  12. {-                   the Runge Kutta Fehlberg method.                       -}
  13. {-                                                                          -}
  14. {-         Unit   : InitVal    procedure RungeKuttaFehlberg                 -}
  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 X  }
  26.   XInitial : Float;                 { Intial value of X at LowerLimit  }
  27.   NumReturn : integer;              { Number of (T, X) pairs  }
  28.   Tolerance : Float;                { Tolerance in units of X/T  }
  29.   TValues : TNvector;               { Values of T  }
  30.   XValues : TNvector;               { Values of X at TValues  }
  31.   Error : byte;                     { Flags if something went wrong  }
  32.  
  33. {$F+}
  34. function TNTargetF(T, X : Float) : Float;
  35.  
  36. {---------------------------------------------------------------}
  37. {-         This is the first order differential equation       -}
  38. {---------------------------------------------------------------}
  39.  
  40. begin
  41.   TNTargetF := X / T + T - 1;
  42. end; { function TNTargetF }
  43. {$F-}
  44.  
  45. procedure Initialize(var LowerLimit : Float;
  46.                      var UpperLimit : Float;
  47.                      var XInitial   : Float;
  48.                      var NumReturn  : integer;
  49.                      var Tolerance  : Float;
  50.                      var Error      : byte);
  51.  
  52. {------------------------------------------------------------------}
  53. {- Output: LowerLimit, UpperLimit, XInitial, NumReturn,           -}
  54. {-         Tolerance, Error                                       -}
  55. {-                                                                -}
  56. {- This procedure initializes the above variables to zero         -}
  57. {------------------------------------------------------------------}
  58.  
  59. begin
  60.   LowerLimit := 0;
  61.   UpperLimit := 0;
  62.   XInitial := 0;
  63.   NumReturn := 0;
  64.   Tolerance := 0;
  65.   Error := 0;
  66. end; { procedure Initialize }
  67.  
  68. procedure GetData(var LowerLimit : Float;
  69.                   var UpperLimit : Float;
  70.                   var XInitial   : Float;
  71.                   var NumReturn  : integer;
  72.                   var Tolerance  : Float);
  73.  
  74. {------------------------------------------------------------}
  75. {- Output: LowerLimit, UpperLimit, XInitial,                -}
  76. {-         NumReturn, Tolerance                             -}
  77. {-                                                          -}
  78. {- This procedure assigns values to the above variables     -}
  79. {- from keyboard input                                      -}
  80. {------------------------------------------------------------}
  81.  
  82. procedure GetLimits(var LowerLimit : Float;
  83.                     var UpperLimit : Float);
  84.  
  85. {------------------------------------------------------------}
  86. {- Output: LowerLimit, UpperLimit                           -}
  87. {-                                                          -}
  88. {- This procedure assigns values to the limits of           -}
  89. {- integration from keyboard input                          -}
  90. {------------------------------------------------------------}
  91.  
  92. begin
  93.   repeat
  94.     repeat
  95.       Write('Lower limit of interval? ');
  96.       Readln(LowerLimit);
  97.       IOCheck;
  98.     until not IOerr;
  99.     Writeln;
  100.     repeat
  101.       Write('Upper limit of interval? ');
  102.       Readln(UpperLimit);
  103.       IOCheck;
  104.     until not IOerr;
  105.     if LowerLimit = UpperLimit then
  106.     begin
  107.       Writeln;
  108.       Writeln('       The limits of integration must be different.');
  109.       Writeln;
  110.     end;
  111.   until LowerLimit <> UpperLimit;
  112. end; { procedure GetLimits }
  113.  
  114. procedure GetXInitial(LowerLimit : Float;
  115.                   var XInitial   : Float);
  116.  
  117. {----------------------------------------------}
  118. {- Input: LowerLimit                          -}
  119. {- Output: XInitial                           -}
  120. {-                                            -}
  121. {- This procedure assigns a value to XInitial -}
  122. {- from keyboard input.                       -}
  123. {----------------------------------------------}
  124.  
  125. begin
  126.   Writeln;
  127.   repeat
  128.     Write('X value at t =', LowerLimit : 14, ': ');
  129.     Readln(XInitial);
  130.   until not IOerr;
  131. end; { procedure GetXInitial }
  132.  
  133. procedure GetTolerance(var Tolerance : Float);
  134.  
  135. {------------------------------------------------------------}
  136. {- Output: Tolerance                                        -}
  137. {-                                                          -}
  138. {- This procedure reads in the tolerance with               -}
  139. {- which to evaluate a solution.                            -}
  140. {------------------------------------------------------------}
  141.  
  142. begin
  143.   Writeln;
  144.   Tolerance := 1E-6;
  145.   repeat
  146.     Write('Tolerance in X/T units (>0)? ');
  147.     ReadFloat(Tolerance);
  148.     IOCheck;
  149.     if Tolerance <= 0 then
  150.     begin
  151.       IOerr := true;
  152.       Tolerance := 1E-6;
  153.     end;
  154.   until not IOerr;
  155. end; { procedure GetTolerance }
  156.  
  157. procedure GetNumReturn(var NumReturn : integer);
  158.  
  159. {------------------------------------------------------------}
  160. {- Output: NumReturn                                        -}
  161. {-                                                          -}
  162. {- This procedure reads in the number of values to be       -}
  163. {- returned from the procedure.                             -}
  164. {------------------------------------------------------------}
  165.  
  166. begin
  167.   Writeln;
  168.   repeat
  169.     Write('Number of values to return (1-', TNArraySize, ')? ');
  170.     Readln(NumReturn);
  171.     IOCheck;
  172.   until not IOerr and (NumReturn <= TNArraySize) and (NumReturn >= 1);
  173. end; { procedure GetNumReturn }
  174.  
  175. begin { procedure GetData }
  176.   GetLimits(LowerLimit, UpperLimit);
  177.   GetXInitial(LowerLimit, XInitial);
  178.   GetNumReturn(NumReturn);
  179.   GetTolerance(Tolerance);
  180.   GetOutputFile(OutFile);
  181. end; { procedure GetData }
  182.  
  183. procedure Results(LowerLimit : Float;
  184.                   UpperLimit : Float;
  185.                   XInitial   : Float;
  186.                   NumReturn  : integer;
  187.                   Tolerance  : Float;
  188.               var TValues    : TNvector;
  189.               var XValues    : TNvector;
  190.                   Error      : byte);
  191.  
  192. {------------------------------------------------------------}
  193. {- This procedure outputs the results to the device OutFile -}
  194. {------------------------------------------------------------}
  195.  
  196. var
  197.   Index : integer;
  198.  
  199. begin
  200.   Writeln(OutFile);
  201.   Writeln(OutFile);
  202.   Writeln(OutFile, 'Lower Limit:' : 29, LowerLimit);
  203.   Writeln(OutFile, 'Upper Limit:' : 29, UpperLimit);
  204.   Writeln(OutFile, 'Value of X at ' : 19, LowerLimit:8:4, ' :' , XInitial);
  205.   Writeln(OutFile, 'Tolerance:' : 29, Tolerance);
  206.   Writeln(OutFile);
  207.   if Error = 4 then
  208.     DisplayWarning;
  209.   if Error in [1, 2, 3] then
  210.     DisplayError;
  211.  
  212.   case Error of
  213.     1 : Writeln(OutFile, 'The tolerance must be greater than ',TNNearlyZero);
  214.  
  215.     2 : Writeln(OutFile,
  216.                 'The number of values returned must be greater than zero.');
  217.  
  218.     3 : Writeln(OutFile, 'The lower limit must be different ',
  219.                          'from the upper limit.');
  220.  
  221.     4 : begin
  222.           Writeln(OutFile, 'Tolerance not reached for all values.');
  223.           Writeln(OutFile, 'The partial results are:');
  224.           Writeln(OutFile);
  225.         end;
  226.   end; { case }
  227.  
  228.   if Error in [0, 4] then
  229.   begin
  230.     Writeln(OutFile, 't' : 15, 'X' : 15);
  231.     for Index := 0 to NumReturn do
  232.       Writeln(OutFile, TValues[Index] : 20 : 8, ' ', XValues[Index]);
  233.   end;
  234. end; { procedure Results }
  235.  
  236. begin { program RungeKuttaFehlberg }
  237.   ClrScr;
  238.   Initialize(LowerLimit, UpperLimit, XInitial, NumReturn, Tolerance, Error);
  239.   GetData(LowerLimit, UpperLimit, XInitial, NumReturn, Tolerance);
  240.   RungeKuttaFehlberg(LowerLimit, UpperLimit, XInitial, Tolerance,
  241.                      NumReturn, TValues, XValues, Error, @TNTargetF);
  242.   Results(LowerLimit, UpperLimit, XInitial, NumReturn, Tolerance,
  243.           TValues, XValues, Error);
  244.   Close(OutFile);
  245. end. { program RungeKuttaFehlberg }
  246.