Formula for n-th iteration of dx/dt=B(x)

Question

Let $B(x)$ be infinitely differentiable with respect to $x$. Drop the use of parentheses on $B$ to delimit the argument $x$ and use them instead to hold the order of the derivative with respect to $x$. i.e. $B(0) = B$, $B(1) = dB/dx$, etc.

Let parentheses on $x$ hold the order of the derivative of $x$ with respect to $t$.

So \begin{align*} x_1 & {}= B(x) = B \\ x_2 & {}= BB_1 \\ x_3 & {}= BB_1^2 + B^2 B_2. \end{align*}

Is there a "nice" formula for the integer coefficient of an arbitrary monomial $B(0)^u(0) \cdot B(1)^u(1) \cdot\dotsb\cdot B(n-1)^u(n-1)$ in $x(n)$?

The first few terms are: \begin{align*} x(1) & {}= B, \\ x(2) & {}= B\cdot B(1), \\ x(3) & {}= B\cdot B(1)^2 + B^2\cdot B(2), \\ x(4) & {}= B\cdot B(1)^3 + 4\cdot B^2\cdot B(1)\cdot B(2) + B^3\cdot B(3). \end{align*}

One of my (many) approaches involved defining $A = 1/B$ so that $A(x)dx/dt = 1$. Then integrating and solving with the Lagrange Inversion Formula yields

$x(n)$ = sum over all sequences $S=(s(2),s(3),\dotsc)$ of nonnegative integers such that $\sum (i-1)\cdot s(i)$ equals $n-1$ of $$(-1)^{T(S)} \cdot (T(S)+n-1)!A^{(-n-T(S))}\cdot \prod_{i=2}^n \frac1{s(i)!}\left(\frac1{i!}\cdot\frac{d^{(i-1)}A}{dx^{(i-1)}}\right)^{s(i)} $$ where $T(S) = \sum_{i=2}^n i\cdot s(i)$.

I know I can simply make the upper limits in these products be infinity, because all but finitely many of the terms in the products are 1, because all but finitely many of the $s(i)$ in any sequence $S$ are zero.

But, for future computational purposes, I want to drag around the finite limits to remind myself when it comes time to implement on a computer.

So then I tried substituting the Faa da Bruno formula for $$ A(n) = d^A/dx^n = \text{sum of $(-1)^k\cdot k!\cdot B^{-1-k}\cdot\sum_{V=(v(1),v(2),\dotsc)}\prod \left(\frac{B(i)}{i!}\right)^{v(i)}\frac1{v(i)!}$} $$ into the equation above and expanding and collecting all similar monomials in the $B$s.

But, I cannot visualize a simple formula for the way all the terms combine.

---

So now I tried computing the terms of this sequence directly.

Homogeneity immediately tells you that any monomial product $B(i)^u(i)$ from $i = 0$ to $n-1$ appearing with nonzero coefficient in $x(n)$ satisfies $u(0) = 1 + \sum_{i = 2}^{n - 1} (i-1)\cdot u(i)$ and $u(1) = n-1 - \sum_{i = 2}^n i\cdot u(i)$. I solved for $u(0)$ and $u(1)$ in terms of the "slack" (?) variables $u(2), \dotsc, u(n-1)$ because the terms whose coefficients I CAN compute are most easily expressed in this form.

So, all I've got so far is \begin{align*} x(n) ={} & B\cdot B(1)^{n-1} + (2^{n-1} - n)\cdot B^2\cdot (B(1))^{n-3}\cdot B(2) +{} \\ & {}(1/4)(3^n - 3 - 2^{n+1}\cdot(n-1) + 2(n-1)^2)\cdot B^3\cdot(B(1))^{n-5}\cdot(B(2))^2 \\ & {}+ (1/4)\cdot(3^{n-1} - 2^{n+1} + 2\cdot n+1)\cdot B^3\cdot (B(1))^{n-4}\cdot B(3) \\ & {}+ \cdots? +{} \\ & {}+ \left(\frac{n^2-3\cdot n+4}2\right)\cdot B^{n-2}\cdot B(1)\cdot B(n-2) + B^{n-1}\cdot B(n-1). \end{align*} I could keep going, deriving longer and longer formulae for more of the terms, and then HOPE that I can guess the general pattern for all of them.

$B(i)^{u(i)}$ means $B(i)$ raised to the $u(i)$-th power.

Answer 1

These expansions can be described in terms of rooted trees. The first few coefficients are easy to derive by hand, but rooted trees provide you a way of generating the coefficients to arbitrary order, see here. You start with the trivial tree, and at each stage of the derivation you add another level to the tree according to specific rules. The coefficients that occur in the tree have various number-theoretic properties.

This is a very interesting part of combinatorics, with applications in numerical analysis and quantum field theory (the only two fields that I understand). As for numerical analysis, the design of RK methods of arbitrary order didn't go anywhere precisely because of the calculatory issues that you experienced. It was actually a pretty big achievement by John Butcher in the 1960s that he was able to describe these prolongations of ODEs in a compact way, and hence provide a description of RK methods of arbitrary order.

Answer 2

Caveat emptor (July 5, 2021, expanded on the 12th), the devil is in the details:

As I point out in the addendum below, the OP's polynomials provide solutions to a (the) classic flow equation, and, although this ODE has Lie written all over it and a history extending beyond him, I have not found these polynomials in any publication prior to 2008 (Edit, Feb 12, 2022: found, see this MO-Q). The OPs polynomials, which I call the refined Eulerian partition polynomials (REPs), presented for the first time as far as I know in my 2008 OEIS entry A145271, are introduced in the addendum below. The REPs are the expansions of $(g(x)D)^n g(x)$ in terms of the $D^k g(x) = \frac{d^k g(x)}{dx^k}$. If you should find them in any publication other than those I've already compiled there, please point it out to me and I will include it in the entry (or you can--it's up to you). (I have included refs in the entry to publications prior to 2008 that present perspectives from those on flow equations, renormalization, and functional iteration dating from Charles Graves in 1850-53, Abel, and Schroeder, but they do not contain the REPs.)

I have found in the literature predating 2008 two other closely related manifestations of iterated derivatives--the Butcher elementary differentials (B-diffs), enumerated by the Connes-Moscovici weights (CMW) of A139002, and the Comtet diff ops, characterized by the A-polynomials of A139605, for the normal-ordering expansion of $(g(x)D)^n$, i.e., in terms of summands with all $D^k$ to the right. 'Naturally grown' forests of nonplanar rooted trees (and other types of trees) can be associated with these three sets of iterated/nested derivatives. 'Naturally grown" is explained in the addendum below. Pages 8 and 10 of "Lie–Butcher series and geometric numerical integration on manifolds" by Lundervold give examples of the B-diffs and associated trees for a vector function.

Hairer, Lubich, and Wanner, in "Geometric Numerical Integration - Structure-Preserving Algorithms for Ordinary Differential Equations" (2nd Ed., Springer, 2006) on p. 51, consider the flow equation

$$\dot{y} = f(y)$$

with initial condition $y(t_0)=y_0$.

The derivatives of $y$ in terms of the derivatives of $f$ are presented in Equation 1.4 on page 52 (decidedly NOT the Bell partition polynomials). Table 1.1 on p. 53 lists the first few corresponding forests of trees; $\alpha(\tau)$, enumerating the non-planar rooted tree types $\tau$; and the corresponding elementary differentials $F(\tau)$, in bijection with the trees. The $n$'th derivative of the solution is given as

$$y^{(n)}(t_0) = \sum_{|\tau|=n} \; \alpha(\tau) \; F(\tau)(y_0),$$

where the cardinality $|\tau|$ is the number of vertices/nodes/knots of the nonplanar rooted tree $\tau$. This is a key identity of the theory of Butcher series.

The tree types $\tau$ are enumerated by the CMW of A139002, i.e., $\alpha(\tau)$, and the trees are in bijection with the B-diffs as denoted in Hairer et al.; however, this bijection dissolves when the tree types are related to the ops $(g(x)D)^n$ of the Comtet iterated Lie ops, with many weights merging for the larger forests. That is, more than one tree type can be associated to the components of the Comtet ops.

This contraction continues when the Comtet ops act on $g(x) = f(x)$, giving the polynomials of A145271, the set of polynomials the OP is trying to characterize. So, there is a refinement proceeding from the array of coefficients of the REPS to those of the Comtet ops to the array of CMW. For example, in the notation of Hairer et al., with the CMW $(1,\;3,\;1,\;1)$,

$$y^{(4)} = f'''(f,f,f) + 3f''(f'f,f) + f'f''(f,f) + f'f'f'f$$

$$ = f^3f''' + 3 f^2f'f'' + f^2f'f''+ (f')^3 f $$

$$ = (fD)^3 \; f = [\;(f(f')^2 + f^2 f'')\; D + 3f^2f' \; D^2 + f^3 \; D^3\;] \;f$$

$$ = f^3f''' + 4 f^2f'f'' + (f')^3 f.$$

The last partition polynomial is a (not so) refined Euler partition polynomial (A145271), the diff op in the line just before it is a Comtet op (A139605), and the first line contains the B-diffs with the Connes-Moscovici weights (A139002), all corresponding to the naturally grown forest of six trees with four vertices.

Page 165 of the 3rd edition of "Numerical Methods for Ordinary Differential Equations" by Butcher contains a table similar to the one in Hairer et al. Theorem 311(B) on p. 167 is equivalent to the equation for $y^{(n)}$ above.

Original posting (2011):

You might also try following the references and links in A145271 ["Coefficients for expansion of $[g(x)d/dx]^n g(x)$; refined Eulerian numbers for calculating compositional inverse of $h(x)= (d/dx)^{-1} 1/g(x)$"] and A139605 ["Weights for expansion of $(f(x)D_x)^n$: coefficients of A-polynomials of Comtet"] of the Online Encyclopedia of Integer Sequences for some ideas. Note that the coefficients of A145271 (your expansion coefficients) are embedded in A139605.

An explicit formula for the coefficients are given in the Comtet reference in A139605 on page 166 as egn. (8) with l=1.

Addendum (July 4, 2021):

The OP is essentially asking for the solution to the classic Lie autonomous ODE, the flow equation discussed in A145271,

$$y' = f'(x) = g(f(x)) = g(y),$$

in terms of derivatives of $g$. The solution is given by the exponential of the iterated vector, or Lie infinitesimal generator (infinigen) $g(x) D_x$ as

$$y = f(x) = e^{(x \; g(w) \; D_w)} \; w \; |_{w=0} = f[ \; x + f^{(-1)}(w)\;] = |_{w=0} $$

with

$$g(w) = \frac{1}{D_w \; f^{(-1)}(w)},$$

where $f^{(-1)}(w)$ is the compositional inverse of $f(z)$ about the point $(z,w)$ and $f(0) =f^{(-1)}(0) =0$.

This is easily derived noting that, with $(z,w) = (f^{(-1)}(w),f(z)),$

$$ g(w) \; D_w = \frac{d}{(f^{(-1)})'(w)dw}= \frac{d}{df^{(-1)}(w)} = \frac{d}{dz} = D_z ,$$

so, acting on a function $H$,

$$e^{(x \; g(w) \; D_w)} \; H(w) = e^{(x \; D_z)} \; H(f(z)) = H[f(x+z)] = H[\;f(x + f^{(-1)}(w))\;].$$

This is a chimera of the inverse function and the Taylor series expansion theorems. (One of the Graves brothers, collaborators of Cayley, published this result in the 1850s.).

Specifically, the OP is seeking methods of constructing the polynomials in $D_x^k g(x)$ of $(g(x)D_x)^{n-1} \; g(x)$. Examples of these refined Eulerian polynomials are given in A145271 along with a matrix formula for computing them and references (see also refs 1 and 2 below). One reference is to the 1857 paper of Cayley on operandators and their depiction as trees via which these polynomials may be generated or represented (see the 'normally grown' trees depicted in the "Mathemagical Forests" link in the OEIS entry, in my blog post "Lagrange a la Lah", and in 4 below). Cayley considers multivariable functions and operators.

The relation to rooted trees follows from the product formula

$$D \; f(x)h(x) = f'(x)h(x) + f(x)h'(x)$$

and pre-Lie action (PLA) of the infinigen

$$g(x) D \; f(x) \; D \to g(x)f'(x) \; D.$$

(This is called a connection or covariant derivation of vector fields in "Butcher series: A story of rooted trees and numerical methods for evolution equations" by McLachlan, Modin, and Verdier.)

This PLA of the infinigen allows the iterative action $(g(x)D)^{n-1} \;g(x)$ to be illustrated diagrammatically as PLA action down the branches of forests of rooted trees with initially a copy of the infinigen attached to each node of each tree except for the root, which has only $g(x)$ attached to it. Each iteration corresponds to growing a grove of $n$ trees with $n+1$ nodes, or vertices, from each tree $T_n$ with $n$ nodes by attaching a branch to each node of $T_n$ with each attachment resulting in a new tree. The resultant total action follows from the flow of PLA down along the branches from the leaves to the root for all the trees of this new 'naturally grown' forest.

For example, the 'seed' or root, $T_1$, with only one node, represents the function $g(x)$. Grow the tree $T_2$ from $T_1$ by attaching a branch to its single node, giving a tree with two nodes-- a tree with a single trunk--with $g(x)$ at the root and the infinigen $g(x)D$ at the top node, or leaf. $T_2$ then represents $g(x)D\;g(x) = g(x)\;g'(x)$. Grow a forest from this single tree by attaching a branch to each node of $T_2$, one at a time, giving two rooted trees with three nodes each, i.e., two $T_3$ trees. One is a tree with two trunks with $g(x)$ attached to the root and an infinigen to each of the other two nodes (leaves). Then the simultaneous independent PLA directed down the trunks gives $g^2(x) \; g''(x)$. The second tree with three nodes is a taller one with a single trunk and a single leaf with again $g(x)$ at the root and an infinigen at the other two nodes. This tree represents the action $g(x)\;(g'(x))^2 $--the PLA of the top infinigen on the next lower one gives $g(x) \;g'(x)D$ and then its action on the root $g(x)$ gives $g(x)\;(g'(x))^2$. Then the forest of two $T_3$ trees represents

$(g(x)D)^2 \; g(x) = g(x)D\; g(x) \; g'(x) = (g(x))^2\; g''(x) + g(x)\; (g'(x))^2$.

The next iteration of this 'natural growth' generates a new forest of six trees each with four nodes, and, in general, the number of trees in a forest of type $T_n$ trees is $(n-1)!$.

To grok more complicated trees, consider the example of the simultaneous PLA of three nodes acting via three branches on the node directly below with $F_1(x) \; D$ attached to one of the upper nodes; $F_2(x) \;D$, to another; and $F_3(x)\; D$, to the last. The three functions $F_n$ are generated by PLA directed along branch paths to these nodes from the leaves. With only these three nodes attached directly to the lower node, the PLA action on that node gives

$F_1(x)\;F_2(x)\;F_3(x)\; g'''(x)\; D$.

Continue growing more forests ad infinitum with an infinigen at each node other than the root and direct PLA down along the branches in the manner described above to generate $(g(x)D)^n g(x)$ for each $n$.

Note that the 'raising/creation' operation generating the next younger larger forest from a older smaller forest is the 'natural growth' operation. The 'lowering/annihilation/destruction' operation generating the the older smaller forest from younger larger forest is to burn down all the trees with more than one trunk, i.e., more than one branch attached to the root, and then remove the root from the remaining planted trees, i.e., trees with only a single trunk.

For additional history and refs, see the relevant contribution to the MO-Q "In 'splendid isolation'".

Related MO-Qs:

1. "Expansions of iterated, or nested, derivatives, or vectors--conjectured matrix computation"

2. "Closed formula for $(g \partial)^n$"

3. "differential operator power coefficients"

4. "Who invented diagrammatic algebra?" (Abdelmalek Abdesselam's answer and my comments to it)