magnus.expansionterms

expansionterms.py

Generate the terms of the Magnus expansion symbolically, to any order.

The numerical core in magnus.magnus evaluates the Magnus expansion with the coefficients of each commutator group written out explicitly in Python, which is fast but fixes the highest order at whatever was typed in. This module derives those same terms from the recursion itself, in exact rational arithmetic, for any order – so the hard-coded ones can be checked rather than trusted, and so an order beyond the implemented ceiling can still be inspected on paper.

The recursion is the standard Bernoulli-number one [1] (in the \(B_1 = -1/2\) convention):

\[\begin{split}\Omega_1(t) &= \int_0^t A(s)\, ds \\ \Omega_n(t) &= \sum_{j=1}^{n-1} \frac{B_j}{j!} \int_0^t S_n^{(j)}(s)\, ds ,\end{split}\]

where the \(S_n^{(j)}\) are themselves defined recursively,

\[\begin{split}S_n^{(1)} &= [\Omega_{n-1}, A] \\ S_n^{(j)} &= \sum_{m=1}^{n-j} [\Omega_m, S_{n-m}^{(j-1)}] , \qquad 2 \leq j \leq n-1 ,\end{split}\]

so that every term of \(\Omega_n\) is a nested commutator of lower-order \(\Omega_m\) with \(A\), carrying a rational coefficient. Because \(B_j = 0\) for every odd \(j \geq 3\), whole groups drop out: only \(j = 1, 2, 4, 6, \ldots\) contribute.

This is the same form the numerical core implements, deliberately. Writing \(\Omega_n\) instead as time-ordered multiple integrals of \(A\) alone is the other common presentation, but it does not correspond to anything the code evaluates, so it would be of no use for checking the implementation.

The number of terms grows quickly – 1 at order 1, 26 at order 6, 211 at order 8, 1918 at order 10 – which is why the implemented ceiling is a deliberate choice rather than an oversight. See Magnus Expansion Terms to Any Order for the derivation, the expansion printed out, and worked examples.

References

Routine listings

bernoulli - Bernoulli number B_n as an exact Fraction

bernoulli_factor - The coefficient B_j / j! of a commutator group

omega_terms - Terms of Omega_n, as (coefficient, nested-commutator) pairs

magnus_terms - omega_terms for every order up to the one requested

format_term - One term as a readable string

print_magnus_terms - The expansion, printed order by order

count_terms - Number of terms in Omega_n, without building them

Attributes

`Word`
`Term`

Functions

`bernoulli`(→ fractions.Fraction)	Returns the Bernoulli number \(B_n\) as an exact `fractions.Fraction`.
`bernoulli_factor`(→ fractions.Fraction)	Returns \(B_j / j!\), the coefficient multiplying the \(j\)-th commutator
`omega_terms`(→ Tuple[Term, Ellipsis])	Returns the terms of \(\Omega_n\) for `n = order`.
`magnus_terms`(→ Dict[int, Tuple[Term, Ellipsis]])	Returns `omega_terms()` for every order from 1 up to `max_order`.
`count_terms`(→ int)	Returns the number of terms in \(\Omega_n\).
`format_term`(→ str)	Renders one term as a readable string.
`print_magnus_terms`(→ None)	Prints the Magnus expansion, order by order, up to `max_order`.

Module Contents

magnus.expansionterms.Word[source]

magnus.expansionterms.Term[source]

magnus.expansionterms.bernoulli(n: int) → fractions.Fraction[source]

Returns the Bernoulli number \(B_n\) as an exact fractions.Fraction.

Uses the \(B_1 = -1/2\) convention, which is the one the Magnus recursion above is written in, and computes from the defining recursion

\[B_m = -\frac{1}{m+1} \sum_{j=0}^{m-1} \binom{m+1}{j} B_j ,\]

in exact rational arithmetic, so the result is not subject to rounding at any order. Every odd \(B_n\) with \(n \geq 3\) comes out exactly zero.

Added in version 1.0.0.

Parameters:: n (int) – Index of the Bernoulli number; must be >= 0.
Returns:: \(B_n\).
Return type:: fractions.Fraction
Raises:: ValueError – If n is negative.

Examples

import magnus.expansionterms as et

[str(et.bernoulli(n)) for n in range(9)]

['1', '-1/2', '1/6', '0', '-1/30', '0', '1/42', '0', '-1/30']

magnus.expansionterms.bernoulli_factor(j: int) → fractions.Fraction[source]

Returns \(B_j / j!\), the coefficient multiplying the \(j\)-th commutator group in the Magnus recursion.

These are the numbers the numerical core hard-codes: \(B_1/1! = -1/2\), \(B_2/2! = 1/12\), \(B_4/4! = -1/720\), \(B_6/6! = 1/30240\).

Added in version 1.0.0.

Parameters:: j (int) – Index of the commutator group; must be >= 0.
Returns:: \(B_j / j!\).
Return type:: fractions.Fraction

Examples

import magnus.expansionterms as et

{j: str(et.bernoulli_factor(j)) for j in (1, 2, 4, 6)}

{1: '-1/2', 2: '1/12', 4: '-1/720', 6: '1/30240'}

magnus.expansionterms.omega_terms(order: int) → Tuple[Term, Ellipsis][source]

Returns the terms of \(\Omega_n\) for n = order.

Each term is a (coefficient, word) pair: an exact fractions.Fraction and a nested commutator built from 'A' and ('Om', m) leaves via ('c', X, Y) nodes. The whole of \(\Omega_n\) is the integral of the sum of these terms, as in the module docstring.

Odd Bernoulli numbers above \(B_1\) vanish, so the \(j = 3, 5, 7, \ldots\) groups contribute nothing and are skipped rather than generated and discarded.

Added in version 1.0.0.

Parameters:: order (int) – Order \(n\) of the term; must be >= 1.
Returns:: The terms of \(\Omega_n\). For order=1 this is the single term \((1, A)\), the integrand of \(\Omega_1 = \int A\).
Return type:: tuple of (fractions.Fraction, tuple)
Raises:: ValueError – If order is less than 1.

Examples

import magnus.expansionterms as et

for coeff, word in et.omega_terms(3):
    print(f"{str(coeff):>6s}  {et.format_term((coeff, word), with_coeff=False)}")

  -1/2  [Om_2, A]
  1/12  [Om_1, [Om_1, A]]

magnus.expansionterms.magnus_terms(max_order: int) → Dict[int, Tuple[Term, Ellipsis]][source]

Returns omega_terms() for every order from 1 up to max_order.

Added in version 1.0.0.

Parameters:: max_order (int) – Highest order to generate; must be >= 1.
Returns:: Maps each order \(n\) to the terms of \(\Omega_n\).
Return type:: dict

Examples

import magnus.expansionterms as et

{n: len(terms) for n, terms in et.magnus_terms(8).items()}

{1: 1, 2: 1, 3: 2, 4: 3, 5: 5, 6: 9, 7: 17, 8: 33}

magnus.expansionterms.count_terms(order: int) → int[source]

Returns the number of terms in \(\Omega_n\).

Builds the terms and counts them, so the count reflects the collection of like terms rather than the raw size of the recursion.

Added in version 1.0.0.

Parameters:: order (int) – Order \(n\); must be >= 1.
Returns:: Number of distinct commutator terms in \(\Omega_n\).
Return type:: int

Examples

import magnus.expansionterms as et

{n: et.count_terms(n) for n in range(1, 11)}

{1: 1, 2: 1, 3: 2, 4: 3, 5: 5, 6: 9, 7: 17, 8: 33, 9: 65, 10: 129}

magnus.expansionterms.format_term(term: Term, with_coeff: bool = True) → str[source]

Renders one term as a readable string.

Added in version 1.0.0.

Parameters:

term ((fractions.Fraction, tuple)) – A term, as returned by omega_terms().
with_coeff (bool, optional) – If True (default), prefix the commutator with its coefficient.

Returns:

For example, '+1/12 [Om_1, [Om_1, A]]'.

Return type:

str

Examples

import magnus.expansionterms as et

[et.format_term(t) for t in et.omega_terms(3)]

['-1/2 [Om_2, A]', '+1/12 [Om_1, [Om_1, A]]']

magnus.expansionterms.print_magnus_terms(max_order: int, file=None) → None[source]

Prints the Magnus expansion, order by order, up to max_order.

Each order is printed as the integrand of \(\Omega_n\), one term per line, with its exact rational coefficient.

Added in version 1.0.0.

Parameters:

max_order (int) – Highest order to print; must be >= 1.
file (file-like, optional) – Destination, forwarded to print(). Default: standard output.

Return type:

None

Examples

import magnus.expansionterms as et

et.print_magnus_terms(4)

Omega_1  (1 term)
    int A

Omega_2  (1 term)
    int -1/2 [Om_1, A]

Omega_3  (2 terms)
    int -1/2 [Om_2, A]
    int +1/12 [Om_1, [Om_1, A]]

Omega_4  (3 terms)
    int -1/2 [Om_3, A]
    int +1/12 [Om_1, [Om_2, A]]
    int +1/12 [Om_2, [Om_1, A]]