home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / l / l042 / 2.ddi / CHAP8.ARC / RUNGE_N.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1987-12-30  |  11.2 KB  |  308 lines
  1. program InitialCondition_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 InitialCondition    -}
  9. {-                 which solves an initial value problem - an nth order     -}
  10. {-                 ordinary differential equation with initial conditions   -}
  11. {-                 specified - using the fourth order, generalized          -}
  12. {-                 Runge-Kutta formula.                                     -}
  13. {-                                                                          -}
  14. {-       Unit   : InitVal    procedure InitialCondition                     -}
  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. const
  25.   TNOrder = TNArraySize;
  26.  
  27. var
  28.   Order : integer;              { The order, n, of the equation  }
  29.   LowerLimit : Float;           { Lower limit of interval  }
  30.   UpperLimit : Float;           { Upper limit of interval  }
  31.   InitialValues : TNvector;     { Initial values at lower limit  }
  32.   NumReturn : integer;          { Number of values to return  }
  33.   NumIntervals : integer;       { Number of intervals  }
  34.   SolutionValues : TNmatrix;    { Values of the n parameters at NumReturn  }
  35.                                 { Mesh points between the intervals  }
  36.   Error : byte;                 { Flags if something went wrong  }
  37.  
  38. {$F+}
  39. function TNTargetF(V : TNvector) : Float;
  40.  
  41. {-------------------------------------------}
  42. {- This is the differential equation.      -}
  43. {-   n                                     -}
  44. {-  d X                               n-1  -}
  45. {-  -----  = TNTargetF(t, X, X', ... X   ) -}
  46. {-     n                                   -}
  47. {-   dt                                    -}
  48. {-                                         -}
  49. {- where n is the order of the equation.   -}
  50. {-                                         -}
  51. {- V[0] corresponds to t                   -}
  52. {- V[1] corresponds to X                   -}
  53. {- V[2] corresponds to 1st derivative of X -}
  54. {- V[3] corresponds to 2nd derivative of X -}
  55. {-                  .                      -}
  56. {-                  .                      -}
  57. {-                  .                      -}
  58. {-------------------------------------------}
  59.  
  60. begin
  61.   TNTargetF := -4 * V[1] * V[4];
  62. end; { function TNTargetF }
  63. {$F-}
  64.  
  65. procedure Initialize(var Order        : integer;
  66.                      var LowerLimit   : Float;
  67.                      var UpperLimit   : Float;
  68.                      var InitialValue : TNvector;
  69.                      var NumIntervals : integer;
  70.                      var NumReturn    : integer;
  71.                      var Error        : byte);
  72.  
  73. {------------------------------------------------------------------}
  74. {- Output: Order, LowerLimit, UpperLimit, LowerInitial,           -}
  75. {-         InitialValues, NumIntervals, NumReturn, Error          -}
  76. {-                                                                -}
  77. {- This procedure initializes the above variables to zero.        -}
  78. {------------------------------------------------------------------}
  79.  
  80. begin
  81.   Order := 0;
  82.   LowerLimit := 0;
  83.   UpperLimit := 0;
  84.   FillChar(InitialValue, SizeOf(InitialValue), 0);
  85.   NumReturn := 0;
  86.   NumIntervals := 0;
  87.   Error := 0;
  88. end; { procedure Initialize }
  89.  
  90. procedure GetData(var Order         : integer;
  91.                   var LowerLimit    : Float;
  92.                   var UpperLimit    : Float;
  93.                   var InitialValues : TNvector;
  94.                   var NumReturn     : integer;
  95.                   var NumIntervals  : integer);
  96.  
  97. {------------------------------------------------------------}
  98. {- Output: Order, LowerLimit, UpperLimit, InitialValues,    -}
  99. {-         NumReturn, NumIntervals                          -}
  100. {-                                                          -}
  101. {- This procedure assigns values to the above variables     -}
  102. {- from keyboard input                                      -}
  103. {------------------------------------------------------------}
  104.  
  105. procedure GetOrder(var Order : integer);
  106.  
  107. {--------------------------------------}
  108. {- Output: Order                      -}
  109. {-                                    -}
  110. {- This procedure reads in the order  -}
  111. {- of the equation from the keyboard. -}
  112. {--------------------------------------}
  113.  
  114. begin
  115.   Writeln;
  116.   repeat
  117.     Write('Order of the equation: (1-', TNorder, ')? ');
  118.     Readln(Order);
  119.     IOCheck;
  120.   until not IOerr and (Order <= TNorder) and (Order >= 1);
  121.   Writeln;
  122. end; { procedure GetOrder }
  123.  
  124. procedure GetLimits(var LowerLimit : Float;
  125.                     var UpperLimit : Float);
  126.  
  127. {------------------------------------------------------------}
  128. {- Output: LowerLimit, UpperLimit                           -}
  129. {-                                                          -}
  130. {- This procedure assigns values to the limits of           -}
  131. {- integration from keyboard input                          -}
  132. {------------------------------------------------------------}
  133.  
  134. begin
  135.   repeat
  136.     repeat
  137.       Write('Lower limit of interval? ');
  138.       Readln(LowerLimit);
  139.       IOCheck;
  140.     until not IOerr;
  141.     Writeln;
  142.     repeat
  143.       Write('Upper limit of interval? ');
  144.       Readln(UpperLimit);
  145.       IOCheck;
  146.     until not IOerr;
  147.     if LowerLimit = UpperLimit then
  148.     begin
  149.       Writeln;
  150.       Writeln('       The limits of integration must be different.');
  151.       Writeln;
  152.     end;
  153.   until LowerLimit <> UpperLimit;
  154. end; { procedure GetLimits }
  155.  
  156. procedure GetInitialValues(Order         : integer;
  157.                            LowerLimit    : Float;
  158.                        var InitialValues : TNvector);
  159.  
  160. {--------------------------------------------------}
  161. {- Input: Order, LowerLimit                       -}
  162. {- Output: InitialValues                          -}
  163. {-                                                -}
  164. {- This procedure assigns values to InitialValues -}
  165. {- from keyboard input.                           -}
  166. {--------------------------------------------------}
  167.  
  168. var
  169.   Term : integer;
  170.  
  171. begin
  172.   InitialValues[0] := LowerLimit;
  173.   Writeln;
  174.   repeat
  175.     Write('Enter X value at t =':24, LowerLimit : 14, ': ');
  176.     Readln(InitialValues[1]);
  177.     IOCheck;
  178.   until not IOerr;
  179.   for Term := 2 to Order do
  180.   repeat
  181.     Write('Derivative ', Term - 1, ' of X at t =', LowerLimit:14,': ');
  182.     Readln(InitialValues[Term]);
  183.     IOCheck;
  184.   until not IOerr;
  185. end; { procedure GetInitialValues }
  186.  
  187. procedure GetNumReturn(var NumReturn : integer);
  188.  
  189. {----------------------------------------------------------}
  190. {- Output: NumReturn                                      -}
  191. {-                                                        -}
  192. {- This procedure reads in the number of values to return -}
  193. {- in the vectors TValues, XValues and XDerivValues.      -}
  194. {----------------------------------------------------------}
  195.  
  196. begin
  197.   Writeln;
  198.   repeat
  199.     Write('Number of values to return (1-', TNArraySize, ')? ');
  200.     Readln(NumReturn);
  201.     IOCheck;
  202.   until not IOerr and (NumReturn <= TNArraySize) and (NumReturn >= 1);
  203. end; { procedure GetNumReturn }
  204.  
  205. procedure GetNumIntervals(NumReturn    : integer;
  206.                       var NumIntervals : integer);
  207.  
  208. {------------------------------------------------------------}
  209. {- Input: NumReturn                                         -}
  210. {- Output: NumIntervals                                     -}
  211. {-                                                          -}
  212. {- This procedure reads in the number of intervals          -}
  213. {- over which to solve the equation.                        -}
  214. {------------------------------------------------------------}
  215.  
  216. begin
  217.   Writeln;
  218.   NumIntervals := NumReturn;
  219.   repeat
  220.     Write('Number of intervals (>= ', NumReturn, ')? ');
  221.     ReadInt(NumIntervals);
  222.     IOCheck;
  223.     if NumIntervals < NumReturn then
  224.     begin
  225.       IOerr := true;
  226.       NumIntervals := NumReturn;
  227.     end;
  228.   until not IOerr;
  229. end; { procedure GetNumIntervals }
  230.  
  231. begin { procedure GetData }
  232.   GetOrder(Order);
  233.   GetLimits(LowerLimit, UpperLimit);
  234.   GetInitialValues(Order, LowerLimit, InitialValues);
  235.   GetNumReturn(NumReturn);
  236.   GetNumIntervals(NumReturn, NumIntervals);
  237.   GetOutputFile(OutFile);
  238. end; { procedure GetData }
  239.  
  240. procedure Results(Order          : integer;
  241.                   LowerLimit     : Float;
  242.                   UpperLimit     : Float;
  243.                   InitialValues  : TNvector;
  244.                   NumIntervals   : integer;
  245.                   NumReturn      : integer;
  246.               var SolutionValues : TNmatrix;
  247.                   Error          : byte);
  248.  
  249. {------------------------------------------------------------}
  250. {- This procedure outputs the results to the device OutFile -}
  251. {------------------------------------------------------------}
  252.  
  253. var
  254.   Term : integer;
  255.   Index : integer;
  256.  
  257. begin
  258.   Writeln(OutFile);
  259.   Writeln(OutFile);
  260.   Writeln(OutFile, 'Lower Limit:' : 35, LowerLimit);
  261.   Writeln(OutFile, 'Upper Limit:' : 35, UpperLimit);
  262.   Writeln(OutFile, 'Number of intervals:' : 35, NumIntervals);
  263.   Writeln(OutFile, 'Initial conditions at lower limit:' : 35);
  264.   for Term := 1 to Order do
  265.     Writeln(OutFile, 'X[' : 21, Term,']=', InitialValues[Term]);
  266.   Writeln(OutFile);
  267.   if Error >= 1 then
  268.     DisplayError;
  269.  
  270.   case Error of
  271.     0 : begin
  272.           for Term := 1 to Order do
  273.           begin
  274.             Writeln(OutFile, 't':10, 'Value X[' : 30, Term, ']');
  275.             for Index := 0 to NumReturn do
  276.               Writeln(OutFile, SolutionValues[Index, 0] : 15 : 8,
  277.                                SolutionValues[Index, Term] : 30);
  278.             Writeln(OutFile);
  279.           end;
  280.         end;
  281.  
  282.     1 : Writeln(OutFile,
  283.                 'The number of values to return must be greater than zero.');
  284.     2 : begin
  285.           Writeln(OutFile, 'The number of intervals must be greater than');
  286.           Writeln(OutFile, 'or equal to the number of values to return.');
  287.         end;
  288.  
  289.     3 : Writeln(OutFile,
  290.                 'The order of the equation must be greater than zero.');
  291.     4 : Writeln(OutFile, 'The lower limit must be different ',
  292.                          'from the upper limit.');
  293.   end; { case }
  294. end; { procedure Results }
  295.  
  296. begin { program InitialCondition }
  297.   ClrScr;
  298.   Initialize(Order, LowerLimit, UpperLimit, InitialValues,
  299.              NumIntervals, NumReturn, Error);
  300.   GetData(Order, LowerLimit, UpperLimit, InitialValues,
  301.           NumReturn, NumIntervals);
  302.   InitialCondition(Order, LowerLimit, UpperLimit, InitialValues,
  303.                    NumReturn, NumIntervals, SolutionValues, Error, @TNTargetF);
  304.   Results(Order, LowerLimit, UpperLimit, InitialValues, NumIntervals,
  305.           NumReturn, SolutionValues, Error);
  306.   Close(OutFile);
  307. end. { program InitialCondition }
  308.