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

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 \eqalign{ & x_1 = B(x) = B \cr & x_2 = BB_1 \cr & x_3 = BB_1^2 + B^2 B_2 \cr}

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

The first few terms are:

$x(1) = B$,

$x(2) = B*B(1)$,

$x(3) = B*B(1)^2 + B^2*B(2)$,

$x(4) = B*B(1)^3 + 4*B^2*B(1)*B(2) + B^3*B(3)$

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),...) of nonnegative integers such that $\sum (i-1)*s(i)$ equals $n-1$ of

$$(-1)^{T(S)} * (T(S)+n-1)!A^{(-n-T(S))}* \prod_{i=2}^n (1/s(i)!)*((1/i!)*d^{(i-1)}A/dx^{(i-1)})^{s(i)}$$

where $T(S) = \sum_{i=2}^n i*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 = sum of (-1)^k * k! * B^(-1-k) * sum over all sequences V=(v(1),v(2),...) of the product of (B(i)/i!)^v(i) /v(i)! into the equation above and expanding and collecting all similar monomials in the Bs.

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 equal 0 to n-1 appearing with nonzero coefficient in x(n) satisfies u(0) = 1 + sum of (i-1)*u(i) from i=2 to n-1 and u(1) = n-1 - sum of i*u(i) from i=2 to n. I solved for u(0) and u(1) in terms of the "slack"(?) variables u(2),...,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

x(n) = B*B(1)^(n-1) + (2^(n-1) - n)B^2(B(1))^(n-3)*B(2) +

(1/4)(3^n - 3 - 2^(n+1)(n-1) + 2(n-1)^2)B^3(B(1))^(n-5)*(B(2))^2

• (1/4)*(3^(n-1) - 2^(n+1) + 2*n+1)B^3(B(1))^(n-4)*B(3) +...? +
• ((n^2-3*n+4)/2) * B^(n-2)*B(1)*B(n-2) + B^(n-1)*B(n-1)

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

-
-1 for (a) Not using LaTeX when the tool is given to you, and when your expressions look a mess without proper typesetting (b) Unnecessary shouting (c) Unnecessary editorializing in the last sentence (d) really bad abuse of notation (really? $B(k) := \partial_x^k B$ and $x(k) := \partial_t^kx$?) –  Willie Wong Oct 4 '10 at 17:31
Ah, multiple equal signs are used by the MarkDown software to indicate "section heading". I thought it was deliberate, I apologize. But do try to fix up your math display. Once you do that I'll remove the down-vote. –  Willie Wong Oct 4 '10 at 17:48
Whoa, it took a while just to understand what you mean. Let me rephrase the question, for the benefit of other readers: "If $dx(t)/dt=B(x(t))$, then what is $(d/dt)^k B(x(t))$, expressed in terms of $B(x)$, $B'(x)$, $B''(x)$, etc.?" –  Hans Lundmark Oct 4 '10 at 18:03
Hm, incidentally I made such calculations when prepared my thesis. I remember that the coefficients were somehow related to the Eulerian numbers. Unfortunately, I didn't find closed expressions. –  Nurdin Takenov Oct 4 '10 at 18:18
@Resolvent: I removed the downvote, since you apparently are not familiar with TeX. I still highly suggest that you learn it, but feel that it is unfair to penalize you for that. Also, perhaps you may re-edit the question to incorporate Hans' remark? His notation is somewhat more standard and easier to understand. –  Willie Wong Oct 4 '10 at 21:39

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.

-
Thank you! Your links led me to en.wikipedia.org/wiki/Butcher_group While it is not a final answer (they give a recursive formula by Cayley of the 1860s using Butcher's notation of the 1960s) perhaps it can lead me to it or will set me on the right path to figure it out for myself. –  resolvent Oct 5 '10 at 4:31
The formula you look for bears the name of Faà di Bruno. –  Denis Serre Oct 5 '10 at 5:26
This theory is described in Chapter 3 of the book "Geometric Numerical Integration" by Hairer, Lubich & Wanner. books.google.com/books?id=T1TaNRLmZv8C&printsec=frontcover –  Hans Lundmark Oct 5 '10 at 7:57
@Denis Serre: To be precise, it's not only Faà di Bruno, but rather what you get if you compute $x''(t)$, $x'''(t)$, etc., using Faà di Bruno, and recursively substitute the previously computed derivatives into the right-hand side at each stage, in order to get the answer in terms of $B$ and its derivatives alone (no $x$'s). –  Hans Lundmark Oct 5 '10 at 8:05
As Hans Lundmark correctly points out, it is more than just Faa di Bruno or just Lagrange Inversion Formula (I use the LIF as presented by G.P. Egorychev). In my original post, I mention that I use Faa di Bruno to express derivatives of A(x) in terms of those of B(x) where A(x)=1/B(x) and substitute into the LIF for the inversion of t=\integral_A(x)dx. I have also considered the generalizations of LIF which express f(x,t) for an arbitrary function, f, in terms of x when inverting a series G(x,t)=0, but all paths lead to massively complex calculations. –  resolvent Oct 5 '10 at 12:44

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.

-