magnus.expansionterms

expansionterms.py

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

The numerical core in magnus.magnus writes the coefficients of each commutator group out explicitly in Python through order 6, and builds orders 7 to 10 from the same recursion when it runs; the ceiling is magnus.magnus.MAGNUS_EXP_ORDER_MAX, set by how far the table of group factors reaches. This module derives those same terms from the recursion itself, in exact rational arithmetic, at any order – so the built-in 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, 9 at order 6, 33 at order 8, 129 at order 10, the counts count_terms() returns – 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

[1]

S. Blanes, F. Casas, J. A. Oteo & J. Ros, “The Magnus expansion and some of its applications”, Phys. Rep. 470, 151 (2009).

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, ...])	Returns the terms of \(\Omega_n\) for n = order.
magnus_terms(→ Dict[int, Tuple[Term, ...]])	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

magnus.expansionterms.Term

magnus.expansionterms.bernoulli(n: int) → fractions.Fraction

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

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

Raises:

ValueError – If j is negative. The check lives in bernoulli(), which this calls.

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, ...]

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\). A word is a nested commutator, written as a tuple, except at order=1, where the single term is (1, 'A') – the word is the bare string, the integrand of \(\Omega_1 = \int A\).

Return type:

tuple of (fractions.Fraction, Word)

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, ...]]

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

Raises:

ValueError – If max_order is less than 1.

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

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

Raises:

ValueError – If order is less than 1. The check lives in omega_terms(), which this calls to build the terms it counts.

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

Renders one term as a readable string.

Added in version 1.0.0.

Parameters:

• term ((fractions.Fraction, Word)) – A term, as returned by omega_terms(). The word is a tuple at every order above 1, and the string 'A' at order 1.

• with_coeff (bool, optional) – If True, prefix the commutator with its coefficient. Default: True.

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

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

Each order opens with a header naming it and counting its terms, then prints the integrand of \(\Omega_n\), one term per line, with its exact rational coefficient. Order 1 is the exception: its single term has coefficient 1 and is printed as int A, without one.

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]]