(* :Name: NumericalMath`Butcher` *)

(* :Title: Runge-Kutta Order Conditions and Butcher Tree Functions *)

(* :Author: Mark Sofroniou (original version by Jerry B. Keiper) *)

(* :Summary:
This package gives the order conditions that a Runge-Kutta method
must satisfy to be of a particular order. It works for both implicit
and explicit methods. The package also calculates and plots Butcher trees
and implements the functions of trees defined in Butcher's book.
Butcher's row and column simplifying conditions assist in the derivation
of high order methods. A more compact and efficient stage-independent
tensor notation has also been implemented. *)

(* :Context: NumericalMath`Butcher` *)

(* :Package Version: 2.1 *)

(* :Copyright: Copyright 1990-1999, Wolfram Research, Inc. *)

(* :History:
 Version 1.2 by Jerry B. Keiper, 1989.
 Updated by Jerry B. Keiper, December 1990.
 Updated by Jerry B. Keiper, August, November 1992. Implemented
 suggestions of Mark Sofroniou, Loughborough University.
 Re-written and updated by Mark Sofroniou, October 1993.
 Updated by Mark Sofroniou July 1994 (acknowledgements to
 A. Harrison and S. N. Papakostas for suggesting the tensor
 notation).
 Updated by Mark Sofroniou January 1998. Performance enhancements to
 numerical functions on trees.
*)

(* :Keywords:
 Butcher trees, Runge-Kutta, order conditions. *)

(* :Source:
 John C. Butcher: The Numerical Analysis of Ordinary Differential
 Equations: Runge-Kutta and General Linear Methods, John Wiley &
 Sons, New York, 1987.
 J. D. Lambert: Numerical Methods for Ordinary Differential Systems,
 The Initial Value Problem, John Wiley & Sons, New York, 1991.
 Other useful references are:
 E. Hairer, S. P. Norsett and G. Wanner, Solving Ordinary Differential
 Equations I: Nonstiff Problems, Second edition, Springer Verlag,
 New York, 1993.
 E. Hairer, and G. Wanner, Solving Ordinary Differential Equations II:
 Stiff and Differential Algebraic Problems, Second edition, Springer Verlag,
 New York, 1996.
 *)

(* :Mathematica Version: 3.0 *)

(* :Limitation:
 This package only derives the order conditions; it does not try
 to solve them for the required coefficients.

 Combinatorial explosion effectively limits the order for which
 order conditions can be found, e.g., 7813 (highly nonlinear)
 conditions are required for a 12th-order method.

 Butcher's simplifying conditions have been added to partly
 address this issue. *)

(* :Discussion:
 The method is due to John C. Butcher and is explained in Chapter
 3 of Butcher's book. It relies heavily on graph theory and
 recursive decomposition of Butcher trees. Pattern matching is
 used extensively. Tensor notation (see Lambert, Chapter 5.12
 for its use in Albrecht's formalism) is more efficient and
 is used for evaluation of the order conditions with numeric stages. *)

BeginPackage["NumericalMath`Butcher`", "DiscreteMath`IntegerPartitions`"];

a::usage =
"a[i,j] is the coefficient of f[Y[j, x]] in the formula for Y[i, x] of the
Runge-Kutta method being considered.  Subscript notation is used to display
components of the matrix a."

b::usage =
"b[j] is the coefficient of f[Y[j, x]] in the formula for Y[x] of the
Runge-Kutta method being considered.  Subscript notation is used to display
components of the vector b."

c::usage =
"c[i] is a notational convenience used to represent Sum[a[i, j], {j, 1, s}] in
the order conditions for the Runge-Kutta method being considered.  Subscript
notation is used to display components of the vector c."

e::usage =
"e represents the vector of length s with unit components."

f::usage =
"f is the elementary symbol used in the representation of Butcher trees."

ButcherAlpha::usage =
"ButcherAlpha[tree] gives the number of ways of labeling the vertices of
the tree tree with a totally ordered set of labels such that if (m, n) is
an edge, then m < n."

ButcherBeta::usage =
"ButcherBeta[tree] gives the number of ways of labelling the tree tree
with ButcherOrder[tree]-1 distinct labels such that the root is not
labelled, but every other vertex is labelled.  ButcherBeta[n, tree]
gives the number of ways of labelling n of the vertices of the tree with
n distinct labels such that every leaf is labelled and the root is not
labelled."

ButcherBetaBar::usage =
"ButcherBetaBar[tree] gives the number of ways of labelling the tree
tree with ButcherOrder[tree] distinct labels such that every vertex
is labelled.  ButcherBetaBar[n, tree] gives the number of ways of
labelling n of the vertices of the tree with n distinct labels such that
every leaf is labelled."

ButcherColumnConditions::usage =
"ButcherColumnConditions[p, s] gives the column simplifying conditions up to
and including order p for s stages.  ButcherColumnConditions[p] gives
the column simplifying conditions in stage-independent tensor notation.
The option RungeKuttaMethod controls the form of the result."

ButcherGamma::usage =
"ButcherGamma[tree] gives the density of the tree tree.  This is also the
reciprocal of the right-hand side of the order condition imposed by tree."

ButcherHeight::usage =
"ButcherHeight[tree] gives the height of the tree tree."

ButcherOrder::usage =
"ButcherOrder[tree] gives the order r of the tree tree."

ButcherPhi::usage =
"ButcherPhi[tree, s] gives the left-hand side of the order condition
imposed by tree of an s-stage Runge-Kutta method.  ButcherPhi[tree] gives
the weight function in stage-independent tensor notation involving dot
products and scalar multiplications.  Using this notation, a is an s by s
matrix, b and c are vectors of length s, and e is the length-s vector with
unit components.  The option RungeKuttaMethod controls the form of the result."

ButcherPlot::usage =
"ButcherPlot[tree,options] gives a plot of the tree tree.\n
ButcherPlot[{tree1, tree2, ...},options] gives an array of plots of
the various trees.  Valid options are ButcherPlotColumns,
ButcherPlotNodeSize, ButcherPlotRootSize and ButcherPlotLabel."

ButcherPlotColumns::usage = "ButcherPlotColumns is an option to ButcherPlot that
specifies the number of columns in the GraphicsArray plot of a list of trees."

ButcherPlotLabel::usage = "ButcherPlotLabel is an option to ButcherPlot that
specifies a list of plot labels, which should be the same length as the
list of trees to be plotted."

ButcherPlotNodeSize::usage = "ButcherPlotNodeSize is an option to ButcherPlot that
specifies a scaling factor to be applied to the default size for nodes in the plot."

ButcherPlotRootSize::usage = "ButcherPlotRootSize is an option to ButcherPlot that
specifies a scaling factor to be applied to the default size for encircling of tree
roots in the plot. The root is not highlighted if the value is zero."

ButcherPrincipalError::usage =
"ButcherPrincipalError[p, s] calculates the principal error for a
method of order p with s stages.  ButcherPrincipalError[p] gives
the principal error using stage-independent tensor notation.  The result
is a list of the order p+1 terms appearing in a Taylor series expansion
of the error.  The option RungeKuttaMethod controls the form of the result."

ButcherQuadratureConditions::usage =
"ButcherQuadratureConditions[p, s] gives the quadrature conditions up to
and including order p for s stages.  ButcherQuadratureConditions[p] gives
the quadrature conditions in stage-independent tensor notation.  The option
RungeKuttaMethod controls the form of the result."

ButcherRowConditions::usage =
"ButcherRowConditions[p, s] gives the row simplifying conditions up to
and including order p for s stages.  ButcherRowConditions[p] gives
the row simplifying conditions in stage-independent tensor notation.
The option RungeKuttaMethod controls the form of the result."

ButcherRowSum::usage =
"ButcherRowSum is an option to RungeKuttaOrderConditions stating whether
the row-sum conditions for the c[i] should be added to the list of order
conditions.  ButcherRowSum may be set to True or False."

ButcherSimplify::usage =
"ButcherSimplify is an option to RungeKuttaOrderConditions specifying
whether to apply Butcher's row and column simplifying conditions.  The
conditions reduce the set of necessary and sufficient trees for an order
p method with s stages.  This is especially useful for deriving high order
methods.  The form of the simplifying assumptions used depends upon whether
the method is implicit or explicit.  For Implicit methods, all trees are
reduced to quadrature type together with subsidiary row and column simplifying
conditions (which are sufficient but by no means necessary for order p).  For
explicit methods, the restriction b[2]==0 for p>=5 implies positivity of the
weights b[i].  If simplification is to be performed the number of stages
must be a number.  Valid settings are True and False (standard order
conditions)."

ButcherSigma::usage =
"ButcherSigma[tree] gives the order of the symmetry group of isomorphisms
of the tree tree with itself."

ButcherTrees::usage =
"ButcherTrees[p] gives a list of lists of the trees for any Runge-Kutta
method of order p.  The inner lists group the trees by order."

ButcherTreeCount::usage = "ButcherTreeCount[p] gives a list of the number of
trees through order p."

ButcherTreeQ::usage =
"ButcherTreeQ[tree] tests for valid syntax of the tree or trees tree."

ButcherTreeSimplify::usage =
"ButcherTreeSimplify[p, eta, xi] gives the set of trees through order p that
are not reduced by Butcher's simplifying assumptions, assuming that the
quadrature conditions through order p, the row simplifying conditions through
order eta, and the column simplifying conditions through order xi all hold.
The result is grouped by order, starting with the first non-vanishing trees."

ButcherWidth::usage =
"ButcherWidth[tree] gives the width of the tree tree."

RungeKuttaOrderConditions::usage =
"RungeKuttaOrderConditions[p, s] gives a list of lists of the order conditions
that any s-stage Runge-Kutta method of order p must satisfy.  The inner lists
group the conditions by order.  The conditions are expressed in terms of the
global variables a, b, c, and the resulting Runge-Kutta method is\n
\n
Y[i,x] = y[x0] + (x - x0) Sum[a[i,j] f[Y[j,x]], {j, 1, s}], i = 1, ..., s\n
\n
Y[x] = y[x0] + (x - x0) Sum[b[j] f[Y[j,x]], {j, 1, s}]\n
\n
where the row-sum conditions c[i] == Sum[a[i, j], {j, 1, s}] are usually
assumed to hold.  RungeKuttaOrderConditions[p] expresses the order conditions
using stage-independent tensor notation.  In this notation, a is an s by s
matrix, b and c are vectors of length s, and e is the length-s vector with
unit components.  The options ButcherRowSum, ButcherSimplify, and
RungeKuttaMethod control the form of the result."

RungeKuttaMethod::usage=
"RungeKuttaMethod is an option to ButcherPhi and related functions
specifying the type of method to be generated.  Valid settings are
Explicit, DiagonallyImplicit and Implicit."

$RungeKuttaMethod::usage=
"$RungeKuttaMethod is a global environment setting, specifying the
type of method to be generated by ButcherPhi and related functions.
Valid settings are Explicit, DiagonallyImplicit and Implicit."

Explicit::usage=
"Explicit is a setting for the option RungeKuttaMethod specifying
the type of Runge-Kutta method to be generated."

DiagonallyImplicit::usage=
"DiagonallyImplicit is a setting for the option RungeKuttaMethod
specifying the type of Runge-Kutta method to be generated."

Implicit::usage=
"Explicit is a setting for the option RungeKuttaMethod specifying
the type of Runge-Kutta method to be generated."

SymbolicTable::usage =
"SymbolicTable is an alternative to Table that allows symbolic ranges."

Unprotect[a, b, c, e, f, ButcherAlpha, ButcherBeta, ButcherBetaBar,
ButcherColumnConditions, ButcherGamma, ButcherHeight, ButcherOrder,
ButcherPhi, ButcherPlot, ButcherPlotColumns, ButcherPlotLabel,
ButcherPlotNodeSize, ButcherPlotRootSize, ButcherPrincipalError,
ButcherQuadratureConditions, ButcherRowConditions, ButcherRowSum,
ButcherSigma, ButcherSimplify, ButcherTrees, ButcherTreeCount,
ButcherTreeQ, ButcherTreeSimplify, ButcherWidth, Explicit,
DiagonallyImplicit, Implicit, RungeKuttaMethod, RungeKuttaOrderConditions,
SymbolicTable];

Begin["`Private`"];

(* Global environment setting specifying the type of method to be
 generated by the weight function. *)

$RungeKuttaMethod = Implicit;

CheckMethod[mtd_,func_]:=
  If[ MemberQ[{Explicit,DiagonallyImplicit,Implicit},mtd],
    True, (* Else *) Message[RungeKuttaMethod::args,func] ];

RungeKuttaMethod::args =
"The option RungeKuttaMethod in `1` did not evaluate to Explicit,
DiagonallyImplicit or Implicit";

(* Subscripted notation for formatting. *)

SetAttributes[{a,b,c},NHoldAll];

a/: Format[a[i_,j_],OutputForm]:= Subscripted[a[i,j]];
a/: Format[a[i_,j_],TeXForm]:= Subscripted[a[i,j]];
a/: MakeBoxes[a[i_,j_], form_]:=
      SubscriptBox["a",RowBox[{MakeBoxes[i,form],",",MakeBoxes[j,form]}]];

(* Reformat double subscripted variable for input. *)

MakeExpression[SubscriptBox["a",indices_RowBox],form_]:=
  MakeExpression[ RowBox[{ "a", "[", indices, "]"}], form ];

b/: Format[b[i_],TeXForm]:= Subscripted[b[i]];
b/: Format[b[i_],OutputForm]:= Subscripted[b[i]];
b/: MakeBoxes[b[i_], form_]:=
      SubscriptBox["b",MakeBoxes[i,form]];

c/: Format[c[i_],OutputForm]:= Subscripted[c[i]];
c/: Format[c[i_],TeXForm]:= Subscripted[c[i]];
c/: MakeBoxes[c[i_], form_]:=
      SubscriptBox["c",MakeBoxes[i,form]];


(* Formatting of indices. *)

charcodebase = ToCharacterCode["i"]-1;
$i[n_]:=
  With[{char = FromCharacterCode[charcodebase+n]},
    $i[n] = ToExpression["NumericalMath`Butcher`Private`$i"<>ToString[n]];
    (Format[#,OutputForm] = char;
      Format[#,TeXForm] = char;
      MakeBoxes[#,form_] = char)& @ $i[n];
    $i[n]
  ];


(* Consider all possible combinations of lower order trees. *)

maketrees[{p_}]:= dups[p]; (* Partition with one element. *)
maketrees[p_]:= Apply[ Outer[Times,##]&, Map[dups,p] ];

(* Use multiplicity to avoid duplication. The number of trees is
 used at each order. *)

dups[{i_,1}]:= trees[i];
dups[{i:(1|2),m_}]:= trees[i]^m;
dups[{i_,m_}]:=
  dups[{i,m}]=
    With[{indx = Array[$i,m], trs = trees[i]},
      Apply[
        Flatten[ Table[ Apply[ Times, Part[ trs, indx ] ], ## ] ]&,
        Thread[ {indx, Join[ {btc[i]}, Drop[indx,-1] ] } ]
      ]
    ];

(* Trees grouped by order. Lower order trees are grafted onto a new
 root by applying f. *)

SetAttributes[{f},Listable];
trees[1] = {f};
trees[p_]:= trees[p] = f[ Flatten[ Map[ maketrees, IntegerPartitions[-1+p] ] ] ];

ButcherTrees[p_Integer?Positive]:=
  Block[{$RecursionLimit = Infinity}, Array[trees,p] ];


(* Test input syntax for a tree. *)

ButcherTreeQ[{}] = False;
ButcherTreeQ[t_]:= Block[{$RecursionLimit=Infinity},forestQ[t]];

forestQ[t_List]:= Apply[And,Map[forestQ,t]];
forestQ[f] = True;
forestQ[t_f]:= treeQ[t];
forestQ[_] = False

treeQ[f]=True;
treeQ[tree:f[t_]]:= treeQ[tree]=treeQ[t];
treeQ[t_^_Integer?Positive]:= treeQ[t];
treeQ[t_Times]:= Apply[And,Map[treeQ,Apply[List,t]]];
treeQ[_]:= False;


(* Enumeration of rooted trees. See for example section 2.3.4.4,
 exercise 2 of D. E. Knuth (1974), Fundamental Algorithms:
 The Art of Computer Programming, Volume 1 (2nd edition),
 Addison Wesley, Massachusetts. The following is a more efficient
 implementation avoiding the use of a double indexed function. *)

ButcherTreeCount[n_Integer?Positive]:= Array[btc,n];

btc[1]=1;
btc[n_]:= btc[n]=
  Quotient[
    Sum[i btc[i] Sum[btc[-i j+n],{j,Quotient[#,i]}],{i,#}],
    #
  ]& @ (-1+n);


(**** Numerical functions on trees ****)

(* Define the order of a tree. *)

SetAttributes[order, Listable];
order[f] = 1;
order[t:f[subt_]]:= order[t] = 1 + order[subt];
order[t_Times]:= Apply[Plus,Map[order,Apply[List,t]]];
order[t_^n_]:= n*order[t];

ButcherOrder[trees_?ButcherTreeQ]:=
  Block[{$RecursionLimit = Infinity}, order[trees]];

(* Define the density of a tree. *)

SetAttributes[gamma, Listable];
gamma[f] = 1;
gamma[t:f[subt_]]:= gamma[t] = order[t]*gamma[subt];
gamma[t_Times]:= Map[gamma,t];
gamma[t_^n_]:= gamma[t]^n;

ButcherGamma[trees_?ButcherTreeQ]:=
  Block[{$RecursionLimit = Infinity}, gamma[trees]];

(* Define the symmetry of a tree *)

SetAttributes[sigma, Listable];
sigma[f] = 1;
sigma[t:f[subt_]]:= sigma[t] = sigma[subt];
sigma[t_Times]:= Map[sigma,t];
sigma[t_^n_]:= n!*sigma[t]^n;

ButcherSigma[trees_?ButcherTreeQ]:=
  Block[{$RecursionLimit = Infinity}, sigma[trees]];

(* Define the number of monotonic labellings *)

(* We use a variant of Hairer's definition for alpha, since it is
 more efficient than Butcher's definition in terms of order, gamma and
 sigma. *)

(* Multinomial coefficient for equal arguments m, repeated n times and
 divided by n! *)

mult[1, _]:= 1;
mult[m_, n_]:= mult[m, n]= Quotient[Pochhammer[n+1, n*(m-1)], m!^n];

SetAttributes[alpha, Listable];
alpha[f] = 1;
alpha[t:f[subt_]]:= alpha[t] = alpha[subt];
alpha[t_Times]:=
  Times[
    Apply[ Multinomial, order[Apply[List, t]] ],
    Map[ alpha, t ]
  ];
alpha[t_^n_]:= Times[ mult[order[t], n], alpha[t]^n ];

ButcherAlpha[trees_?ButcherTreeQ]:=
  Block[{$RecursionLimit = Infinity}, alpha[trees] ];

(* Define the height f a tree *)

SetAttributes[height, Listable];
height[f] = 1;
height[t:f[subt_]]:= height[t] = 1 + height[subt];
height[t_Times]:= Max[height[Apply[List,t]]];
height[t_^_]:= height[t];

ButcherHeight[trees_?ButcherTreeQ]:=
  Block[{$RecursionLimit = Infinity}, height[trees]];

(* Define the width of a tree *)

SetAttributes[{width,widthbar}, Listable];
width[f] = 0;
width[t_f]:= widthbar[t];
widthbar[f] = 1;
widthbar[t:f[subt_]]:= widthbar[t] = widthbar[subt];
widthbar[t_Times]:= Apply[Plus,widthbar[Apply[List,t]]];
widthbar[t_^n_]:= n*widthbar[t];

ButcherWidth[trees_?ButcherTreeQ]:=
  Block[{$RecursionLimit = Infinity}, widthbar[trees]];

(* Define the number of ways of labelling the tree *)

SetAttributes[beta, Listable];
beta[f] = 1;
beta[t_]:= beta[t] = Quotient[ (order[t]-1)!, sigma[t] ];
beta[0, f] = 1;
beta[0, _] = 0;
beta[n_, t_]:=
  beta[n, t] =
    Quotient[
      n!*(order[t]-1-width[t])!,
      sigma[t]*(n-width[t])!*(order[t]-1-n)!
    ];

ButcherBeta[tree_?ButcherTreeQ]:=
  Block[{$RecursionLimit = Infinity}, beta[tree]];

ButcherBeta[n:(0|_Integer?Positive), tree_?ButcherTreeQ]:=
  Block[{$RecursionLimit = Infinity}, beta[n, tree]];


SetAttributes[betabar, Listable];
betabar[f] = 1;
betabar[t_]:= order[t]*beta[t];
betabar[0, _] = 0;
betabar[n_, t_]:=
  betabar[n, t] =
    Quotient[
      n!*(order[t]-widthbar[t])!,
      sigma[t]*(n-widthbar[t])!*(order[t]-n)!
    ];

ButcherBetaBar[tree_?ButcherTreeQ]:=
  Block[{$RecursionLimit = Infinity}, betabar[tree]];

ButcherBetaBar[n:(0|_Integer?Positive), tree_?ButcherTreeQ]:=
  Block[{$RecursionLimit = Infinity}, betabar[n, tree]];


(**** Order conditions ****)

(* Butcher's index notation. *)

(* The function phi requires bookkeeping for the various indices of summation.
 It is more efficient to collect components in individual sums, since the
 conditions are maintained in factored form when evaluated numerically. *)

(* Indexing Function. *)

indx[i_, symb, Explicit]:= {i,2,symb};
indx[i_,  j_,  Explicit]:= {i,2,j-1};
indx[i_,  j_,  DiagonallyImplicit]:= {i,j};
indx[i_,   _,  Implicit]:= {i,symb};

(* Terminal vertices. *)

sumterms[depth_,_,f^n_.] := c[$i[depth]]^n;

(* Non-terminal vertices. Current depth of the tree is the ith index in
 the summation. *)

sumterms[depth_,mtd_,f[subt_]^n_.]:=
  With[{sumindex = $i[++index]},
    Sum[Evaluate[a[$i[depth],sumindex] sumterms[index,mtd,subt],
          indx[sumindex,$i[depth],mtd]]]^n
  ];

(* Multiple branching. *)

sumterms[depth_,mtd_,subt_Times]:= Map[sumterms[depth,mtd,#]&, subt];

SetAttributes[iphi, Listable];

(* Special case for the tree with one vertex. *)

iphi[f,_] := Sum[Evaluate[b[$i[1]],{$i[1],symb}]];

(* General case for a rooted tree. index is the jth index in the
 summation. *)

iphi[f[subt_],mtd_]:=
  Block[{index=1},
    Sum[Evaluate[b[$i[1]] sumterms[1,mtd,subt],
          indx[$i[1],symb,mtd]]]
  ];

(* Butcher's index notation. Definition enables the storage of symbolic
 sums and substitutes for summation values *)

ButcherPhi[trees_?ButcherTreeQ, stages_Symbol,opts___?OptionQ]:=
  Block[{$RecursionLimit=Infinity},
    Module[{mtd},
      mtd = RungeKuttaMethod /. {opts} /. Options[ButcherPhi];
      (iphi[trees,mtd] /. symb->stages) /; CheckMethod[mtd,ButcherPhi]
    ]
  ];


(* Tensor notation. *)

(* Terminal vertices. *)

dotterms[f] = c;

(* Non-terminal vertices. *)

dotterms[t:f[subt_]]:= dotterms[t]=a.dotterms[subt];

(* Multiple branching. *)

dotterms[t_Times]:= Map[dotterms, t];
dotterms[t_^n_]:= dotterms[t]^n;

SetAttributes[tphi, Listable];

(* Special case for the tree with one vertex. *)

tphi[f] = b.e;

(* General case for a rooted tree. *)

tphi[f[subt_]]:= b.dotterms[subt];

ButcherPhi[trees_?ButcherTreeQ]:=
  Block[{$RecursionLimit=Infinity}, tphi[trees]];


(* Create local elements for phi recursion with numeric stages. *)

nphi[trees_,s_,mtd:(Implicit|DiagonallyImplicit|Explicit)]:=
  Module[{amat,bvec = Array[b,s],cvec,dotterms,phi},

    Switch[mtd,
      Explicit,
      amat = Table[If[i<=j,0,a[i,j]],{i,s},{j,s}];
      cvec = Join[{0},Array[c,s-1,2]],
      DiagonallyImplicit,
      amat = Table[If[i<j,0,a[i,j]],{i,s},{j,s}]; cvec = Array[c,s],
      Implicit,
      amat = Array[a,{s,s}]; cvec = Array[c,s]
    ];

    dotterms[f] = cvec;
    dotterms[t:f[subt_]]:= dotterms[t]=amat.dotterms[subt];
    dotterms[t_Times]:= Map[dotterms,t];
    dotterms[t_^n_]:= dotterms[t]^n;

    SetAttributes[phi, Listable];
    phi[f] = bvec.Table[1,{s}];
    phi[f[subt_]]:= bvec.dotterms[subt];

    phi[trees]
  ];

(* Numeric stages (uses tensor formulism for efficiency). *)

Options[ButcherPhi] = {RungeKuttaMethod:>$RungeKuttaMethod};

ButcherPhi[trees_?ButcherTreeQ, s_Integer?Positive,opts___?OptionQ]:=
  Block[{$RecursionLimit=Infinity},
    Module[{mtd},
      mtd = RungeKuttaMethod /. {opts} /. Options[ButcherPhi];
      nphi[trees,s,mtd] /; CheckMethod[mtd,ButcherPhi]
    ]
  ];


(* Principal error function. *)

Options[ButcherPrincipalError] = {RungeKuttaMethod:>$RungeKuttaMethod};

ButcherPrincipalError[p_Integer?Positive,
    s:(_Integer?Positive|_Symbol),opts___?OptionQ]:=
  Module[{mtd},
    mtd = RungeKuttaMethod /. {opts} /. Options[ButcherPrincipalError];
    ((ButcherPhi[#,s,RungeKuttaMethod->mtd] - 1/ButcherGamma[#])/
       ButcherSigma[#])& @ Last[ButcherTrees[p+1]] /;
         CheckMethod[mtd,ButcherPrincipalError]
  ];

ButcherPrincipalError[p_Integer?Positive]:=
  Module[{trees = Last[ButcherTrees[p+1]]},
    (ButcherPhi[trees] - 1/ButcherGamma[trees])/ButcherSigma[trees]
  ];


(* Display a Table with symbolic range. *)

SymbolicTable[expr_, {i_, i0_:1, i1_Integer}]:= Table[expr, {i, i0, i1}];

Options[RungeKuttaOrderConditions] =
  {ButcherRowSum->False,ButcherSimplify->False,
     RungeKuttaMethod:>$RungeKuttaMethod};

RungeKuttaOrderConditions::opts =
"The option `1` in RungeKuttaOrderConditions did not evaluate to `2`";

RKOCtest =
  If[#1,
    True,
    Message[RungeKuttaOrderConditions::opts,Apply[Sequence,#2]]; False
  ]&;

RKOCmessages = {
  {ButcherRowSum, "True or False"},
  {ButcherSimplify, "False or True with s a positive integer"},
  {RungeKuttaMethod, "Explicit, DiagonallyImplicit or Implicit"}};

RKOCOptionsTest[s_,opts___]:=
  Module[{butchersimpQ, datatypes, mtd, rowsumQ},
    {rowsumQ, butchersimpQ, mtd} =
      {ButcherRowSum, ButcherSimplify, RungeKuttaMethod} /. {opts} /.
         Options[RungeKuttaOrderConditions];
    datatypes = {
      MemberQ[{True, False}, rowsumQ],
      butchersimpQ===False ||
		(IntegerQ[s] && TrueQ[Positive[s]] && butchersimpQ===True),
      MemberQ[{Explicit, DiagonallyImplicit, Implicit}, mtd]};
    If[Apply[And, MapThread[RKOCtest, {datatypes, RKOCmessages}]],
      {rowsumQ, butchersimpQ, mtd},
      $Failed
    ]
  ];

RungeKuttaOrderConditions[p_Integer?Positive,
    s:(_Symbol|_Integer?Positive), opts___?OptionQ]:=
  Module[{ans},
    ans /; ((ans = RKOCOptionsTest[s,opts]) =!= $Failed) &&
              (ans = rkconds[p, s, ans]; True)
  ];

RungeKuttaOrderConditions[p_Integer?Positive, opts___?OptionQ]:=
  Module[{ans},
    ans /; ((ans = RKOCOptionsTest[symb, opts]) =!= $Failed) &&
              (ans = rkconds[p, ans]; True)
  ];

(* Thread Equal through equation lists - avoids problems with
 True and False. *)

eq[a_List,b_List]:= MapThread[eq,{a,b}];
eq[a_,b_]:= Equal[a,b];

(* Order conditions in tensor notation. *)

rkconds[p_, {rowsumQ_ , False, _}]:=
  With[{trees = ButcherTrees[p]},
    If[rowsumQ, Append[#, {a.e==c}]&, (* else *) #& ] @
      eq[ ButcherPhi[trees], 1/ButcherGamma[trees] ]
  ];

(* Numeric stages. *)

rkconds[p_, s_Integer?Positive, {rowsumQ_ , False, mtd_}]:=
  With[{trees = ButcherTrees[p]},
    If[rowsumQ,
	Append[#, First[ButcherRowConditions[1,s,RungeKuttaMethod->mtd]]]&,
	#& ] @
      eq[ ButcherPhi[trees,s,RungeKuttaMethod->mtd], 1/ButcherGamma[trees] ]
  ];

(* Order conditions in Butcher's index notation. *)

rkconds[p_, s_Symbol, {rowsumQ_ , False, mtd_}]:=
  With[{trees = ButcherTrees[p]},
    If[rowsumQ,
      Append[#, ButcherRowConditions[1,s,RungeKuttaMethod->mtd]]&,
      #& (* else *)
    ] @ eq[ ButcherPhi[trees,s,RungeKuttaMethod->mtd], 1/ButcherGamma[trees] ]
  ];


(* Implicit row simplifying conditions A(eta). *)

Options[ButcherRowConditions] = {RungeKuttaMethod :> $RungeKuttaMethod};

ButcherRowConditions[p_Integer?Positive, s:(_Symbol|_Integer?Positive),
    opts___?OptionQ]:=
  Module[{mtd},
    mtd = RungeKuttaMethod /. {opts} /. Options[ButcherRowConditions];
    Table[rowsimpc[j,s,mtd],{j,p}] /; CheckMethod[mtd,ButcherRowConditions]
  ];

rindx[i_,j_,_,Explicit]:= {i,2,j-1};
rindx[i_,j_,_,DiagonallyImplicit]:= {i,j};
rindx[i_,_,s_,Implicit]:= {i,s};

ritab[i_,s_,Explicit]:= {i,2,s};
ritab[i_,s_,_]:= {i,s};

rowsimpc[1, s_, Explicit]:=
  SymbolicTable[
    Sum[Evaluate[a[$i[1],$i[2]],{$i[2],$i[1]-1}]] == c[$i[1]],
  {$i[1],2,s}];

rowsimpc[p_, s_, mtd_]:=
  SymbolicTable[
    Sum[Evaluate[a[$i[1],$i[2]] c[$i[2]]^(p-1),
      rindx[$i[2],$i[1],s,mtd]]] == c[$i[1]]^p/p,
  ritab[$i[1],s,mtd]];


ButcherRowConditions[p_Integer?Positive, opts___?OptionQ]:= Array[brc,p];

brc[1] = a.e == c;
brc[q_]:= brc[q] = a.(c^(q-1)) == c^q/q;


(* Implicit column simplifying conditions D(xi). *)

Options[ButcherColumnConditions] = {RungeKuttaMethod:>$RungeKuttaMethod};

ButcherColumnConditions[p_Integer?Positive, s:(_Symbol|_Integer?Positive),
    opts___?OptionQ]:=
  Module[{mtd},
    mtd = RungeKuttaMethod /. {opts} /. Options[ButcherColumnConditions];
    Table[colsimpc[j,s,mtd],{j,p}] /; CheckMethod[mtd,ButcherColumnConditions]
  ];

cindx[i_,j_,s_,Explicit]:= {i,j+1,s};
cindx[i_,j_,s_,DiagonallyImplicit]:= {i,j,s};
cindx[i_,_,s_,Implicit]:= {i,s};

colsimpc[p_, s_, Explicit]:=
  SymbolicTable[
    Sum[
      Evaluate[a[$i[1],$i[2]] b[$i[1]] c[$i[1]]^(p-1),
        cindx[$i[1],$i[2],s,Explicit]]] == 
          (Apply[If[#1==1,#2,#2 #3]&,{$i[2],b[$i[2]]/p,1-c[$i[2]]^p}])
  ,{$i[2],s}];

colsimpc[p_, s_, mtd_]:=
  SymbolicTable[
    Sum[
      Evaluate[a[$i[1],$i[2]] b[$i[1]] c[$i[1]]^(p-1),
        cindx[$i[1],$i[2],s,mtd]]] == b[$i[2]] (1-c[$i[2]]^p)/p
  ,{$i[2],s}];

ButcherColumnConditions[p_Integer?Positive, opts___?OptionQ]:= Array[bcc,p];

bcc[1] = b.a == b (e-c);
bcc[q_]:= bcc[q] = (b c^(q-1)).a == b (e-c^q)/q;


(* Quadrature conditions (trees of type B(p)). *)

Options[ButcherQuadratureConditions] = {RungeKuttaMethod:>$RungeKuttaMethod};

ButcherQuadratureConditions[p_Integer?Positive,
    s:(_Integer?Positive|_Symbol), opts___?OptionQ]:=
  Module[{mtd},
    mtd = RungeKuttaMethod /. {opts} /. Options[ButcherQuadratureConditions];
    eq[ ButcherPhi[#,s,RungeKuttaMethod->mtd], 1/ButcherGamma[#] ]& @
      Join[{{f}},Table[{f[f^(i-1)]},{i,2,p}]] /;
        CheckMethod[mtd,ButcherQuadratureConditions]
  ];

ButcherQuadratureConditions[p_Integer?Positive, opts___?OptionQ]:=
  eq[ ButcherPhi[#], 1/ButcherGamma[#] ]& @
    Join[{{f}},Table[{f[f^(i-1)]},{i,2,p}]];

(* Simplifying conditions for implicit schemes. *)

rkconds[p_, s_Integer?Positive,
		{rowsumQ_, True, mtd:(DiagonallyImplicit|Implicit)}]:=
  Module[{eta, xi},
    eta = Max[1,Quotient[p-1,2]];
    xi = Max[0,p-eta-1];
    Join[
      ButcherQuadratureConditions[p,s,RungeKuttaMethod->mtd],
      If[rowsumQ, #&, Drop[#,1]&] @
	    ButcherRowConditions[eta,s,RungeKuttaMethod->mtd],
      If[xi==0,{},ButcherColumnConditions[xi,s,RungeKuttaMethod->mtd]]
    ]
  ];


(* Derive trees not reduced by simplifying conditions. *)

VanishTree[t_,t_]:=t;
VanishTree[__]:=Sequence[];

SetAttributes[isimptree,Listable];

(* Trees of type B(p). *)

isimptree[f|f[f^x_.],__]:= Sequence[];

(* One-leg trees, reduced by D(1). *)

isimptree[t:f[f[x_]],_,xi_]:= VanishTree[t] /; xi>=1;

(* D(xi). *)

isimptree[t:f[f^x_. f[y_]],_,xi_]:= VanishTree[t] /; x<=xi-1;

(* A(eta). *)

isimptree[t:f[subt_],eta_,_]:=
  VanishTree[t,f[ subt /. f[f^x_.]:>f^(x+1) /; x<=eta-1 ]];

ButcherTreeSimplify[p_Integer?Positive,eta:(0|_Integer?Positive),
    xi:(0|_Integer?Positive)]:=
  DeleteCases[ isimptree[ ButcherTrees[p] , eta, xi ], {}];


(* Explicit row simplifying conditions A(eta). *)

(* Include explicit row-sum conditions. *)

ExplicitRowSimplify[s_, True, eta_]:=
  Join[
    ButcherRowConditions[1,s,RungeKuttaMethod->Explicit],
    ExplicitRowSimplify[s,False,eta]
  ];

ExplicitRowSimplify[1, False, _]:= Sequence[]; (* One stage case. *)

(* No additional row conditions for p<=4. *)

ExplicitRowSimplify[_, False, 1]:= Sequence[];

(* Special cases for explicit schemes of order p>=5, Butcher's equations
 330e and 330f (b[2]==0 and remaining 2nd row conditions redundant). *)

rowsimpc[i_,q_]:= Sum[a[i,j] c[j]^(q-1), {j, 2, i-1}] == c[i]^q/q;

ExplicitRowSimplify[s_, False, eta_]:=
  Insert[Table[rowsimpc[i,q],{q,2,eta},{i,3,s}],b[2]==0,{1,1}];


(* Explicit column simplifying conditions D(1).
 The first relation s=1 is neglected, since it is implied by imposing
 the row sum conditions, the order one and two trees and the remaining
 column conditions (e.g. Hairer, Norsett and Wanner's proof in lemma 1.4).
 The second relation can be neglected similarly using b[2]=0 for order p>=5.
 Higher order row simplifying conditions (xi>1) are not considered for
 explicit methods since they contradict order p>=1 (Butcher section 331). *)

dropecs[_,  1, _]:= Sequence[];
dropecs[3|4, _,  conds_]:= {Drop[Flatten[conds],1]};
dropecs[_,  2, _]:= Sequence[];
dropecs[_,  s_,  conds_]:= {Drop[Flatten[conds],2] /. b[2]->0};

ExplicitColumnSimplify[p_,s_]:=
  dropecs[p,s,ButcherColumnConditions[1,s,RungeKuttaMethod->Explicit]];


(* Additional trees not reduced by simplifying conditions for
 explicit schemes (Hairer, Norsett and Wanner II.6). *)

SetAttributes[esimptree, Listable];

(* Trees of type B(p) are already defined by ButcherQuadratureConditions. *)

esimptree[f[f^x_.],_]:= Sequence[];

(* One-leg trees are reduced by D(1). *)

esimptree[f[f[x_]],_]:= Sequence[];

(* Reduce trees with one vertex leading to at most
 (eta-1) terminal vertices. *)

esimptree[t:f[f[f^x_.]^y_. z_.],eta_]:= VanishTree[t] /; x<=(eta-1);

(* All remaining trees are not reduced. *)

esimptree[t_,_]:= t;

AdditionalConditions[p_,s_,eta_]:=
  addconds[ DeleteCases[ esimptree[Drop[ButcherTrees[p],3],eta], {} ], s ];

addconds[{}, _]:= {};
addconds[trees_,s_]:=
  eq[ ButcherPhi[trees,s,RungeKuttaMethod->Explicit], 1/ButcherGamma[trees] ];

(* Simplifying conditions for explicit schemes. *)

(* Use standard order conditions for p=1,2. *)

rkconds[p:(1|2), s_Integer?Positive, {rowsumQ_, True, Explicit}]:=
  rkconds[p,s,{rowsumQ,False,Explicit}];

rkconds[p_, s_Integer?Positive, {rowsumQ_, True, Explicit}]:=
  With[{eta = Quotient[p-1,2]},
    Join[
      If[p>=5, (# /. b[2]->0)&, #& ] @
        ButcherQuadratureConditions[p,s,RungeKuttaMethod->Explicit],
      ExplicitRowSimplify[s,rowsumQ,eta],
      ExplicitColumnSimplify[p,s],
      If[p>=5, (# /. b[2]->0)&, #& ] @ AdditionalConditions[p,s,eta]
    ]
  ];

(**** Tree plotting ****)

(* Find the list of subtrees (one level up) of a tree. *)

SetAttributes[splitpower,Listable];
splitpower[subt_^n_.]:= Table[subt,{n}];
split[f[subt_Times]]:= Flatten[splitpower[Apply[List,subt]]];
split[f[subt_]]:= splitpower[subt];

buildtree[f, _, _, x_, y_]:= {Point[{x, y}]};
buildtree[t_f, h_, w_, x_, y_]:=
  Module[{nx, ny=y+h, nw, sp = split[t]},
    nw = w/Length[sp];
    nx = x - w/2 - nw/2;
    Flatten[
      {Point[{x, y}],
       Map[ {Line[{{x,y}, {nx+=nw,ny}}], buildtree[#, h, nw, nx, ny]}&, sp]}
    ]
  ];

(* Encircle root. *)

root[0|0., _]:= Point[{0,0}];
root[r_, ar_]:= {Point[{0,0}], Circle[{0,0}, {r, r/ar}]};

Options[ButcherPlot] =
  {ButcherPlotColumns -> Automatic, ButcherPlotNodeSize -> 1,
   ButcherPlotRootSize -> 1, ButcherPlotLabel->{}};

ButcherPlot::opts =
"The option `1` in ButcherPlot did not evaluate to `2`"

BTPmessages = {
  {"ButcherPlotNodeSize","a positive number"},
  {"ButcherPlotRootSize","zero or a positive number"},
  {"ButcherPlotColumns","a positive integer or Automatic"},
  {"ButcherPlotLabel","a list of plot labels"}};

BTPtest =
  If[#1,
    True,
    Message[ButcherPlot::opts,Apply[Sequence,#2]]; False
  ]&;

BTPOptionsTest[ff_,opts___]:=
  Module[{ns, rs, col, lbl, datatypes},
    {ns,rs,col,lbl} = {ButcherPlotNodeSize,ButcherPlotRootSize,
        ButcherPlotColumns,ButcherPlotLabel} /.
            {opts} /. Options[ButcherPlot];
    lbl = Flatten[lbl];
    datatypes = {
      TrueQ[Positive[ns]],
      MatchQ[rs,0|0.|_?Positive],
      MatchQ[col, Automatic | (_Integer?Positive)],
      MatchQ[lbl,{}|_?(VectorQ[#]&&Length[#]===Length[ff]&)]};
    If[Apply[And, MapThread[BTPtest, {datatypes, BTPmessages}]],
      {ns,rs,col,lbl},
      $Failed
    ]
  ];

ButcherPlot[forest_?ButcherTreeQ, opts___?OptionQ]:=
  Module[{ans,ff = Flatten[{forest}]},
    ans /; ((ans = BTPOptionsTest[ff,opts]) =!= $Failed &&
                (ans = BTPlot[ff,ans]; True))
  ];

BTPlot[ff_, {ns_, rs_, col_, lbl_}]:=
  Block[{$RecursionLimit=Infinity},
    Module[{ar, br, h, lenff, m, n, pf},
      lenff = Length[ff];
      h = Max[1, Map[height,ff] - 1];
      n = If[col === Automatic, Ceiling[N[Sqrt[Sqrt[h] lenff]]], col];
      m = Ceiling[lenff/n];
      ar = Sqrt[N[h]];
      br = Max[rs, rs ns] 0.03 Sqrt[ar];
      pf = MapThread[
             Graphics[
               {PointSize[ns (.0001 m n)^(.4)], root[br ar, ar],
                buildtree[#1, 1/h, .95, 0, 0]},
               PlotRange -> {{-.52-br, .52+br}, {-.02-br, 1.02+br}},
               AspectRatio -> ar, PlotLabel -> #2
             ]&,
             {ff, If[lbl==={}, Table[None,{lenff}], lbl]}
           ];
      h = Partition[pf, n];
      If[Length[h] < m, AppendTo[h, Drop[pf, n(m-1)]]];
      Show[GraphicsArray[h], AspectRatio -> 1]
    ]
  ];

End[ ];		(* End `Private` Context. *)

SetAttributes[{ButcherAlpha, ButcherBeta, ButcherBetaBar,
ButcherColumnConditions, ButcherGamma, ButcherHeight, ButcherOrder,
ButcherPhi, ButcherPlot, ButcherPrincipalError,
ButcherQuadratureConditions, ButcherRowConditions, ButcherRowSum,
ButcherSigma, ButcherSimplify, ButcherTrees, ButcherTreeCount,
ButcherTreeQ, ButcherTreeSimplify, ButcherWidth,
RungeKuttaOrderConditions, SymbolicTable}, ReadProtected];

Protect[a, b, c, e, f, ButcherAlpha, ButcherBeta, ButcherBetaBar,
ButcherColumnConditions, ButcherGamma, ButcherHeight, ButcherOrder,
ButcherPhi, ButcherPlot, ButcherPlotColumns, ButcherPlotLabel,
ButcherPlotNodeSize, ButcherPlotRootSize, ButcherPrincipalError,
ButcherQuadratureConditions, ButcherRowConditions, ButcherRowSum,
ButcherSigma, ButcherSimplify, ButcherTrees, ButcherTreeCount,
ButcherTreeQ, ButcherTreeSimplify, ButcherWidth, Explicit,
DiagonallyImplicit, Implicit, RungeKuttaMethod, RungeKuttaOrderConditions,
SymbolicTable];

EndPackage[ ];		(* End package Context. *)
