home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e032 / 3.ddi / FILES / NUMERICA.PAK / BUTCHER.M < prev    next >
Encoding:
Text File  |  1992-07-29  |  6.2 KB  |  206 lines
  1.  
  2. (* :Name: NumericalMath`Butcher` *)
  3.  
  4. (* :Title: Runge-Kutta Order Conditions *)
  5.  
  6. (* :Author: Jerry B. Keiper *)
  7.  
  8. (* :Summary:
  9.     Gives the order conditions which any Runge-Kutta method must
  10.     satisfy to be of any particular order.  Works for both implicit
  11.     and explicit methods.
  12. *)
  13.  
  14. (* :Context: NumericalMath`Butcher` *)
  15.  
  16. (* :Package Version: 2.0 *)
  17.  
  18. (* :Copyright: Copyright 1990,  Wolfram Research, Inc.
  19.  
  20.     Permission is hereby granted to modify and/or make copies of
  21.     this file for any purpose other than direct profit, or as part
  22.     of a commercial product, provided this copyright notice is left
  23.     intact.  Sale, other than for the cost of media, is prohibited.
  24.  
  25.     Permission is hereby granted to reproduce part or all of
  26.     this file, provided that the source is acknowledged.
  27. *)
  28.  
  29. (* :History:
  30.     Updated by Jerry B. Keiper, December, 1990.
  31.     Version 1.2 by Jerry B. Keiper, 1989.
  32. *)
  33.  
  34. (* :Keywords:
  35.     Butcher trees, Runge-Kutta, order conditions
  36. *)
  37.  
  38. (* :Source:
  39.     John C. Butcher: The Numerical Analysis of Ordinary Differential
  40.     Equations: Runge-Kutta and General Linear Methods, John Wiley &
  41.     Sons, New York, 1987
  42. *)
  43.  
  44. (* :Mathematica Version: 2.0 *)
  45.  
  46. (* :Limitation:
  47.     This package only finds the order conditions; it does not try
  48.     to solve them for the required coefficients.
  49.  
  50.     Combinatorial explosion effectively limits the order for which
  51.     order conditions can be found,  e.g., 7813 (highly nonlinear)
  52.     conditions are required for a 12th-order method.
  53.  
  54.     This is a rather rudimentary package and could be improved
  55.     considerably.  The functions r (order), gamma (density),
  56.     and phi are defined as private functions and could be exposed
  57.     to the user.  Furthermore, tools for aiding in the solution
  58.     of the order conditions could be useful.
  59. *)
  60.  
  61. (* :Warnings:
  62.     Derivative is unprotected so that rules can be attached to it.
  63.     It is protected after.
  64. *)
  65.  
  66. (* :Discussion:
  67.     The method is due to John C. Butcher and is explained in chapter
  68.     3 of Butcher's book.  It relies heavily on graph theory and
  69.     recursive decomposition of Butcher trees.  Pattern matching is
  70.     used extensively.
  71. *)
  72.  
  73. BeginPackage["NumericalMath`Butcher`"];
  74.  
  75. RungeKuttaOrderConditions::usage =
  76. "RungeKuttaOrderConditions[p, s] gives a list of the order conditions that
  77. any s-stage Runge-Kutta method of order p must satisfy.  These conditions are
  78. expressed in terms of the global variables a, b, c, and the resulting
  79. Runge-Kutta method is\n
  80. \n
  81. Y[i,x] = y[x0] + (x - x0) Sum[a[i,j] f[Y[j,x]], {j, 1, s}], i = 1, ..., s\n
  82. \n
  83. Y[x] = y[x0] + (x - x0) Sum[b[j] f[Y[j,x]], {j, 1, s}]\n
  84. \n
  85. where c[i] == Sum[a[i, j], {j, 1, s}]."
  86.  
  87. a::usage =
  88. "a[i,j] is the coefficient of f[Y[j, x]] in the formula for Y[i, x] for the
  89. Runge-Kutta method being considered."
  90.  
  91. b::usage =
  92. "b[j] is the coefficient of f[Y[j, x]] in the formula for Y[x] for the
  93. Runge-Kutta method being considered."
  94.  
  95. c::usage =
  96. "c[i] is a notational convenience representing Sum[a[i, j], {j, 1, s}] in
  97. the order conditions for the Runge-Kutta method being considered."
  98.  
  99. Begin["NumericalMath`Butcher`Private`"];
  100.  
  101. (* Attach rules to Derivative must unprotect it *)
  102. Unprotect[ Derivative]
  103.  
  104. (* general rules for ' *)
  105. (a_Integer)' := 0;
  106. (a_ + b_)' := a' + b';
  107. (a_ b_)' := Expand[a' b + a b'];
  108. (a_^n_)' := Expand[n a^(n-1) a'];
  109.  
  110. (* rules for recursively evaluating y', y'', y''', ...
  111. note: derivatives of f are implicit in the number of arguments it has. *)
  112. y' = f;
  113. f' = f[f];
  114. (f[g_])' := (f[g])' = f[Expand[g f]] + f[g'];
  115.  
  116. (* make f multi-linear. *)
  117. f[a_ + b_] := f[a] + f[b];
  118. f[n_Integer a_] := n f[a];
  119.  
  120. (* rewrite a derivative of y as the list of trees appearing in the expression. *)
  121. trees[f] = {f};
  122. trees[f[f]] = {f[f]};
  123. trees[a_Plus] := Apply[List, Map[trees, a]]
  124. trees[n_. a_f] := a
  125.  
  126. (* define the order of a tree. *)
  127. r[t_f] := (r[t] = 1 + Map[r, Apply[Plus, t[[1]]]]) /; Head[t[[1]]] === Times;
  128. r[t_f] := r[t] = 1 + r[t[[1]]];
  129. r[f^n_.] := n;
  130. r[(t_f)^n_Integer] := r[t^n] = n r[t];
  131.  
  132. (* define the density of a tree. *)
  133. gamma[t_f] := (gamma[t] = r[t] Map[gamma, t[[1]]]) /; Head[t[[1]]] === Times;
  134. gamma[t_f] := gamma[t] = r[t] gamma[t[[1]]];
  135. gamma[f^n_.] = 1;
  136. gamma[(t_f)^n_Integer] := gamma[t]^n;
  137.  
  138. (* the function phi is more complicated, since we need bookkeeping for the
  139. various indices of summation.  The argument s must be an integer and is
  140. the upper bound of the summation in the order conditions.
  141. phi[l, s, forest] recursively produces a list containing the largest
  142. summation variable needed for that forest and those factors associated
  143. with the subtree rooted at vertex l and having as its own subtrees those
  144. listed in forest.  phi[l, s, tree] does the same except it is for a 
  145. single subtree that appears in the forest. *)
  146.  
  147. (* find the list of subtrees (one level up) of a tree. *)
  148. split[tree_f] :=
  149.     Module[{i, j, sub = tree[[1]]},
  150.         sub = If[Head[sub] === Times, Apply[List, sub], {sub}];
  151.         Do[If[Head[sub[[i]]] === Power, 
  152.             sub[[i]] = Table[sub[[i, 1]], {j, sub[[i, 2]]}]],
  153.             {i, Length[sub]}];
  154.         Flatten[sub]
  155.         ];
  156.  
  157. (* special case for the tree with one vertex. *)
  158. phi[s_, f] = Sum[b[i[0]], {i[0], s}];
  159.  
  160. (* special case for the root of the tree. *)
  161. phi[s_, tree_f] := 
  162.     Module[{prod, j, k, ranges},
  163.         prod = phi[0, s, split[tree]];
  164.         k = prod[[1]];
  165.         prod = b[i[0]] prod[[2]];
  166.         ranges = Table[{i[j], s}, {j, 0, k}];
  167.         Apply[Sum, Prepend[ranges, prod]]
  168.         ];
  169.  
  170. (* a forest of sub-trees. *)
  171. phi[l_, s_, forest_List] :=
  172.     Module[{j, k=l, prod = 1, temp},
  173.         Do[    temp = phi[k, s, forest[[j]]];
  174.             prod *= temp[[2]];
  175.             k = temp[[1]],
  176.             {j, Length[forest]}];
  177.         {k, prod}
  178.         ];
  179.  
  180. (* special case for a leaf. *)
  181. phi[l_, s_, f] := {l, c[i[l]]};
  182.  
  183. (* a single subtree. *)
  184. phi[l_, s_, tree_f] :=
  185.     Module[{temp},
  186.         temp = phi[l+1, s, split[tree]];
  187.         {temp[[1]], a[i[l], i[l+1]] temp[[2]]}
  188.         ];
  189.  
  190. RungeKuttaOrderConditions[p_, s_] :=
  191.     Module[{treelist, gammalist, i, yder},
  192.         yder = y;
  193.         treelist = {};
  194.         Do[treelist = Join[treelist, trees[yder = yder']], {i,p}];
  195.         gammalist = 1/Map[gamma, treelist];
  196.         Map[Apply[Equal, #]&,
  197.             Transpose[{Map[phi[s, #]&, treelist], gammalist}]]
  198.         ]
  199.  
  200. End[ ];        (* "NumericalMath`Butcher`Private`" *)
  201.  
  202. Protect[RungeKuttaOrderConditions, a, b, c];
  203.  
  204. EndPackage[ ];        (* "NumericalMath`Butcher`" *)
  205.  
  206.