(*
The Package!
*)
BeginPackage["RungeKutta`"]
(*
First define the usage messages.
*)
RungeKutta::usage =
"This packages is designed to compute RungeKutta order conditions. The main functions are:\n\nRKConditions which computes the order conditions along with conditions for row and column simplifying assumptions. \n\nShowTree which displays the trees corresponding to the conditions.\n\nOrder, Density, Symmetry, Height, Width, Alpha, Beta and BetaBar give the corresponding value of a certain tree. See Butcher: The Numerical Analysis of Ordinary Differential Equations (Runga-Kutta and general linear methods), John Wiley & Sons, New York, 1987.\n\nReferences to trees is primarily done via the tree's number (that is it's position with the list of conditions/trees). To aid the search of a certain tree the commands TreeToNumber, NumberToTree and FindDensity have been developed."
RKConditions::usage =
"RKConditions[p, c, a, b, opts] creates the conditions that the Runge-Kutta method with Butcher array\n
\t c|a \n
\t --- \n
\t |b \n
must meet to be of order p. The conditions are named CondName[1] through CondName[?]. Continuous order conditions are possible using the option Continuous. a and c are assumed to be lists of the same length and b should be a square matrix with dimension the same as a and c. a, b and c need not be defined before evaluating RKConditions. In this case only the conditions are evaluated and are given in terms of multiplications and dot products. Row and column simplifications are also computed and named. The quadrature conditions, as well as being part of the order conditions, are given special names.\nThe options are: CondName, RowName, ColumnName, QuadName, Continuous, EndPointOrder, RowOrder, ColumnOrder."
CondName::usage =
"The conditions will be named CondName[1] through CondName[?]. The default name is cond."
QuadName::usage =
"The quadrature conditions will be named QuadName[1] through QuadName[p] where p is the order requested in RKConditions. The default name is qcond."
RowName::usage =
"The row simplifying assumptions will be named RowName[1] through RowName[RowOrder]. The default name is rcond."
ColumnName::usage =
"The column simplifying assumptions will be named ColumnName[1] through ColumnName[ColumnOrder]. The default name is ccond."
Continuous::usage =
"If set to anything other than False then the value is assumed to be the variable used for continuous order conditions. The default is False."
EndPointOrder::usage =
"specifies the order that continuous order conditions should attain at their end points. The default value is p+1, where p is the order requested in RKConditions."
RowOrder::usage =
"specifies the order of the row simplifying assumptions. The default is p, where p is the order requested in RKConditions."
ColumnOrder::usage =
"specifies the order of the column simplifying assumptions. The default is p, where p is the order requested in RKConditions."
Symmetry::usage =
"Symmetry[treenumber] gives the symmetry of the given tree."
Alpha::usage =
"Alpha[treenumber] gives the number of ways of labelling the given tree with a totally ordered set V (with #V = Order[treenumber]) in such a way that if (m,n) is an edge then mLength[#2],
True,
If[Length[#1]==Length[#2], #1<#2, False]
]&
];
Do[treelist =
Join[treelist, PartToTrees[ ordpart[[i]] ] ],
{i, Length[ordpart]}];
startorder = Append[startorder, Length[treelist]+1],
{ord, currentorder+1, p}];
treeno = startorder[[currentorder+1]];
Do[
Do[ den = 1;
list = treelist[[treeno]];
Do[ den *= densitylist[[ list[[j]] ]], {j,Length[list]}];
orderlist = Append[orderlist, ord];
densitylist = Append[densitylist, den*ord];
treeno++,
{startorder[[ord+1]] - startorder[[ord]]}],
{ord, currentorder+1, Length[startorder]-1}]
];
(*
ExtendTreeList[treeno] extends the list of trees by repeated calling MakeTreeList until the length of treelist is at least treeno.
*)
ExtendTreeList[treeno_] :=
While[ treeno >= startorder[[-1]],
MakeTreeList[Length[startorder]]
]
(*
RKConditions[p,a,b,w,opts] is the main function in the package and produces the conditions for order p. First MakeTreeList is called. Then the conditions are constructed recursively as described in Butcher section 144 and somewhat more explicitly described in Hairer/Norsett/Wanner Problem II.2.1. We follow this construction by defining the condition for order 1 and then using this recursively for higher orders. The conditions themselves are first defined using dot products and multiplications of a and b. When all the conditions are constructed the dependence on w is then added with a dot product. The conditions may eeasily be viewed by executing RKConditions with a, b and w undefined. This way they are expressed in terms of the dot products and multiplications. In the case when a,b and w are undefined it is necessary to leave some of the conditions in unevaluated for (eg condition 1 reads cond[1] = Sum[w[[i]], {i,Length[w]}]) the reason is the dependence on the length of w. Once the conditions are derived the quadrature conditions, row simplification assumptions and column simplification assumptions and determined, agian with dot products and multiplications. name is used for the naming of the conditions, qname for the quadrature conditions, rname for the row simplification conditions and cname for the column simplifications. var is defined to be either 1 (in the case of non-continuous conditions) or the variable for continuous conditions (in the case of a continuous method). In the case of a continuous method it is possible that the end point order of the method (which should have the variable replaced by 1) hence the code with var->1.
*)
Options[RKConditions] =
{ CondName->Global`cond,
QuadName->Global`qcond,
RowName->Global`rcond,
ColumnName->Global`ccond,
Continuous->False
};
RKConditions[p_,a_,b_,w_,opts___] :=
With[{ name = (CondName /. {opts} /. Options[RKConditions]),
qname = (QuadName /. {opts} /. Options[RKConditions]),
rname = (RowName /. {opts} /. Options[RKConditions]),
cname = (ColumnName /. {opts} /. Options[RKConditions]),
cts = (Continuous /. {opts} /. Options[RKConditions]),
roword = (RowOrder /. {opts} /. RowOrder->p),
colord = (ColumnOrder /. {opts} /. ColumnOrder->p)
},
Module[{partitions, numtrees, tree, cond, epord, var},
If[ cts =!= False,
var = cts;
epord = (EndPointOrder /. {opts} /. EndPointOrder->p+1),
var = 1;
epord = p
];
Clear[name,qname,rname,cname];
MakeTreeList[Max[p,epord]];
numtrees = Length[treelist];
name[1] = ones;
Do[ cond = 1;
tree = treelist[[i]];
Do[ cond *= If[ tree[[j]]==1, a, b.name[tree[[j]]] ],
{j, Length[tree]}];
name[i] = cond,
{i, numtrees}];
Do[name[i] =
w.name[i] == var^orderlist[[i]]/densitylist[[i]],
{i,numtrees}];
Do[ name[i] = (name[i] /. var->1), {i, startorder[[p+1]], numtrees}];
Attributes[name] = {Listable};
If[VectorQ[a] && MatrixQ[b] && VectorQ[w],
name[1] = Sum[w[[i]], {i,Length[w]}] == var;
If[qname =!= None,
qname[i_] := name[startorder[[i]]];
Attributes[qname] = {Listable}
];
If[rname =!= None,
rname[1] = Thread[
Table[Sum[b[[i,j]], {j,Length[b]}],
{i,Length[b]}] - a == 0
];
Do[ rname[i] = Thread[b.(a^(i-1)) - (a^i)/i == 0], {i,2,roword}];
rname[i_,j_] := rname[i][[j]];
Attributes[rname] = {Listable}
];
If[cname =!= None,
cname[1] = Thread[w.b -
w * (Table[1, {j,Length[w]}] - a) == 0];
Do[ cname[i] = Thread[(w*a^(i-1)) . b -
w * (Table[1, {j,Length[w]}] - a^i) / i == 0],
{i, 2, colord}];
cname[i_,j_] := cname[i][[j]];
Attributes[cname] = {Listable};
],
name[1] = Unevaluated[Sum[w[[Global`i]], {Global`i,Length[w]}]] == var;
]
]
];
(*
FindDensity[density] produces a list of all those trees currently defined that have density density.
*)
FindDensity[density_] :=
Module[{trees={}},
Do[
If[densitylist[[i]]==density,
trees = Append[trees, i]
],
{i,Length[densitylist]}];
trees
]
(*
TreeToNumber[tree] converts a tree as described in usage (different from the version in treelist) to the number of that tree (ie its position in treelist). The definition is a self-calling recursive one. The idea is to take each subtree and find its number then, with the list of numbers that results, search treelist for this tree. (The search is reduced by first determining the order of the tree).
*)
TreeToNumber[tree_List] :=
Module[{joins={}, treeno, order=1},
Do[joins =
Join[joins,
If[NumberQ[tree[[i]]],
Table[1, {tree[[i]]}],
{TreeToNumber[tree[[i]]]}
]
],
{i, Length[tree]}];
joins = Sort[joins, Greater];
Do[order += orderlist[[joins[[i]]]],
{i,Length[joins]}];
MakeTreeList[order];
treeno = startorder[[order]];
While[treelist[[treeno]] != joins, treeno++];
treeno]
(*
NumberToTree[treeno] converts a tree number to a description of the tree. The description is given in the usage notes. The construction is basically done with ConstructTree which is self-calling. From the list of trees to join (coming from treelist) one converts each of these tree numbers to its tree (hence the self-calling mechanism) and then joins each of these subtree constructions. The only thing to look out for is that any 1's at the end of the list of tree numbers to join (which always appear at the end because of the non-increasing order) should be grouped together.
*)
NumberToTree[treeno_] :=
Module[{},
ExtendTreeList[treeno];
ConstructTree[treelist[[treeno]]]
]
ConstructTree[tree_List] :=
Module[{treejoin={}, i=1},
Do[
If[tree[[i]]==1,
treejoin = Append[treejoin, Length[tree]-i+1];
Break[],
treejoin = Append[treejoin,
ConstructTree[treelist[[tree[[i]] ]] ] ]
],
{i,Length[tree]}];
treejoin
]
(*
Order[treeno] first extends the list of trees and then uses orderlist to give the order of a tree.
*)
Attributes[Order] = {Listable};
Order[treeno_] :=
Module[{},
ExtendTreeList[treeno];
orderlist[[treeno]]
]
(*
Density[treeno] first extends the list of trees and then uses densitylist to give the order of a tree.
*)
Attributes[Density] = {Listable};
Density[treeno_] :=
Module[{},
ExtendTreeList[treeno];
densitylist[[treeno]]
]
(*
Symmetry[treeno] first extends the list of trees and then uses the recursive definition of symmetry given in Butcher formula 144b&c to determine the symmetry of a tree.
*)
Attributes[Symmetry] = {Listable};
Symmetry[1] = 1;
Symmetry[treeno_] :=
Module[{tree, subtree, subtreepower, symmetry=1},
ExtendTreeList[treeno];
tree = treelist[[treeno]];
subtreepower = 1;
Do[ subtree = tree[[i]];
If[ i != Length[tree] && subtree == tree[[i+1]],
subtreepower++,
symmetry *= Factorial[subtreepower] * Symmetry[subtree]^subtreepower;
subtreepower = 1
],
{i, Length[tree]}];
symmetry
]
(*
Height[treeno] first extends the list of trees and then uses the recursive definition of height given in Butcher section 144 to determine the height of a tree.
*)
Attributes[Height] = {Listable};
Height[treeno_] :=
Module[{},
ExtendTreeList[treeno];
FindHeight[treeno]
]
FindHeight[1] = 1;
FindHeight[treeno_] :=
Max[Table[FindHeight[treelist[[treeno,i]]]+1, {i,Length[treelist[[treeno]]]}]]
(*
Width[treeno] first extends the list of trees and then uses the recursive definition of width given in Butcher section 144 to determine the width of a tree. This is the widthbar version of Butcher.
*)
Attributes[Width] = {Listable};
Width[treeno_] :=
Module[{},
ExtendTreeList[treeno];
FindWidth[treeno]
]
FindWidth[1] = 1;
FindWidth[treeno_] :=
Sum[FindWidth[treelist[[treeno,i]]], {i,Length[treelist[[treeno]]]}]
(*
Alpha[treeno] first extends the list of trees and then uses the definition of alpha given in Butcher theorem 145E to determine alpha for a tree.
*)
Attributes[Alpha] = {Listable};
Alpha[treeno_] :=
orderlist[[treeno]]! / (densitylist[[treeno]] * Symmetry[treeno])
(*
Beta[treeno] first extends the list of trees and then uses the definition of beta given in Butcher theorem 145B to determine beta for a tree.
Beta[treeno, n] uses the theorem 145C in Butcher.
*)
Attributes[Beta] = {Listable};
Beta[treeno_] :=
(Order[treeno] - 1)! / Symmetry[treeno]
Beta[treeno_, n_] :=
Module[{r, w, s},
r = Order[treeno];
w = If[treeno==1, 0, Width[treeno]];
s = Symmetry[treeno];
n! (r-1-w)! / (s (n-w)! (r-1-n)!)
]
(*
BetaBar[treeno] first extends the list of trees and then uses the definition of betabar given in Butcher theorem 145B to determine betabar for a tree.
BetaBar[treeno, n] uses the theorem 145C in Butcher.
*)
Attributes[BetaBar] = {Listable};
BetaBar[treeno_] :=
Order[treeno]! / Symmetry[treeno]
BetaBar[treeno_, n_] :=
Module[{r, w, s},
r = Order[treeno];
w = Width[treeno];
s = Symmetry[treeno];
n! (r-w)! / (s (n-w)! (r-n)!)
]
(*
NumberOfTrees[p] produces two lists. The i-th element of the first list gives the number of trees/conditions of order i, whilst the i-th element of the second list gives the number of trees/conditions of order at most i. The definition follows the construction in Butcher theorem 145A.
*)
NumberOfTrees[p_] :=
Module[{a, series},
series =
Sum[a[n] x^(n-1), {n,p}] -
Normal[Series[Product[(1-x^n)^-a[n], {n,p}], {x,0,p}]];
Do[ a[n] = Solve[Coefficient[series, x, n-1]==0, a[n]][[1,1,2]], {n,p}];
{Table[a[n], {n,p}],
Table[Sum[a[i], {i, n}], {n,p}]}
]
(*
MakePointsAndLines[tree, root, width] given a tree (that is a list of tres to join as in treelist), coordinates of its root and the width that the tree may occupy this module determines all the points and lines that are needed to draw that tree. The definition is a self-calling recursive one in which the root of a tree is appended to the end of the list of points and then each of the bramches sprouting from the root gives a line which is added to the list of lines. The branches from each root are speard evenly over the width interval. A result of the construction is that trees never actually reach their horizontal bounds but bushy trees come close. The main idea is that to find the points and lines for a tree one calls MakePointsAndLines with the tree to be converted, with root set to (0,0) and with a width of 1 (so that all trees will be bounded between -1 and +1 horizontally). The various vertical levels of a tree are set at increments of 1 so that a tree of height h will be bounded vertically by 0 and h-1.
*)
MakePointsAndLines[tree_, root_, width_] :=
Module[{branches, newroot, newwidth},
pointsandlines[[1]] = Append[pointsandlines[[1]], Point[root]];
branches = Length[tree];
If[ branches != 0,
newwidth = width/branches;
Do[ newroot = root + {newwidth*(2i-1)-width, 1};
pointsandlines[[2]] =
Append[pointsandlines[[2]], Line[{root, newroot}]];
MakePointsAndLines[treelist[[ tree[[i]] ]], newroot, newwidth],
{i, Length[tree]}]
];
pointsandlines
]
(*
MakeTreePic[treeno, height, opts] constructs a picture of the given tree and includes some text for that tree (probably things like tree number, density etc.) After calling MakePointsAndLines the graphical text for the picture is constcuted using the name provided by ShowText (this should be a function of the form name[treenumber, height, fontsize] and should prooduce a graphics output. After this is added to the picture, two invisible points are then added. These points ensure that PlotRange->All doesn't rescale the picture. The points are at the upper-right and lower-left corners. Picture drawing is suppressed at this point.
*)
MakeTreePic[treeno_, height_, opts___] :=
With[
{fontsize = FontSize /. {opts} /. Options[ShowTree],
nodesize = NodeSize /. {opts} /. Options[ShowTree],
showtext = ShowText /. {opts} /. Options[ShowTree],
aspectratio = AspectRatio /. {opts} /. Options[ShowTree]},
Module[{},
pointsandlines={{},{}};
MakePointsAndLines[treelist[[treeno]], {0,0}, 1];
Show[
Graphics[ {PointSize[nodesize], Flatten[pointsandlines]}],
If[showtext =!= None,
showtext[treeno, height, fontsize],
Graphics[Point[{0,0}]]
],
Graphics[ {GrayLevel[1],
{Point[{1,height}], Point[{-1,0}]}}
],
AspectRatio->aspectratio,
PlotRange->All,
DisplayFunction->Identity,
AspectRatio->aspectratio
]
]
]
(*
NumberAndDensity[treeno, height, fontsize] produces graphical text showing the number of a tree and its density (the two main numbers used in solving Runge-Kutta conditions). The text is placed below the picture.
*)
NumberAndDensity[treeno_, height_, fontsize_] :=
Graphics[
Text[
FontForm[
StringForm["`` g=``", treeno, densitylist[[treeno]]],
{"Symbol", fontsize}],
{0,-0.15}, {0,1}]
]
(*
Everything[treeno, height, fontsize] was intened to be an alternative to NumberAndDensity that shoed the trees number, height, width, density, symmetry, order, alpha, beta and betabar. I haven't resolved the problem of aligning and mixing fonts within a graphical text output. The intention is that this would be the default text shown when only one tree is to bo drawn.
*)
Everything[treeno_, height_, fontsize_] :=
Graphics[
]
(*
ShowTree[treeno_Integer, opt] draws a picture of a single tree.
*)
ShowTree[treeno_Integer, opts___] :=
Module[{},
Options[ShowTree] =
{ NodeSize->.035,
FontSize->18,
ShowText->NumberAndDensity,
AspectRatio->11/8.5
};
Show[MakeTreePic[treeno, Height[treeno]-1, opts],
DisplayFunction->$DisplayFunction,
AspectRatio->(AspectRatio /. {opts} /. Options[ShowTree]),
PlotRange->All
];
];
(*
ShowTree[treeno_List, opt] draws a picture of a list of trees. GraphicsArray is used. First the number of rows is determined with q. A list of lists of zeros is then created and the graphics elements will take the place of the zeros.
*)
ShowTree[list_List, opts___] :=
Module[{pic, height, piclist, columns, q, r, i, j},
Options[ShowTree] =
{ NodeSize->.035,
FontSize->10,
Columns->5,
ShowText->NumberAndDensity,
AspectRatio->11/8.5
};
height = Order[Max[list]];
columns = Columns /. {opts} /. Options[ShowTree];
q = Floor[Length[list]/columns];
r = Length[list] - columns * q;
If[r==0, q--; r=columns];
piclist = Append[Table[0, {q}, {columns}], Table[0, {r}]];
i=1; j=1;
Do[ pic =
MakeTreePic[list[[treeno]], height-1, opts];
piclist[[i,j]] = pic;
If[++j > columns, j=1;i++],
{treeno, Length[list]}]
Show[ GraphicsArray[ piclist,
GraphicsSpacing->0 ],
AspectRatio->(AspectRatio/.{opts}/.Options[ShowTree])*(q+1)/columns ];
];
End[];
EndPackage[];