home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / l / l042 / 2.ddi / CHAP8.ARC / INITVAL1.INC < prev    next >
Encoding:
Text File  |  1987-12-30  |  37.1 KB  |  1,077 lines
  1. {----------------------------------------------------------------------------}
  2. {-                                                                          -}
  3. {-     Turbo Pascal Numerical Methods Toolbox                               -}
  4. {-     Copyright (c) 1986, 87 by Borland International, Inc.                -}
  5. {-                                                                          -}
  6. {----------------------------------------------------------------------------}
  7.  
  8. procedure InitialCond1stOrder{(LowerLimit   : Float;
  9.                               UpperLimit    : Float;
  10.                               XInitial      : Float;
  11.                               NumReturn     : integer;
  12.                               NumIntervals  : integer;
  13.                           var TValues       : TNvector;
  14.                           var XValues       : TNvector;
  15.                           var Error         : byte;
  16.                               FuncPtr       : Pointer)};
  17.  
  18. type
  19.   Ptr = ^Data;
  20.  
  21.   Data = record
  22.            T, X : Float;
  23.            Next : Ptr;
  24.          end;
  25.  
  26.   Stack = Ptr;
  27.  
  28. var
  29.   Values : Stack;
  30.   Spacing, HalfSpacing : Float;        { Size of each subinterval  }
  31.   T, X, F1, F2, F3, F4 : Float;        { Step values  }
  32.   Index : integer;                    { A counter  }
  33.  
  34. procedure InitializeStack(var Values : Stack);
  35. begin
  36.   Values := nil;
  37. end; { procedure InitializeStack }
  38.  
  39. procedure Push(var Values : Stack;
  40.                    TValue : Float;
  41.                    XValue : Float);
  42.  
  43. {----------------------------------}
  44. {- Input: Values, TValue, XValue  -}
  45. {- Output: Values                 -}
  46. {-                                -}
  47. {- Push data onto the Stack       -}
  48. {----------------------------------}
  49.  
  50. var
  51.   NewNode : Ptr;
  52.  
  53. begin
  54.   New(NewNode);
  55.   NewNode^.T := TValue;
  56.   NewNode^.X := XValue;
  57.   NewNode^.Next := Values;
  58.   Values := NewNode;
  59. end; { procedure Push }
  60.  
  61. procedure Pop(var Values : Stack;
  62.               var TValue : Float;
  63.               var XValue : Float);
  64.  
  65. {----------------------------------}
  66. {- Input: Values, TValue, XValue  -}
  67. {- Output: Values                 -}
  68. {-                                -}
  69. {- Pop Data from the Stack        -}
  70. {----------------------------------}
  71.  
  72. var
  73.   OldNode : Ptr;
  74.  
  75. begin
  76.   OldNode := Values;
  77.   Values := OldNode^.Next;
  78.   TValue := OldNode^.T;
  79.   XValue := OldNode^.X;
  80.   Dispose(OldNode);
  81. end; { procedure Pop }
  82.  
  83. procedure TestAndInitialize(LowerLimit   : Float;
  84.                             UpperLimit   : Float;
  85.                             XInitial     : Float;
  86.                         var NumIntervals : integer;
  87.                         var NumReturn    : integer;
  88.                         var Values       : Stack;
  89.                         var Spacing      : Float;
  90.                         var HalfSpacing  : Float;
  91.                         var T            : Float;
  92.                         var X            : Float;
  93.                         var Error        : byte);
  94.  
  95. {---------------------------------------------------------------}
  96. {- Input: LowerLimit, UpperLimit, XInitial, NumIntervals,      -}
  97. {-        NumReturn                                            -}
  98. {- Output: Values, Spacing, HalfSpacing, T, X, Error           -}
  99. {-                                                             -}
  100. {- This procedure initializes the above variables.  Values is  -}
  101. {- initialized to NIL;  Spacing is initialized to              -}
  102. {- (UpperLimit - LowerLimit)/NumIntervals, HalfSpacing is half -}
  103. {- this; T is initialized to LowerLimit; X is initialized to   -}
  104. {- XInitial.  Also NumIntervals and NumReturn are checked for  -}
  105. {- errors (it must be greater than zero).                      -}
  106. {---------------------------------------------------------------}
  107.  
  108. begin
  109.   Error := 0;
  110.   if NumReturn < 1 then
  111.     Error := 1;
  112.   if NumIntervals < NumReturn then
  113.     Error := 2;
  114.   if LowerLimit = UpperLimit then
  115.     Error := 3;
  116.   if Error = 0 then
  117.   begin
  118.     InitializeStack(Values);
  119.     Spacing := (UpperLimit - LowerLimit) / NumIntervals;
  120.     HalfSpacing := Spacing / 2;
  121.     T := LowerLimit;
  122.     X := XInitial;
  123.     Push(Values, T, X);
  124.   end;
  125. end; { procedure TestAndInitialize }
  126.  
  127. procedure GetValues(NumReturn    : integer;
  128.                     NumIntervals : integer;
  129.                 var Values       : Stack;
  130.                 var TValues      : TNvector;
  131.                 var XValues      : TNvector);
  132.  
  133. {-------------------------------------------------------------}
  134. {- Input: NumReturn, NumIntervals, Values                    -}
  135. {- Output: TValues, XValues                                  -}
  136. {-                                                           -}
  137. {- This procedure fills in the arrays TValues and XValues    -}
  138. {- with data from the Values stack. For example, if there    -}
  139. {- are five times as many items on the stack (NumIntervals)  -}
  140. {- as should be returned (NumReturn) then every fifth value  -}
  141. {- from the stack will be returned.                          -}
  142. {-------------------------------------------------------------}
  143.  
  144. var
  145.   Index, Term : integer;
  146.   T, X : Float;
  147.  
  148. begin
  149.   Term := NumIntervals;
  150.   for Index := NumReturn downto 0 do
  151.   begin
  152.     Pop(Values, T, X);
  153.     Term := Pred(Term);
  154.     while (Term / NumIntervals >= Index / NumReturn) and (Term >= 0) do
  155.     begin
  156.       Pop(Values, T, X);
  157.       Term := Pred(Term);
  158.     end;
  159.     TValues[Index] := T;
  160.     XValues[Index] := X;
  161.   end;
  162. end; { procedure GetValues }
  163.  
  164. begin  { procedure InitialCond1stOrder }
  165.   TestAndInitialize(LowerLimit, UpperLimit, XInitial, NumIntervals,
  166.                     NumReturn, Values, Spacing, HalfSpacing, T, X, Error);
  167.   if Error = 0 then
  168.   begin
  169.     for Index := 1 to NumIntervals do
  170.     begin
  171.       F1 := Spacing * UserFunction1(T, X, FuncPtr);
  172.       F2 := Spacing * UserFunction1(T + HalfSpacing, X + F1 / 2, FuncPtr);
  173.       F3 := Spacing * UserFunction1(T + HalfSpacing, X + F2 / 2, FuncPtr);
  174.       F4 := Spacing * UserFunction1(T + Spacing, X + F3, FuncPtr);
  175.       X := X + (F1 + 2 * F2 + 2 * F3 + F4) / 6;
  176.       T := T + Spacing;
  177.       Push(Values, T, X);
  178.     end;
  179.     GetValues(NumReturn, NumIntervals, Values, TValues, XValues);
  180.   end;
  181. end; { procedure InitialCond1stOrder }
  182.  
  183. procedure RungeKuttaFehlberg{(LowerLimit : Float;
  184.                               UpperLimit : Float;
  185.                               XInitial   : Float;
  186.                               Tolerance  : Float;
  187.                               NumReturn  : integer;
  188.                           var TValues    : TNvector;
  189.                           var XValues    : TNvector;
  190.                           var Error      : byte;
  191.                               FuncPtr    : Pointer)};
  192.  
  193. type
  194.   Ptr = ^Data;
  195.  
  196.   Data = record
  197.            T, X : Float;
  198.            Next : Ptr;
  199.          end;
  200.  
  201.   Stack = Ptr;
  202.  
  203. var
  204.   NumValues : integer;        { Number of values on the stack  }
  205.   Values : Stack;             { Pointer to stack of all computed values  }
  206.   MaxSpacing : Float;          { Maximum size of each subinterval  }
  207.   Spacing : Float;             { Size of each subinterval  }
  208.   T, X, TempX : Float;         { Step values  }
  209.   Difference : Float;          { Fractional difference between steps  }
  210.  
  211. procedure InitializeStack(var Values : Stack);
  212. begin
  213.   Values := nil;
  214. end; { procedure InitializeStack }
  215.  
  216. procedure Push(var Values : Stack;
  217.                    TValue : Float;
  218.                    XValue : Float);
  219.  
  220. {----------------------------------}
  221. {- Input: Values, TValue, XValues -}
  222. {- Output: Values                 -}
  223. {-                                -}
  224. {- Push data onto the Stack       -}
  225. {----------------------------------}
  226.  
  227. var
  228.   NewNode : Ptr;
  229.  
  230. begin
  231.   New(NewNode);
  232.   NewNode^.T := TValue;
  233.   NewNode^.X := XValue;
  234.   NewNode^.Next := Values;
  235.   Values := NewNode;
  236. end; { procedure Push }
  237.  
  238. procedure Pop(var Values : Stack;
  239.               var TValue : Float;
  240.               var XValue : Float);
  241.  
  242. {----------------------------------}
  243. {- Input: Values, TValue, XValues -}
  244. {- Output: Values                 -}
  245. {-                                -}
  246. {- Pop Data from the Stack        -}
  247. {----------------------------------}
  248.  
  249. var
  250.   OldNode : Ptr;
  251.  
  252. begin
  253.   OldNode := Values;
  254.   Values := OldNode^.Next;
  255.   TValue := OldNode^.T;
  256.   XValue := OldNode^.X;
  257.   Dispose(OldNode);
  258. end; { procedure Pop }
  259.  
  260.  
  261. procedure TestAndInitialize(LowerLimit : Float;
  262.                             UpperLimit : Float;
  263.                             XInitial   : Float;
  264.                             Tolerance  : Float;
  265.                             NumReturn  : integer;
  266.                         var NumValues  : integer;
  267.                         var TValues    : TNvector;
  268.                         var XValues    : TNvector;
  269.                         var Values     : Stack;
  270.                         var MaxSpacing : Float;
  271.                         var Spacing    : Float;
  272.                         var T          : Float;
  273.                         var X          : Float;
  274.                         var TempX      : Float;
  275.                         var Error      : byte);
  276.  
  277. {------------------------------------------------------------------------}
  278. {- Input: LowerLimit, UpperLimit, XInitial, Tolerance, NumReturn        -}
  279. {- Output: NumValues, TValues, XValues, Values, MaxSpacing, Spacing,    -}
  280. {-         T, X, TempX, Error                                           -}
  281. {-                                                                      -}
  282. {- This procedure initializes the above variables.  NumValues, TValues, -}
  283. {- XValues and are initialized to zero, Values is initialized to NIL,   -}
  284. {- X and TempX are initialized to XInitial, T is initialized to         -}
  285. {- LowerLimit; MaxSpacing and  Spacing are initialized to               -}
  286. {- (LowerLimit - UpperLimit)/NumReturn.  Also, Tolerance and NumReturn  -}
  287. {- check for errors (they must be greater than TNNearlyZero).           -}
  288. {------------------------------------------------------------------------}
  289.  
  290. begin
  291.   Error := 0;
  292.   if Tolerance <= TNNearlyZero then
  293.     Error := 1;      { Tolerance less than zero }
  294.   if NumReturn < 1 then
  295.     Error := 2;      { NumReturn less than one  }
  296.   if LowerLimit = UpperLimit then
  297.     Error := 3;      { End points are identical  }
  298.   NumValues := 0;
  299.   FillChar(XValues, SizeOf(XValues), 0);
  300.   TValues[0] := LowerLimit;
  301.   XValues[0] := XInitial;
  302.   InitializeStack(Values);
  303.   T := LowerLimit;
  304.   X := XInitial;
  305.   TempX := XInitial;
  306.   if Error = 0 then
  307.   begin
  308.     MaxSpacing := (UpperLimit - LowerLimit) / NumReturn;
  309.     Spacing := MaxSpacing;
  310.   end;
  311. end; { procedure TestAndInitialize }
  312.  
  313. procedure RungeKuttaSix(Spacing    : Float;
  314.                         T          : Float;
  315.                     var X          : Float;
  316.                     var Difference : Float);
  317.  
  318. {-------------------------------------------------------------------}
  319. {- Input: Spacing, T                                               -}
  320. {- Output: X, Difference                                           -}
  321. {-                                                                 -}
  322. {- This procedure applies the six-stage, fifth-order Runge-Kutta   -}
  323. {- formula and the four-stage, fourth-order Runge-Kutta formula to -}
  324. {- approximate the value of X at T + Spacing given the value of X  -}
  325. {- at T.  The difference between the fourth order and fifth order  -}
  326. {- formulas is returned as Difference.                             -}
  327. {-------------------------------------------------------------------}
  328.  
  329. var
  330.   F1, F2, F3, F4, F5, F6 : Float;
  331.   DummyT, DummyX : Float;
  332.  
  333. begin
  334.   F1 := Spacing * UserFunction1(T, X, FuncPtr);
  335.   DummyT := T + Spacing / 4;
  336.   DummyX := X + F1 / 4;
  337.   F2 := Spacing * UserFunction1(DummyT, DummyX, FuncPtr);
  338.   DummyT := T + 3 * Spacing / 8;
  339.   DummyX := X + (3 * F1 + 9 * F2) / 32;
  340.   F3 := Spacing * UserFunction1(DummyT, DummyX, FuncPtr);
  341.   DummyT := T + 12 * Spacing / 13;
  342.   DummyX := X + (1932 * F1 - 7200 * F2 + 7296 * F3) / 2197;
  343.   F4 := Spacing * UserFunction1(DummyT, DummyX, FuncPtr);
  344.   DummyT := T + Spacing;
  345.   DummyX := X + (8341 * F1 - 32832.0 * F2 + 29440 * F3 - 845 * F4) / 4104;
  346.   F5 := Spacing * UserFunction1(DummyT, DummyX, FuncPtr);
  347.   DummyT := T + Spacing / 2;
  348.   DummyX := X + (-6080 * F1 + 41040.0 * F2 - 28352 * F3 +
  349.                   9295 * F4 - 5643 * F5) / 20520;
  350.   F6 := Spacing * UserFunction1(DummyT, DummyX, FuncPtr);
  351.   X := X + (2375 * F1 + 11264 * F3 + 10985 * F4 - 4104 * F5) / 20520;
  352.   Difference := ABS(3135.0 * F1 - 33792.0 * F3 - 32955.0 * F4 +
  353.                 22572.0 * F5 + 41040.0 * F6) / 1128600.0;
  354. end; { procedure RungeKuttaSix }
  355.  
  356. procedure CalculateNewSpacing(Tolerance  : Float;
  357.                               Difference : Float;
  358.                               MaxSpacing : Float;
  359.                           var Spacing    : Float);
  360.  
  361. {----------------------------------------------------------------}
  362. {- Input: Tolerance, Difference, MaxSpacing                     -}
  363. {- Output: Spacing                                              -}
  364. {-                                                              -}
  365. {- If the Difference between to two Runge-Kutta evaluations     -}
  366. {- is larger or smaller than the Tolerance then reduce or       -}
  367. {- increase the size of the Spacing appropriately.  However,    -}
  368. {- the Spacing should not be reduced smaller than 0.1 - Spacing -}
  369. {- nor should it be increased larger than MaxSpacing.           -}
  370. {----------------------------------------------------------------}
  371.  
  372. var
  373.   DeltaDif : Float;
  374.  
  375. begin
  376.   if ABS(Difference) > TNNearlyZero then
  377.   begin
  378.     DeltaDif := Sqrt(Sqrt(ABS(Spacing * Tolerance/(2 * Difference))));
  379.     if DeltaDif < 0.1 then
  380.       Spacing := 0.1 * Spacing
  381.     else
  382.       begin
  383.         Spacing := DeltaDif * Spacing;
  384.         if ABS(Spacing) > ABS(MaxSpacing) then
  385.           Spacing := MaxSpacing;
  386.       end
  387.   end;
  388. end; { procedure CalculateNewSpacing }
  389.  
  390. procedure GetValues(NumReturn : integer;
  391.                 var NumValues : integer;
  392.                 var Values    : Stack;
  393.                 var TValues   : TNvector;
  394.                 var XValues   : TNvector);
  395.  
  396. {-------------------------------------------------------------}
  397. {- Input: NumReturn, NumValues                               -}
  398. {- Output: NumValues, Values, TValues, XValues               -}
  399. {-                                                           -}
  400. {- This procedure fills in the arrays TValues and XValues    -}
  401. {- with data from the Values stack. If NumReturn (the        -}
  402. {- number of values to be returned) is less than             -}
  403. {- NumValues (the number of items on the stack) then only    -}
  404. {- some of the values on the stack will be returned.  For    -}
  405. {- example, if there are five times as many items on the     -}
  406. {- stack (NumValues) as should be returned (NumReturn) then  -}
  407. {- every fifth value from the stack will be returned.        -}
  408. {-------------------------------------------------------------}
  409.  
  410. var
  411.   Index, Term : integer;
  412.   TVal, XVal : Float;
  413.  
  414. begin
  415.   if NumValues > NumReturn then
  416.     begin
  417.       Term := NumValues;
  418.       for Index := NumReturn downto 1 do
  419.       begin
  420.         Pop(Values, TVal, XVal);
  421.         Term := Pred(Term);
  422.         while Term/NumValues > Index/NumReturn do
  423.         begin
  424.           Pop(Values, TVal, XVal);
  425.           Term := Pred(Term);
  426.         end;
  427.         TValues[Index] := TVal;
  428.         XValues[Index] := XVal;
  429.       end;
  430.     end
  431.   else
  432.     for Index := NumValues downto 1 do
  433.       Pop(Values, TValues[Index], XValues[Index]);
  434. end; { procedure GetValues }
  435.  
  436. begin { procedure RungeKuttaFehlberg }
  437.   TestAndInitialize(LowerLimit, UpperLimit, XInitial, Tolerance, NumReturn,
  438.                     NumValues, TValues, XValues, Values,
  439.                     MaxSpacing, Spacing, T, X, TempX, Error);
  440.   if Error = 0 then
  441.   begin
  442.     while (ABS((T - LowerLimit)) < ABS((UpperLimit - LowerLimit))) and
  443.           (ABS(Spacing) > TNNearlyZero) do
  444.     begin
  445.       RungeKuttaSix(Spacing, T, TempX, Difference);
  446.       if Difference > ABS(Tolerance * Spacing) then
  447.         TempX := X
  448.       else
  449.         begin
  450.           T := T + Spacing;
  451.           X := TempX;
  452.           Push(Values, T, X);
  453.           NumValues := Succ(NumValues);
  454.         end;
  455.       CalculateNewSpacing(Tolerance, Difference, MaxSpacing, Spacing);
  456.     end; { while }
  457.  
  458.     if (ABS(Spacing) <= TNNearlyZero) and (T < UpperLimit) then
  459.       Error := 4;     { Tolerance not reached  }
  460.     GetValues(NumReturn, NumValues, Values, TValues, XValues);
  461.   end;
  462. end; { procedure RungeKuttaFehlberg }
  463.  
  464. procedure Adams{(LowerLimit  : Float;
  465.                 UpperLimit   : Float;
  466.                 XInitial     : Float;
  467.                 NumReturn    : integer;
  468.                 NumIntervals : integer;
  469.             var TValues      : TNvector;
  470.             var XValues      : TNvector;
  471.             var Error        : byte;
  472.                 FuncPtr      : Pointer)};
  473.  
  474. type
  475.   TNFourVector = array[1..4] of Float;
  476.  
  477.   Ptr = ^Data;
  478.  
  479.   Data = record
  480.            T, X : Float;
  481.            Next : Ptr;
  482.          end;
  483.  
  484.   Stack = Ptr;
  485.  
  486. var
  487.   Values : Stack;
  488.   Spacing : Float;                 { Size of each subinterval }
  489.   T : Float;                       { Value of T }
  490.   X, NewX : Float;                 { Iteration variables }
  491.   Index : integer;                { A counter }
  492.   DerivValues : TNFourVector;     { Derivative at last 4 X values }
  493.  
  494. procedure InitializeStack(var Values : Stack);
  495. begin
  496.   Values := nil;
  497. end; { procedure InitializeStack }
  498.  
  499. procedure Push(var Values : Stack;
  500.                    TValue : Float;
  501.                    XValue : Float);
  502.  
  503. {----------------------------------}
  504. {- Input: Values, TValue, XValue  -}
  505. {- Output: Values                 -}
  506. {-                                -}
  507. {- Push data onto the Stack       -}
  508. {----------------------------------}
  509.  
  510. var
  511.   NewNode : Ptr;
  512.  
  513. begin
  514.   New(NewNode);
  515.   NewNode^.T := TValue;
  516.   NewNode^.X := XValue;
  517.   NewNode^.Next := Values;
  518.   Values := NewNode;
  519. end; { procedure Push }
  520.  
  521. procedure Pop(var Values : Stack;
  522.               var TValue : Float;
  523.               var XValue : Float);
  524.  
  525. {----------------------------------}
  526. {- Input: Values, TValue, XValue  -}
  527. {- Output: Values                 -}
  528. {-                                -}
  529. {- Pop Data from the Stack        -}
  530. {----------------------------------}
  531.  
  532. var
  533.   OldNode : Ptr;
  534.  
  535. begin
  536.   OldNode := Values;
  537.   Values := OldNode^.Next;
  538.   TValue := OldNode^.T;
  539.   XValue := OldNode^.X;
  540.   Dispose(OldNode);
  541. end; { procedure Pop }
  542.  
  543. procedure TestAndInitialize(LowerLimit   : Float;
  544.                             UpperLimit   : Float;
  545.                             XInitial     : Float;
  546.                         var NumIntervals : integer;
  547.                             NumReturn    : integer;
  548.                         var Values       : Stack;
  549.                         var DerivValues  : TNFourVector;
  550.                         var Spacing      : Float;
  551.                         var T            : Float;
  552.                         var X            : Float;
  553.                         var Error        : byte);
  554.  
  555. {--------------------------------------------------------------------}
  556. {- Input: LowerLimit, UpperLimit, XInitial, NumIntervals            -}
  557. {-        NumReturn                                                 -}
  558. {- Output: Values, DerivValues, Spacing, T, X, Error                -}
  559. {-                                                                  -}
  560. {- This procedure initializes the above variables. Values is set to -}
  561. {- Nil; DerivValues is initialized to zero;  Spacing is initialized -}
  562. {- to (UpperLimit - LowerLimit)/NumIntervals; T is initialized to   -}
  563. {- LowerLimit; NumIntervals and NumReturn are checked for errors    -}
  564. {- (it must be greater than zero).                                  -}
  565. {--------------------------------------------------------------------}
  566.  
  567. begin
  568.   Error := 0;
  569.   if NumReturn <= 0 then
  570.     Error := 1;
  571.   if NumIntervals < NumReturn then
  572.     Error := 2;
  573.   if LowerLimit = UpperLimit then
  574.     Error := 3;
  575.   if Error = 0 then
  576.   begin
  577.     InitializeStack(Values);
  578.     FillChar(DerivValues, SizeOf(DerivValues), 0);
  579.     DerivValues[1] := UserFunction1(LowerLimit, XInitial, FuncPtr);
  580.     Spacing := (UpperLimit - LowerLimit) / NumIntervals;
  581.     T := LowerLimit;
  582.     X := XInitial;
  583.     Push(Values, T, X);
  584.   end;
  585. end; { procedure TestAndInitialize }
  586.  
  587. procedure GetFirstThreeValues(Spacing     : Float;
  588.                           var T           : Float;
  589.                           var X           : Float;
  590.                           var Values      : Stack;
  591.                           var DerivValues : TNFourVector);
  592.  
  593. {---------------------------------------------------------------}
  594. {- Input: Spacing, T, X, Values, DerivValues                   -}
  595. {- Output: T, X, Values, DerivValues                           -}
  596. {-                                                             -}
  597. {- This procedure uses the Runge-Kutta one step method to      -}
  598. {- solve for the first three (T + i-Spacing, 1<=i<=3) Values.  -}
  599. {- The derivative at these points are also calculated          -}
  600. {- (DerivValues).                                              -}
  601. {---------------------------------------------------------------}
  602.  
  603. var
  604.   Index : integer;
  605.   HalfSpacing : Float;
  606.   F1, F2, F3, F4 : Float;
  607.  
  608. begin
  609.   HalfSpacing := Spacing / 2;
  610.   for Index := 1 to 3 do
  611.   begin
  612.     F1 := Spacing * UserFunction1(T, X, FuncPtr);
  613.     F2 := Spacing * UserFunction1(T + HalfSpacing, X + F1 / 2, FuncPtr);
  614.     F3 := Spacing * UserFunction1(T + HalfSpacing, X + F2 / 2, FuncPtr);
  615.     F4 := Spacing * UserFunction1(T + Spacing, X + F3, FuncPtr);
  616.     X := X + (F1 + 2 * F2 + 2 * F3 + F4) / 6;
  617.     T := T + Spacing;
  618.     Push(Values, T, X);
  619.     DerivValues[Index + 1] := UserFunction1(T, X, FuncPtr);
  620.   end;
  621. end; { procedure GetFirstThreeValues }
  622.  
  623. procedure GetValues(NumReturn    : integer;
  624.                     NumIntervals : integer;
  625.                 var Values       : Stack;
  626.                 var TValues      : TNvector;
  627.                 var XValues      : TNvector);
  628.  
  629. {-------------------------------------------------------------}
  630. {- Input: NumReturn, NumIntervals, Values                    -}
  631. {- Output: TValues, XValues                                  -}
  632. {-                                                           -}
  633. {- This procedure fills in the arrays TValues and XValues    -}
  634. {- with data from the Values stack. For example, if there    -}
  635. {- are five times as many items on the stack (NumIntervals)  -}
  636. {- as should be returned (NumReturn) then every fifth value  -}
  637. {- from the stack will be returned.                          -}
  638. {-------------------------------------------------------------}
  639.  
  640. var
  641.   Index, Term : integer;
  642.   T, X : Float;
  643.  
  644. begin
  645.   Term := NumIntervals;
  646.   for Index := NumReturn downto 0 do
  647.   begin
  648.     Pop(Values, T, X);
  649.     Term := Pred(Term);
  650.     while (Term / NumIntervals >= Index/NumReturn) and (Term >= 0) do
  651.     begin
  652.       Pop(Values, T, X);
  653.       Term := Pred(Term);
  654.     end;
  655.     TValues[Index] := T;
  656.     XValues[Index] := X;
  657.   end;
  658. end; { procedure GetValues }
  659.  
  660. begin { procedure Adams }
  661.   TestAndInitialize(LowerLimit, UpperLimit, XInitial, NumIntervals,
  662.                     NumReturn, Values, DerivValues, Spacing, T, X, Error);
  663.   if Error = 0 then
  664.   begin
  665.     GetFirstThreeValues(Spacing, T, X, Values, DerivValues);
  666.     for Index := 4 to NumIntervals do
  667.     begin
  668.       { Apply predictor }
  669.       NewX := X + Spacing / 24 * (55 * DerivValues[4] - 59 * DerivValues[3] +
  670.                   37 * DerivValues[2] - 9 * DerivValues[1]);
  671.  
  672.       T := T + Spacing;
  673.       DerivValues[1] := DerivValues[2];
  674.       DerivValues[2] := DerivValues[3];
  675.       DerivValues[3] := DerivValues[4];
  676.       DerivValues[4] := UserFunction1(T, NewX, FuncPtr);
  677.  
  678.       { Apply corrector }
  679.       NewX := X + Spacing / 24 * (9 * DerivValues[4] + 19 * DerivValues[3] -
  680.               5 * DerivValues[2] + DerivValues[1]);
  681.       X := NewX;
  682.       Push(Values, T, X);
  683.     end;
  684.  
  685.     GetValues(NumReturn, NumIntervals, Values, TValues, XValues);
  686.   end;
  687. end; { procedure Adams }
  688.  
  689. procedure InitialCond2ndOrder{(LowerLimit  : Float;
  690.                               UpperLimit   : Float;
  691.                               InitialValue : Float;
  692.                               InitialDeriv : Float;
  693.                               NumReturn    : integer;
  694.                               NumIntervals : integer;
  695.                           var TValues      : TNvector;
  696.                           var XValues      : TNvector;
  697.                           var XDerivValues : TNvector;
  698.                           var Error        : byte;
  699.                               FuncPtr      : Pointer)};
  700.  
  701. type
  702.   Ptr = ^Data;
  703.  
  704.   Data = record
  705.            T, X, xPrime : Float;
  706.            Next : Ptr;
  707.          end;
  708.  
  709. var
  710.   Values : Ptr;                          { Pointer to the stack  }
  711.   Spacing, HalfSpacing : Float;           { Size of each subinterval  }
  712.   Index : integer;
  713.   F1x, F2x, F3x, F4x,
  714.   F1xPrime, F2xPrime,
  715.   F3xPrime, F4xPrime : Float;             { Iteration variables  }
  716.   T, X, xPrime : Float;                   { Values pushed and pulled  }
  717.                                          { from the stack.           }
  718.  
  719.  
  720. procedure InitializeStack(var Values : Ptr);
  721. begin
  722.   Values := nil;
  723. end; { procedure InitializeStack }
  724.  
  725. procedure Push(var Values      : Ptr;
  726.                    TValue      : Float;
  727.                    XValue      : Float;
  728.                    XPrimeValue : Float);
  729.  
  730. {----------------------------------}
  731. {- Input: Values, TValue, XValue  -}
  732. {-        XPrimeValue             -}
  733. {- Output: Values                 -}
  734. {-                                -}
  735. {- Push data onto the Stack       -}
  736. {----------------------------------}
  737.  
  738. var
  739.   NewNode : Ptr;
  740.  
  741. begin
  742.   New(NewNode);
  743.   NewNode^.T := TValue;
  744.   NewNode^.X := XValue;
  745.   NewNode^.XPrime := XPrimeValue;
  746.   NewNode^.Next := Values;
  747.   Values := NewNode;
  748. end; { procedure Push }
  749.  
  750. procedure Pop(var Values      : Ptr;
  751.               var TValue      : Float;
  752.               var XValue      : Float;
  753.               var XPrimeValue : Float);
  754.  
  755. {----------------------------------}
  756. {- Input: Values, TValue, XValue  -}
  757. {-        XPrimeValue             -}
  758. {- Output: Values                 -}
  759. {-                                -}
  760. {- Pop Data from the Stack        -}
  761. {----------------------------------}
  762.  
  763. var
  764.   OldNode : Ptr;
  765.  
  766. begin
  767.   OldNode := Values;
  768.   Values := OldNode^.Next;
  769.   TValue := OldNode^.T;
  770.   XValue := OldNode^.X;
  771.   XPrimeValue := OldNode^.xPrime;
  772.   Dispose(OldNode);
  773. end; { procedure Pop }
  774.  
  775. procedure TestAndInitialize(LowerLimit   : Float;
  776.                             UpperLimit   : Float;
  777.                         var NumIntervals : integer;
  778.                             NumReturn    : integer;
  779.                         var Values       : Ptr;
  780.                         var Spacing      : Float;
  781.                         var Error        : byte);
  782.  
  783. {---------------------------------------------------------------}
  784. {- Input: LowerLimit, UpperLimit, NumIntervals, NumReturn      -}
  785. {- Output: Values, Spacing, Error                              -}
  786. {-                                                             -}
  787. {- This procedure initializes the above variables.  Values     -}
  788. {- is initialized to NIL. Spacing is initialized to            -}
  789. {-        (UpperLimit - LowerLimit)/NumIntervals.              -}
  790. {- NumIntervals and NumReturn are checked for errors (they     -}
  791. {- must be greater than zero).                                 -}
  792. {---------------------------------------------------------------}
  793.  
  794. begin
  795.   Error := 0;
  796.   if NumReturn <= 0 then
  797.     Error := 1;
  798.   if NumIntervals < NumReturn then
  799.     Error := 2;
  800.   if LowerLimit = UpperLimit then
  801.     Error := 3;
  802.   if Error = 0 then
  803.   begin
  804.     InitializeStack(Values);
  805.     Spacing := (UpperLimit - LowerLimit)/NumIntervals;
  806.   end;
  807. end; { procedure TestAndInitialize }
  808.  
  809. procedure GetValues(NumReturn    : integer;
  810.                     NumIntervals : integer;
  811.                 var Values       : Ptr;
  812.                 var TValues      : TNvector;
  813.                 var XValues      : TNvector;
  814.                 var XDerivValues : TNvector);
  815.  
  816. {-------------------------------------------------------------}
  817. {- Input: NumReturn, NumIntervals, Values                    -}
  818. {- Output: TValues, XValues, XDerivValues                    -}
  819. {-                                                           -}
  820. {- This procedure fills in the arrays TValues, XValues and   -}
  821. {- XDeriValues with data from the Values stack. For example, -}
  822. {- if there are five times as many items on the stack        -}
  823. {- (NumIntervals) as should be returned (NumReturn) then     -}
  824. {- every fifth value from the stack will be returned.        -}
  825. {-------------------------------------------------------------}
  826.  
  827. var
  828.   Index, Term : integer;
  829.   T, X, xPrime : Float;
  830.  
  831. begin
  832.   Term := NumIntervals;
  833.   for Index := NumReturn downto 0 do
  834.   begin
  835.     Pop(Values, T, X, xPrime);
  836.     Term := Pred(Term);
  837.     while (Term / NumIntervals >= Index / NumReturn) and (Term >= 0) do
  838.     begin
  839.       Pop(Values, T, X, xPrime);
  840.       Term := Pred(Term);
  841.     end;
  842.     TValues[Index] := T;
  843.     XValues[Index] := X;
  844.     XDerivValues[Index] := xPrime;
  845.   end;
  846. end; { procedure GetValues }
  847.  
  848. begin  { procedure InitialCond2ndOrder }
  849.   TestAndInitialize(LowerLimit, UpperLimit,
  850.                     NumIntervals, NumReturn, Values, Spacing, Error);
  851.   if Error = 0 then
  852.   begin
  853.     T := LowerLimit;
  854.     X := InitialValue;
  855.     xPrime := InitialDeriv;
  856.     Push(Values, T, X, xPrime);
  857.     HalfSpacing := Spacing / 2;
  858.     for Index := 1 to NumIntervals do
  859.     begin
  860.       F1x := Spacing * xPrime;
  861.       F1xPrime := Spacing * UserFunction2(T, X, xPrime, FuncPtr);
  862.       F2x := Spacing * (xPrime + 0.5 * F1xPrime);
  863.       F2xPrime := Spacing * UserFunction2(T + HalfSpacing, X + 0.5 * F1x,
  864.                   xPrime + 0.5 * F1xPrime, FuncPtr);
  865.       F3x := Spacing * (xPrime + 0.5 * F2xPrime);
  866.       F3xPrime := Spacing * UserFunction2(T + HalfSpacing, X + 0.5 * F2x,
  867.                   xPrime + 0.5 * F2xPrime, FuncPtr);
  868.       F4x := Spacing * (xPrime + F3xPrime);
  869.       F4xPrime := Spacing * UserFunction2(T + Spacing, X + F3x, xPrime +
  870.                                           F3xPrime, FuncPtr);
  871.       X := X + (F1x + 2 * F2x + 2 * F3x + F4x) / 6;
  872.       xPrime := xPrime + (F1xPrime + 2 * F2xPrime + 2 * F3xPrime + F4xPrime) / 6;
  873.       T := T + Spacing;
  874.       Push(Values, T, X, xPrime);
  875.     end;
  876.     GetValues(NumReturn, NumIntervals, Values, TValues, XValues, XDerivValues);
  877.   end;
  878. end; { procedure InitialCond2ndOrder }
  879.  
  880. procedure InitialCondition{(Order         : integer;
  881.                            LowerLimit     : Float;
  882.                            UpperLimit     : Float;
  883.                            InitialValues  : TNvector;
  884.                            NumReturn      : integer;
  885.                            NumIntervals   : integer;
  886.                        var SolutionValues : TNmatrix;
  887.                        var Error          : byte;
  888.                            FuncPtr        : Pointer)};
  889.  
  890. type
  891.   Ptr = ^Data;
  892.  
  893.   Data = record
  894.            Values : TNvector;
  895.            Next : Ptr;
  896.          end;
  897.  
  898. var
  899.   ValuesStack : Ptr;                     { Pointer to the stack  }
  900.   Spacing, HalfSpacing : Float;           { Size of each subinterval  }
  901.   Index : integer;
  902.   Term : integer;
  903.   F1 : TNvector;
  904.   F2 : TNvector;
  905.   F3 : TNvector;
  906.   F4 : TNvector;
  907.   CurrentValues : TNvector;
  908.   TempValues : TNvector;
  909.  
  910. procedure InitializeStack(var ValuesStack : Ptr);
  911. begin
  912.   ValuesStack := nil;
  913. end; { procedure InitializeStack }
  914.  
  915. procedure Push(var ValuesStack   : Ptr;
  916.                var CurrentValues : TNvector);
  917.  
  918. {--------------------------------------}
  919. {- Input: ValuesStack, CurrentValues  -}
  920. {- Output: ValuesStack                -}
  921. {-                                    -}
  922. {- Push data onto the Stack           -}
  923. {--------------------------------------}
  924.  
  925. var
  926.   NewNode : Ptr;
  927.  
  928. begin
  929.   New(NewNode);
  930.   NewNode^.Values := CurrentValues;
  931.   NewNode^.Next := ValuesStack;
  932.   ValuesStack := NewNode;
  933. end; { procedure Push }
  934.  
  935. procedure Pop(var ValuesStack   : Ptr;
  936.               var CurrentValues : TNvector);
  937.  
  938. {--------------------------------------}
  939. {- Input: ValuesStack, CurrentValues  -}
  940. {- Output: ValuesStack                -}
  941. {-                                    -}
  942. {- Pop data from the Stack            -}
  943. {--------------------------------------}
  944.  
  945. var
  946.   OldNode : Ptr;
  947.  
  948. begin
  949.   OldNode := ValuesStack;
  950.   ValuesStack := OldNode^.Next;
  951.   CurrentValues := OldNode^.Values;
  952. end; { procedure Pop }
  953.  
  954. procedure TestAndInitialize(Order        : integer;
  955.                             LowerLimit   : Float;
  956.                             UpperLimit   : Float;
  957.                         var NumIntervals : integer;
  958.                             NumReturn    : integer;
  959.                         var ValuesStack  : Ptr;
  960.                         var Spacing      : Float;
  961.                         var Error        : byte);
  962.  
  963. {-----------------------------------------------------------------}
  964. {- Input: Order, LowerLimit, UpperLimit, NumIntervals, NumReturn -}
  965. {- Output: ValuesStack, Spacing, Error                           -}
  966. {-                                                               -}
  967. {- This procedure initializes the above variables. ValuesStack   -}
  968. {- is initialized to NIL. Spacing is initialized to              -}
  969. {-        (UpperLimit - LowerLimit)/NumIntervals.                -}
  970. {- Order, NumIntervals, and NumReturn are checked for errors;    -}
  971. {- they must be greater than zero.                               -}
  972. {-----------------------------------------------------------------}
  973.  
  974. begin
  975.   Error := 0;
  976.   if NumReturn <= 0 then
  977.     Error := 1;
  978.   if NumIntervals < NumReturn then
  979.     Error := 2;
  980.   if Order < 1 then
  981.     Error := 3;
  982.   if LowerLimit = UpperLimit then
  983.     Error := 4;
  984.   if Error = 0 then
  985.   begin
  986.     InitializeStack(ValuesStack);
  987.     Spacing := (UpperLimit - LowerLimit) / NumIntervals;
  988.   end;
  989. end; { procedure TestAndInitialize }
  990.  
  991. procedure GetValues(NumReturn      : integer;
  992.                     NumIntervals   : integer;
  993.                 var ValuesStack    : Ptr;
  994.                 var SolutionValues : TNmatrix);
  995.  
  996. {----------------------------------------------------------}
  997. {- Input: NumReturn, NumIntervals, ValuesStack            -}
  998. {- Output: SolutionValues                                 -}
  999. {-                                                        -}
  1000. {- This procedure fills in the matrix SolutionValues      -}
  1001. {- with data from the values stack. For example,          -}
  1002. {- if there are five times as many items on the stack     -}
  1003. {- (NumIntervals) as should be returned (NumReturn) then  -}
  1004. {- every fifth value from the stack will be returned.     -}
  1005. {----------------------------------------------------------}
  1006.  
  1007. var
  1008.   Index, Term : integer;
  1009.   CurrentValues : TNvector;
  1010.  
  1011. begin
  1012.   Term := NumIntervals;
  1013.   for Index := NumReturn downto 0 do
  1014.   begin
  1015.     Pop(ValuesStack, CurrentValues);
  1016.     Term := Pred(Term);
  1017.     while (Term / NumIntervals >= Index / NumReturn) and (Term >= 0) do
  1018.     begin
  1019.       Pop(ValuesStack, CurrentValues);
  1020.       Term := Pred(Term);
  1021.     end;
  1022.     SolutionValues[Index] := CurrentValues;
  1023.   end;
  1024. end; { procedure GetValues }
  1025.  
  1026. begin { procedure InitialCondition }
  1027.   TestAndInitialize(Order, LowerLimit, UpperLimit,
  1028.                     NumIntervals, NumReturn, ValuesStack, Spacing, Error);
  1029.   if Error = 0 then
  1030.   begin
  1031.     CurrentValues := InitialValues;
  1032.     Push(ValuesStack, CurrentValues);
  1033.     HalfSpacing := Spacing / 2;
  1034.     for Index := 1 to NumIntervals do
  1035.     begin
  1036.       { First step - calculate F1  }
  1037.       for Term := 1 to Order - 1 do
  1038.         F1[Term] := Spacing * CurrentValues[Term + 1];
  1039.       F1[Order] := Spacing * UserFunction3(CurrentValues, FuncPtr);
  1040.  
  1041.       { 2nd step - calculate F2  }
  1042.       TempValues[0] := CurrentValues[0] + HalfSpacing;
  1043.       for Term := 1 to Order do
  1044.         TempValues[Term] := CurrentValues[Term] + 0.5 * F1[Term];
  1045.       for Term := 1 to Order - 1 do
  1046.         F2[Term] := Spacing * TempValues[Term + 1];
  1047.       F2[Order] := Spacing * UserFunction3(TempValues, FuncPtr);
  1048.  
  1049.       { Third step - calculate F3  }
  1050.       for Term := 1 to Order do
  1051.         TempValues[Term] := CurrentValues[Term] + 0.5 * F2[Term];
  1052.       for Term := 1 to Order - 1 do
  1053.         F3[Term] := Spacing * TempValues[Term + 1];
  1054.       F3[Order] := Spacing * UserFunction3(TempValues, FuncPtr);
  1055.  
  1056.       { Fourth step - calculate F4  }
  1057.       TempValues[0] := CurrentValues[0] + Spacing;
  1058.       for Term := 1 to Order do
  1059.         TempValues[Term] := CurrentValues[Term] + F3[Term];
  1060.       for Term := 1 to Order - 1 do
  1061.         F4[Term] := Spacing * TempValues[Term + 1];
  1062.       F4[Order] := Spacing * UserFunction3(TempValues, FuncPtr);
  1063.  
  1064.       { Combine F1, F2, F3, and F4 to get  }
  1065.       { the solution at this mesh point    }
  1066.       CurrentValues[0] := CurrentValues[0] + Spacing;
  1067.       for Term := 1 to Order do
  1068.         CurrentValues[Term] := CurrentValues[Term] +
  1069.                                (F1[Term] + 2 * F2[Term] +
  1070.                                2 * F3[Term] + F4[Term]) / 6;
  1071.  
  1072.       Push(ValuesStack, CurrentValues);
  1073.     end;
  1074.     GetValues(NumReturn, NumIntervals, ValuesStack, SolutionValues);
  1075.   end;
  1076. end; { procedure InitialCondition }
  1077.