home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / l / l042 / 2.ddi / CHAP8.ARC / RUNGE_S1.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1987-12-30  |  21.2 KB  |  534 lines
  1. program InitialConditionSystem_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 InitialConditionSystem. -}
  9. {-             This procedure solves a system of n first order ordinary     -}
  10. {-             differential equations with initial conditions specified     -}
  11. {-             using the fourth order  Runge-Kutta formula.                 -}
  12. {-                                                                          -}
  13. {-   Unit   : InitVal    procedure InitialConditionSystem                   -}
  14. {-                                                                          -}
  15. {----------------------------------------------------------------------------}
  16.  
  17. {$I-}      { Disable I/O error trapping  }
  18. {$R+}      { Enable range checking  }
  19.  
  20. uses
  21.   InitVal, Dos, Crt, Common;
  22.  
  23. const
  24.   TNRowSize    = TNArraySize;    { Maximum number of diff. eqs.  }
  25.   TNColumnSize = TNArraySize;    { Maximum size of vectors  }
  26.  
  27. var
  28.   NumEquations : integer;        { The number of first order equations  }
  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.   Vector : FuncVect;             { Pointers to user defined functions }
  38.  
  39. {$F+}
  40. function TNTargetF1(V : TNvector) : Float;
  41.  
  42. {--------------------------------------------------}
  43. {- This is the first differential equation.       -}
  44. {-                                                -}
  45. {-  dx[1]                                         -}
  46. {-  -----  = TNTargetF1(T, X[1], X[2], ... X[m])  -}
  47. {-    dt                                          -}
  48. {-                                                -}
  49. {- The vector V is defined:                       -}
  50. {-      V[0] = T                                  -}
  51. {-      V[1] = X[1]                               -}
  52. {-      V[2] = X[2]                               -}
  53. {-           .                                    -}
  54. {-           .                                    -}
  55. {-           .                                    -}
  56. {-      V[m] = X[m]                               -}
  57. {-                                                -}
  58. {- where m is the number of coupled equations.    -}
  59. {--------------------------------------------------}
  60.  
  61. begin
  62.   TNTargetF1 := V[2];
  63. end; { function TNTargetF1 }
  64.  
  65. function TNTargetF2(V : TNvector) : Float;
  66.  
  67. {--------------------------------------------------}
  68. {- This is the 2nd differential equation.         -}
  69. {-                                                -}
  70. {-  dx[2]                                         -}
  71. {-  -----  = TNTargetF2(T, X[1], X[2], ... X[m])  -}
  72. {-    dt                                          -}
  73. {-                                                -}
  74. {- The vector V is defined:                       -}
  75. {-      V[0] = T                                  -}
  76. {-      V[1] = X[1]                               -}
  77. {-      V[2] = X[2]                               -}
  78. {-           .                                    -}
  79. {-           .                                    -}
  80. {-           .                                    -}
  81. {-      V[m] = X[m]                               -}
  82. {-                                                -}
  83. {- where m is the number of coupled equations.    -}
  84. {--------------------------------------------------}
  85.  
  86. begin
  87.   TNTargetF2 := 9/2 * Sin(5 * V[0]) - 32/2 * V[1];
  88. end; { function TNTargetF2 }
  89.  
  90. function TNTargetF3(V : TNvector) : Float;
  91.  
  92. {--------------------------------------------------}
  93. {- This is the third differential equation.       -}
  94. {-                                                -}
  95. {-  dx[3]                                         -}
  96. {-  -----  = TNTargetF3(T, X[1], X[2], ... X[m])  -}
  97. {-    dt                                          -}
  98. {-                                                -}
  99. {- The vector V is defined:                       -}
  100. {-      V[0] = T                                  -}
  101. {-      V[1] = X[1]                               -}
  102. {-      V[2] = X[2]                               -}
  103. {-           .                                    -}
  104. {-           .                                    -}
  105. {-           .                                    -}
  106. {-      V[m] = X[m]                               -}
  107. {-                                                -}
  108. {- where m is the number of coupled equations.    -}
  109. {--------------------------------------------------}
  110.  
  111. begin
  112. end; { function TNTargetF3 }
  113.  
  114. function TNTargetF4(V : TNvector) : Float;
  115.  
  116. {--------------------------------------------------}
  117. {- This is the fourth differential equation.      -}
  118. {-                                                -}
  119. {-  dx[4]                                         -}
  120. {-  -----  = TNTargetF4(T, X[1], X[2], ... X[m])  -}
  121. {-    dt                                          -}
  122. {-                                                -}
  123. {- The vector V is defined:                       -}
  124. {-      V[0] = T                                  -}
  125. {-      V[1] = X[1]                               -}
  126. {-      V[2] = X[2]                               -}
  127. {-           .                                    -}
  128. {-           .                                    -}
  129. {-           .                                    -}
  130. {-      V[m] = X[m]                               -}
  131. {-                                                -}
  132. {- where m is the number of coupled equations.    -}
  133. {--------------------------------------------------}
  134.  
  135. begin
  136. end; { function TNTargetF4 }
  137.  
  138. function TNTargetF5(V : TNvector) : Float;
  139.  
  140. {--------------------------------------------------}
  141. {- This is the fifth differential equation.       -}
  142. {-                                                -}
  143. {-                                                -}
  144. {-  dx[5]                                         -}
  145. {-  -----  = TNTargetF5(T, X[1], X[2], ... X[m])  -}
  146. {-    dt                                          -}
  147. {-                                                -}
  148. {- The vector V is defined:                       -}
  149. {-      V[0] = T                                  -}
  150. {-      V[1] = X[1]                               -}
  151. {-      V[2] = X[2]                               -}
  152. {-           .                                    -}
  153. {-           .                                    -}
  154. {-           .                                    -}
  155. {-      V[m] = X[m]                               -}
  156. {-                                                -}
  157. {- where m is the number of coupled equations.    -}
  158. {--------------------------------------------------}
  159.  
  160. begin
  161. end; { function TNTargetF5 }
  162.  
  163. function TNTargetF6(V : TNvector) : Float;
  164.  
  165. {--------------------------------------------------}
  166. {- This is the sixth differential equation.       -}
  167. {-                                                -}
  168. {-  dx[6]                                         -}
  169. {-  -----  = TNTargetF6(T, X[1], X[2], ... X[m])  -}
  170. {-    dt                                          -}
  171. {-                                                -}
  172. {- The vector V is defined:                       -}
  173. {-      V[0] = T                                  -}
  174. {-      V[1] = X[1]                               -}
  175. {-      V[2] = X[2]                               -}
  176. {-           .                                    -}
  177. {-           .                                    -}
  178. {-           .                                    -}
  179. {-      V[m] = X[m]                               -}
  180. {-                                                -}
  181. {- where m is the number of coupled equations.    -}
  182. {--------------------------------------------------}
  183.  
  184. begin
  185. end; { function TNTargetF6 }
  186.  
  187. function TNTargetF7(V : TNvector) : Float;
  188.  
  189. {--------------------------------------------------}
  190. {- This is the seventh differential equation.     -}
  191. {-                                                -}
  192. {-                                                -}
  193. {-  dx[7]                                         -}
  194. {-  -----  = TNTargetF7(T, X[1], X[2], ... X[m])  -}
  195. {-    dt                                          -}
  196. {-                                                -}
  197. {- The vector V is defined:                       -}
  198. {-      V[0] = T                                  -}
  199. {-      V[1] = X[1]                               -}
  200. {-      V[2] = X[2]                               -}
  201. {-           .                                    -}
  202. {-           .                                    -}
  203. {-           .                                    -}
  204. {-      V[m] = X[m]                               -}
  205. {-                                                -}
  206. {- where m is the number of coupled equations.    -}
  207. {--------------------------------------------------}
  208.  
  209. begin
  210. end; { function TNTargetF7 }
  211.  
  212. function TNTargetF8(V : TNvector) : Float;
  213.  
  214. {--------------------------------------------------}
  215. {- This is the eighth differential equation.      -}
  216. {-                                                -}
  217. {-  dx[8]                                         -}
  218. {-  -----  = TNTargetF8(T, X[1], X[2], ... X[m])  -}
  219. {-    dt                                          -}
  220. {-                                                -}
  221. {- The vector V is defined:                       -}
  222. {-      V[0] = T                                  -}
  223. {-      V[1] = X[1]                               -}
  224. {-      V[2] = X[2]                               -}
  225. {-           .                                    -}
  226. {-           .                                    -}
  227. {-           .                                    -}
  228. {-      V[m] = X[m]                               -}
  229. {-                                                -}
  230. {- where m is the number of coupled equations.    -}
  231. {--------------------------------------------------}
  232.  
  233. begin
  234. end; { function TNTargetF8 }
  235.  
  236. function TNTargetF9(V : TNvector) : Float;
  237.  
  238. {--------------------------------------------------}
  239. {- This is the ninth differential equation.       -}
  240. {-                                                -}
  241. {-  dx[9]                                         -}
  242. {-  -----  = TNTargetF9(T, X[1], X[2], ... X[m])  -}
  243. {-    dt                                          -}
  244. {-                                                -}
  245. {- The vector V is defined:                       -}
  246. {-      V[0] = T                                  -}
  247. {-      V[1] = X[1]                               -}
  248. {-      V[2] = X[2]                               -}
  249. {-           .                                    -}
  250. {-           .                                    -}
  251. {-           .                                    -}
  252. {-      V[m] = X[m]                               -}
  253. {-                                                -}
  254. {- where m is the number of coupled equations.    -}
  255. {--------------------------------------------------}
  256.  
  257. begin
  258. end; { function TNTargetF9 }
  259.  
  260. function TNTargetF10(V : TNvector) : Float;
  261.  
  262. {--------------------------------------------------}
  263. {- This is the tenth differential equation.       -}
  264. {-                                                -}
  265. {-  dx[10]                                        -}
  266. {-  -----  = TNTargetF10(T, X[1], X[2], ... X[m]) -}
  267. {-    dt                                          -}
  268. {-                                                -}
  269. {- The vector V is defined:                       -}
  270. {-      V[0] = T                                  -}
  271. {-      V[1] = X[1]                               -}
  272. {-      V[2] = X[2]                               -}
  273. {-           .                                    -}
  274. {-           .                                    -}
  275. {-           .                                    -}
  276. {-      V[m] = X[m]                               -}
  277. {-                                                -}
  278. {- where m is the number of coupled equations.    -}
  279. {--------------------------------------------------}
  280.  
  281. begin
  282. end; { function TNTargetF10 }
  283. {$F-}
  284.  
  285. procedure Initialize(var NumEquations : integer;
  286.                      var LowerLimit   : Float;
  287.                      var UpperLimit   : Float;
  288.                      var InitialValue : TNvector;
  289.                      var NumIntervals : integer;
  290.                      var NumReturn    : integer;
  291.                      var Error        : byte);
  292.  
  293. {------------------------------------------------------------------}
  294. {- Output: NumEquations, LowerLimit, UpperLimit, LowerInitial,    -}
  295. {-         InitialValues, NumIntervals, NumReturn, Error          -}
  296. {-                                                                -}
  297. {- This procedure initializes the above variables to zero.        -}
  298. {------------------------------------------------------------------}
  299.  
  300. begin
  301.   NumEquations := 0;
  302.   LowerLimit := 0;
  303.   UpperLimit := 0;
  304.   FillChar(InitialValue, SizeOf(InitialValue), 0);
  305.   NumReturn := 0;
  306.   NumIntervals := 0;
  307.   Error := 0;
  308.   Vector[1] := @TNTargetF1;
  309.   Vector[2] := @TNTargetF2;
  310.   Vector[3] := @TNTargetF3;
  311.   Vector[4] := @TNTargetF4;
  312.   Vector[5] := @TNTargetF5;
  313.   Vector[6] := @TNTargetF6;
  314.   Vector[7] := @TNTargetF7;
  315.   Vector[8] := @TNTargetF8;
  316.   Vector[9] := @TNTargetF9;
  317.   Vector[10] := @TNTargetF10;
  318. end; { procedure Initialize }
  319.  
  320. procedure GetData(var NumEquations  : integer;
  321.                   var LowerLimit    : Float;
  322.                   var UpperLimit    : Float;
  323.                   var InitialValues : TNvector;
  324.                   var NumReturn     : integer;
  325.                   var NumIntervals  : integer);
  326.  
  327. {------------------------------------------------------------}
  328. {- Output: NumEquations, LowerLimit, UpperLimit,            -}
  329. {-         InitialValues, NumReturn, NumIntervals           -}
  330. {-                                                          -}
  331. {- This procedure assigns values to the above variables     -}
  332. {- from keyboard input                                      -}
  333. {------------------------------------------------------------}
  334.  
  335. procedure GetNumEquations(var NumEquations : integer);
  336.  
  337. {--------------------------------------}
  338. {- Output: NumEquations               -}
  339. {-                                    -}
  340. {- This procedure reads in the number -}
  341. {- of equations from the keyboard.    -}
  342. {--------------------------------------}
  343.  
  344. begin
  345.   Writeln;
  346.   repeat
  347.     Write('Number of first order equations: (1-', TNRowSize, ')? ');
  348.     Readln(NumEquations);
  349.     IOCheck;
  350.   until not IOerr and (NumEquations <= TNRowSize) and (NumEquations >= 1);
  351. end; { procedure GetNumEquations }
  352.  
  353. procedure GetLimits(var LowerLimit : Float;
  354.                     var UpperLimit : Float);
  355.  
  356. {------------------------------------------------------------}
  357. {- Output: LowerLimit, UpperLimit                           -}
  358. {-                                                          -}
  359. {- This procedure assigns values to the limits of           -}
  360. {- integration from keyboard input                          -}
  361. {------------------------------------------------------------}
  362.  
  363. begin
  364.   Writeln;
  365.   repeat
  366.     repeat
  367.       Write('Lower limit of interval? ');
  368.       Readln(LowerLimit);
  369.       IOCheck;
  370.     until not IOerr;
  371.     Writeln;
  372.     repeat
  373.       Write('Upper limit of interval? ');
  374.       Readln(UpperLimit);
  375.       IOCheck;
  376.     until not IOerr;
  377.     if LowerLimit = UpperLimit then
  378.     begin
  379.       Writeln;
  380.       Writeln('       The limits of integration must be different.');
  381.       Writeln;
  382.     end;
  383.   until LowerLimit <> UpperLimit;
  384. end; { procedure GetLimits }
  385.  
  386.   procedure GetInitialValues(NumEquations  : integer;
  387.                              LowerLimit    : Float;
  388.                          var InitialValues : TNvector);
  389.  
  390. {--------------------------------------------------}
  391. {- Input: NumEquations, LowerLimit                -}
  392. {- Output: InitialValues                          -}
  393. {-                                                -}
  394. {- This procedure assigns values to InitialValues -}
  395. {- from keyboard input.                           -}
  396. {--------------------------------------------------}
  397.  
  398. var
  399.   Term : integer;
  400.  
  401. begin
  402.   InitialValues[0] := LowerLimit;
  403.   Writeln;
  404.   for Term := 1 to NumEquations do
  405.   repeat
  406.     Write('Enter X[', Term, '] value at t =', LowerLimit : 14, ': ');
  407.     Readln(InitialValues[Term]);
  408.     IOCheck;
  409.   until not IOerr;
  410. end; { procedure GetInitialValues }
  411.  
  412. procedure GetNumReturn(var NumReturn : integer);
  413.  
  414. {----------------------------------------------------------}
  415. {- Output: NumReturn                                      -}
  416. {-                                                        -}
  417. {- This procedure reads in the number of values to return -}
  418. {- in the vectors TValues, XValues and XDerivValues.      -}
  419. {----------------------------------------------------------}
  420.  
  421. begin
  422.   Writeln;
  423.   repeat
  424.     Write('Number of values to return (1-', TNColumnSize, ')? ');
  425.     Readln(NumReturn);
  426.     IOCheck;
  427.   until not IOerr and (NumReturn <= TNColumnSize) and (NumReturn >= 1);
  428. end; { procedure GetNumReturn }
  429.  
  430. procedure GetNumIntervals(NumReturn    : integer;
  431.                       var NumIntervals : integer);
  432.  
  433. {------------------------------------------------------------}
  434. {- Input: NumReturn                                         -}
  435. {- Output: NumIntervals                                     -}
  436. {-                                                          -}
  437. {- This procedure reads in the number of intervals          -}
  438. {- over which to solve the equation.                        -}
  439. {------------------------------------------------------------}
  440.  
  441. begin
  442.   Writeln;
  443.   NumIntervals := NumReturn;
  444.   repeat
  445.     Write('Number of intervals (>= ', NumReturn, ')? ');
  446.     ReadInt(NumIntervals);
  447.     IOCheck;
  448.     if NumIntervals < NumReturn then
  449.     begin
  450.       IOerr := true;
  451.       NumIntervals := NumReturn;
  452.     end;
  453.   until not IOerr;
  454. end; { procedure GetNumIntervals }
  455.  
  456. begin { procedure GetData }
  457.   GetNumEquations(NumEquations);
  458.   GetLimits(LowerLimit, UpperLimit);
  459.   GetInitialValues(NumEquations, LowerLimit, InitialValues);
  460.   GetNumReturn(NumReturn);
  461.   GetNumIntervals(NumReturn, NumIntervals);
  462.   GetOutputFile(OutFile);
  463. end; { procedure GetData }
  464.  
  465. procedure Results(NumEquations   : integer;
  466.                   LowerLimit     : Float;
  467.                   UpperLimit     : Float;
  468.                   InitialValues  : TNvector;
  469.                   NumIntervals   : integer;
  470.                   NumReturn      : integer;
  471.               var SolutionValues : TNmatrix;
  472.                   Error          : byte);
  473.  
  474. {------------------------------------------------------------}
  475. {- This procedure outputs the results to the device OutFile -}
  476. {------------------------------------------------------------}
  477.  
  478. var
  479.   Term : integer;
  480.   Index : integer;
  481.  
  482. begin
  483.   Writeln(OutFile);
  484.   Writeln(OutFile);
  485.   Writeln(OutFile, 'Lower Limit:' : 35, LowerLimit);
  486.   Writeln(OutFile, 'Upper Limit:' : 35, UpperLimit);
  487.   Writeln(OutFile, 'Number of intervals:' : 35, NumIntervals);
  488.   Writeln(OutFile, 'Initial conditions at lower limit:' : 35);
  489.   for Term := 1 to NumEquations do
  490.     Writeln(OutFile, 'X[' : 21, Term,']= ', InitialValues[Term]);
  491.   Writeln(OutFile);
  492.   if Error >= 1 then
  493.     DisplayError;
  494.  
  495.   case Error of
  496.     0 : begin
  497.           for Term := 1 to NumEquations do
  498.           begin
  499.             Writeln(OutFile, 't':10, 'Value X[' : 30, Term, ']');
  500.             for Index := 0 to NumReturn do
  501.               Writeln(OutFile, SolutionValues[Index, 0] : 15 : 8,
  502.                                SolutionValues[Index, Term] : 30);
  503.             Writeln(OutFile);
  504.           end;
  505.         end;
  506.  
  507.     1 : Writeln(OutFile,
  508.                 'The number of values to return must be greater than zero.');
  509.  
  510.     2 : begin
  511.           Writeln(OutFile, 'The number of intervals must be greater than');
  512.           Writeln(OutFile, 'or equal to the number of values to return.');
  513.         end;
  514.  
  515.     3 : Writeln(OutFile,
  516.                 'The number of equations must be greater than zero.');
  517.     4 : Writeln(OutFile, 'The lower limit must be different ',
  518.                          'from the upper limit.');
  519.   end; { case }
  520. end; { procedure Results }
  521.  
  522. begin { program InitialConditionSystem }
  523.   ClrScr;
  524.   Initialize(NumEquations, LowerLimit, UpperLimit, InitialValues,
  525.              NumIntervals, NumReturn, Error);
  526.   GetData(NumEquations, LowerLimit, UpperLimit, InitialValues,
  527.           NumReturn, NumIntervals);
  528.   InitialConditionSystem(NumEquations, LowerLimit, UpperLimit, InitialValues,
  529.                          NumReturn, NumIntervals, SolutionValues, Error, Vector);
  530.   Results(NumEquations, LowerLimit, UpperLimit, InitialValues, NumIntervals,
  531.           NumReturn, SolutionValues, Error);
  532.   Close(OutFile);
  533. end. { program InitialConditionSystem }
  534.