Runge-Kutta methods, trees, and Mathematica
Author: Folkmar Bornemann

A Mathematica-notebook accompanying e-print arXiv:math.NA/0211049 
Last change: November 4, 2002, 15:00

5. Generating order conditions with Mathematica.

Attributes[Tree] = {Orderless, 
          Listable}; 
ToTree[f_] := Tree @@ ToTree /@ f;  
    ToTree[f_Symbol] := Tree[];
\(Format[β_Tree] := 
      StringReplace[
        ToString[
          List @@ β /. {} -> "⊙"], {"{" -> "[", "}" -> "]"}]\)

(Subscript[β, 1]) = 
      ToTree[(f^′)′[(f^′)′[
            f, (f^′)[f]], f]]; (Subscript[β, 2]) = 
      ToTree[(f^′)′[
          f, (f^′)′[(f^′)[f], 
            f]]];
If[(Subscript[β, 1]) == (Subscript[β, 2]), (Subscript[β, 1]), , ]

TreeOrder[β_] := 
        1 + Plus @@ TreeOrder /@ β;

β = 
      ToTree[(f^′)′[(f^′)′′[(f^′)[f], (f^′)[f], f], f]];  TreeOrder[β]

TreeFactorial[β_] := 
      TreeOrder[β] Times @@ TreeFactorial /@ β

TreeFactorial[β]

TreeAlpha[β_] := 
      Length [Permutations[β]]/Length[β]! 
        Times @@ TreeAlpha /@ β

TreeAlpha[β]

TreeA[β_, 
        n_] := (
vars = {i, j, k, l, m, p, q, r, u, v, 
            w}; 
Times @@ (Sum[(Subscript[a, vars[[n]], 
                          vars[[n + 1]]]) TreeA[#, n + 1] // 
                    Evaluate, 
{vars[[n + 1]], s} // 
                    Evaluate] & /@ β))
TreeOrderCondition[β_] := 
      Sum[(Subscript[b, i]) TreeA[β, 1] // Evaluate, {i, s}] == 
        1/TreeFactorial[β]

TreeOrderCondition[β]

AddLeaf[β_] := AddLeaf[β] = 
        Fold[Union, {Prepend[β, 
              Tree[]]}, 
ReplacePart[β, 
                AddLeaf[β[[#]]], #] & /@ 
            Range[Length[β]]]
Trees[order_] := 
      NestList[Union @@ AddLeaf /@ # &, {Tree[]}, order - 1]

Trees[4]

Length /@ Trees[10]

Plus @@ %

explicit = {(Subscript[a, i_, j_]) :> 0 /; i ≤ j, 
          (Subscript[c, 1]) -> 0};
OrderConditions[order_, 
        stages_] := (
autonomous = 
          Table[Sum[(Subscript[a, i, j]), {j, stages}] == (Subscript[c, i]), {i, 
              stages}]; 
{(TreeOrderCondition /@ Union @@ Trees[order]) /. 
                s -> stages /. 
ToRules[
                And @@ autonomous], autonomous} /. explicit)

First[OrderConditions[4, 4]]

6. Examples of usage.

RungeKuttaMethod[p_, s_, 
        pre_] := (
rkTemplate = {Table[
                (Subscript[a, i, j]), {i, s}, {j, s}], 
              Table[(Subscript[b, i]), {j, 1}, {i, s}], 
Table[
                (Subscript[c, i]), {j, 1}, {i, s}]} /. 
            explicit; 
conditions = 
          pre ⋃ 
            Union @@ OrderConditions[p, s]; 
solveVars = 
          Complement[Flatten[rkTemplate], {0}]; 
sol = 
          Solve[conditions, 
            solveVars]; 
Thread[{A, b, c} == 
                  MatrixForm /@ rkTemplate /. #] & /@ sol)

RungeKuttaMethod[4, 
      4, {(Subscript[c, 2]) == (Subscript[c, 3]), (Subscript[b, 2]) == (Subscript[b, 3])}]

pre = {(Subscript[c, 2]) == u, (Subscript[c, 3]) == 1/4, (Subscript[c, 4]) == 1/2, 
          (Subscript[c, 5]) == 3/4, (Subscript[c, 6]) == 1, (Subscript[b, 2]) == 0, 
          (Subscript[a, 4, 3]) == v};
First[RungeKuttaMethod[5, 6, pre]]