Source code for pybs.series.Bseries

from fractions import Fraction

from pybs.utils import LinearCombination
from pybs.combinations import empty_tree, Forest
from pybs.unordered_tree import UnorderedTree, leaf


[docs]class BseriesRule(object): """Objects of this class are B-series rules for numerical methods. They treat forests as characters of the Hopf algebra. For other B-series rules use :Class:`VectorfieldRule` or :class:`ForestRule`.""" def __init__(self, arg=None): if arg is None: self._call = lambda x: 0 elif isinstance(arg, LinearCombination): self._call = lambda tree: arg[tree] # TODO: Check that there are no non-trees. elif callable(arg): self._call = arg def __call__(self, arg): if isinstance(arg, UnorderedTree) or arg == empty_tree: return self._call(arg) elif isinstance(arg, Forest): result = 1 for tree, multiplicity in arg.items(): result *= self._call(tree) ** multiplicity # TODO: Use reduce() or something? return result # else: # return 0 elif isinstance(arg, LinearCombination): result = 0 for elem, multiplicity in arg.items(): result += self(elem) * multiplicity return result
[docs]class VectorfieldRule(object): """Objects of this class are B-series rules that correspond to modified equations. That means they act on forests as infinitesimal characters. """ def __init__(self, arg=None): if arg is None: self._call = lambda x: 0 elif isinstance(arg, LinearCombination): self._call = lambda tree: arg[tree] # TODO: Check that this is reasonable. elif callable(arg): self._call = arg def __call__(self, arg): if isinstance(arg, UnorderedTree) or arg == empty_tree: return self._call(arg) elif isinstance(arg, Forest): if arg.number_of_trees() == 1: # TODO: Do nicer. for elem in arg: return self._call(elem) else: return 0 elif isinstance(arg, LinearCombination): result = 0 for elem, multiplicity in arg.items(): result += self(elem) * multiplicity return result
[docs]class ForestRule(object): """General rule with arbitrary behavior for forest. Use this if the two others are unsuitable. """ # Results on forests are not deducable from results on trees. def __init__(self, arg=None): if arg is None: self._call = lambda x: 0 elif isinstance(arg, LinearCombination): self._call = lambda forest: arg[forest] elif callable(arg): self._call = arg def __call__(self, arg): if isinstance(arg, (UnorderedTree, Forest)): return self._call(arg) elif isinstance(arg, LinearCombination): result = 0 for elem, multiplicity in arg.items(): result += self(elem) * multiplicity return result
def _unit(tree): if tree == empty_tree: return 1 else: return 0 def _exact(tree): if tree == empty_tree: return 1 return Fraction(1, tree.density()) def _kahan(tree): if tree == empty_tree: return 1 if tree.is_tall(): return Fraction(1, (2 ** (tree.order()-1)) * tree.symmetry()) else: return 0 # TODO: Necessary?? def _AVF_old(a, tree): if tree == empty_tree: return 1 if tree == leaf: return 1 elif not tree.is_binary(): return 0 else: if tree.number_of_children() == 1: return Fraction(_AVF(a, tree.keys()[0]), 2) elif tree.number_of_children() == 2: alpha = Fraction(2*a + 1, 4) if len(tree._ms) == 2: return alpha * _AVF(a, tree.keys()[0]) * \ _AVF(a, tree.keys()[1]) else: return alpha * _AVF(a, tree.keys()[0]) ** 2 def _AVF(tree): "According to `Energy-Preserving Runge-Kutta Methods`, Celledoni et al." if tree == empty_tree: return 1 if tree == leaf: return 1 else: result = Fraction(1, tree.number_of_children() + 1) for child_tree, multiplicity in tree.items(): result *= _AVF(child_tree) ** multiplicity return result def _unit_field(tree): if tree == leaf: return 1 return 0 exponential = BseriesRule(_exact) unit = BseriesRule(_unit) unit_field = VectorfieldRule(_unit_field) AVF = BseriesRule(_AVF)