Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

I would like to find the MacLaurin expansion of an iterated function. Finding the first few terms is not hard, but it doesn't take long before Mathematica runs out of memory using the straightforward program.

Is there some good way to find terms with reasonable time and space?

My problem allows some flexibility with the number of iterations and terms—but obviously more of both would be better. The goal is to be able to compute the terms at some high precision without loosing precision. This will only work when the input happens to be small, of course. The more iterations the more work is done per iteration but the smaller the range to which it applies; the more terms the larger the range.

At the moment I'm using 100 iterations and 150 terms (though this function is odd so it's only half that many coefficients). The eventual precision will be tens of thousands of digits, which causes some difficulties since the coefficients are large (on the order of $10^{133}$) even out to $x^{149}.$

share|improve this question
I'm guessing it will be quite complicated. Faa di Bruno for the case of $f(f(x))$ already looks quite complicated for me. –  Guess who it is. Sep 23 '11 at 16:39
@J.M.: Indeed, your intuition is spot-on. Mine is a simple formula, before iteration, but it takes over two gigabytes of memory and about an hour to get the terms mentioned above. –  Charles Sep 23 '11 at 16:42
In any event, you'd probably be interested in the matrix Bell polynomials. What version of Mathematica are you using? –  Guess who it is. Sep 23 '11 at 16:46
@J.M.: 8.0.1. What is this about matrix Bell polynomials? Haven't heard of those before. –  Charles Sep 23 '11 at 16:53
To use Mathematica notation: D[Nest[f, z, k], {z, n}] == BellY[Table[Map[Derivative[j][f], Reverse[NestList[f, z, k - 1]]], {j, n}]]. The only reason I know these is that I had been fiddling around with Faa di Bruno lately. I'm new to this myself, that's why I upvoted your question. :) –  Guess who it is. Sep 23 '11 at 17:03

1 Answer 1

up vote 4 down vote accepted

I don't really know about the general parameters of your problems at hand; but possibly my method using Pari/GP allows rather flexible work with up to 128 (or:why not 160) terms of a series and my standard precision is 200 dec digits. I've built a set of routines to express functional iteration as matrix-operations, and a flexible handling for that iterations (sometimes even fractional) by diagonalization, if applicable. If I encounter series with coefficients of 1e100, but signs are alternating then I can -in the same environment- use Euler-summation to arrive at results anyway.
In principle, what I'm doing is to construct a vector A1 containing the leading, say 128, coefficients of a power series of a function under consideration. Then I create the same way the vectors $\small A_2,A_3, \ldots $ having the coefficients of the according powers of the function, easily done by a standard call in Pari/Gp, but can also easily programmed. For the zeroth power it is just the vector $\small A_0= [ 1 ,0,0,...] $ up to the selected dimension. Then all the vectors are concatenated to a matrix $\small A = [A_0, A_1,A_2, \ldots A_n] $ (Such matrices are called Carleman-Matrix btw)
With a vector $V(x)$ which contains the consecutive powers of x (or its current value) I get then by the call $ \small V(x) \cdot A $ the result $\small V(f(x)) $ to a certain approximation.

Then function composition reduces down do matrix-multiplication $\small V(x) \cdot A \cdot B$ gives $\small V(g(f(x))$ if $\small B$ is the Carlemanmatrix of the function $\small g(x)$ , and iteration to matrix-powers $\small f(f(x)) = V(x) \cdot A^2 $ . In Pari/GP this can be done to a certain simple extent also with symbolic coefficients (but then not big matrix-size, say up to 64 x 64 at most.)

Unfortunately I don't have mathematica, so I can't be helpful with this, but I could post my couple of basic routines to show how this is done in principle.

There is a lot of finetuning, whereabouts, types of possible functions etc so I really do not know, how far this all is meaningful for your current problem at all.

[added] With a small set of basic matrices and the concept of matrix-decomposition this formalism allows a very handy upn-like environment for manipulating/analyzing powerseries. Having the basic algebraic operations at hand:
Operation of increment, or when repeated of addition, equivalents recentering of power series. It uses the Pascalmatrix P and its powers as Carlemanoperator

$ \qquad \small \begin{eqnarray} V(x) \cdot P^\tau &=& V(x+1) \\ V(x) \cdot (P^\tau)^{-1} &=& V(x-1) \\ V(x) \cdot (P^\tau)^y &=& V(x+y) \\ \end{eqnarray}$

Operation of multiplication, equvalents rescaling of powerseries, when a small prefix d indicates the diagonal-matrix use of a vector

$ \qquad \small \begin{eqnarray} V(x) \cdot dV(y) &=& V(x*y) \\ V(x) \cdot dV(y)^h &=& V(x*y^h) \\ \end{eqnarray}$

Operation of exponentiation, basically uses the matrix of Bell-polynomials, and logarithmizing using the factorially rescaled matrix of Stirling numbers second and first kind respectively
$ \qquad \small \begin{eqnarray} V(x) \cdot fS2F &=& V(\exp(x)-1) \\ V(x) \cdot fS1F &=& V(\log(1+x)) \\ \end{eqnarray}$

one can then - just like with an upn-calculator - define the operation for the "square-root" Carleman-matrix for $\small f(x) = \sqrt{x} $ by

$ \qquad \small \begin{eqnarray} V(x) &\cdot& {P^\tau} ^{-1} &=& V(x-1) \\ V(x-1) &\cdot& fS1F &=& V(\log(x)) \\ V(\log(x)) &\cdot& dV(1/2) &=& V(\log(x)/2) \\ V(\log(x)/2) &\cdot& fS2F &=& V(\exp(\log(x)/2)-1) \\ V(\exp(\log(x)/2)-1) &\cdot& P^\tau &=& V(\exp(\log(x)/2)) \\ \end{eqnarray} $

all together

$\qquad \small(V(x) \cdot {P^\tau}^{-1}) \cdot (fS1F \cdot dV(1/2)\cdot fS2F \cdot P^\tau )= V(\sqrt x ) $

somehow like $\small x \circ \text{DEC } \circ \text{LOGPLUS } \circ \text{MUL } 1/2 \circ \text{EXPMINUS } \circ \text{INC } \equiv x \circ \text{SQRT } $

and one shall see, that extracting the middle part into a separate matrix $\small SQ = fS1F \cdot dV(1/2) \cdot fS2F $ provides the carleman-matrix for $\small V(x)\cdot SQ = V(\sqrt{1+x}-1) $ with the expected coefficients for that function.

So to have a set of basic matrix-constants and basic procedures that is a toolbox for a very flexible handling of composition and iteration of functions which allow representation as power series.

Even more: if one needs fractional arguments, like with addition of fractional y in the formula above, one can introduce the general method of fractional powers of that matrices using the matrix-logarithm or diagonalization; especially the fractional power of P for fractionally-often increments is easy to define/implement.

All the above gives only an overview of the idea; there are some caveats, since we deal with matrix-products of matrices, which are ideally of infinite size, and for instance the partial product in the above formula for the sqrt of x, $\small {P\tau}^{-1} fS1F$ cannot be taken out and be made explicite (the associativity "breaks") due to occuring singularities. However, taking $\small V(x) \cdot {P^\tau}^{-1} =V(x-1) $ first allows to proceed - which means nothing else than that for the definition of power series for certain operations we need the recentering essentially .

[end additions]

For advanced handling with that problematic there is online-available material of R.P.Brent, who has in-depth analyzed composition of power series in its algorithmical implementation. I found for instance R.P.Brent "Fast Algorithms for Manipulating formal Power series" in "Journal of the Association for Computing Machinery, Vol 25, No 4, October 1978, pp 581-595" or "on the complexity of composition (...)" in "Siam J. Comput. Vol. 9, No.1 Feb 1980" and a handful related articles, but Brent's books should be present in many math-dept. libraries.

share|improve this answer
I doubted at first but it looks like PARI/GP vastly outperforms Mathematica in this calculation. For some reason Mathematica gives the wrong answer, and PARI the right answer, when intermediate precision is the same as final precision. Math'ca gives the right answer when it's given the result in symbolic mode but it takes forever. –  Charles Sep 23 '11 at 19:48
Yes, that may be, as long as we work numerically. With symbolic coefficients a 64x64-Carlemanmatrix in Pari/GP needs also quite a time. Moreover, Mathematica has the ability to sometimes find analytical closed-form expressions when only finitely many leading coefficients of a powerseries are given which is impossible in Pari/GP, thus allows exact results in a wider class of functions. –  Gottfried Helms Sep 23 '11 at 20:34
Concerning coefficients of huge size there is sometimes the possibility to handle (and then evaluate) using methods of Euler- or Noerlund summation. This can easily introduced as simple matrix-multiplications in the same environment as sketched above. A worked example dealing with a question in MO is here go.helms-net.de/math/tetdocs/_mo/mo_101217/… ; in that example the powerseries/Carlemanmatrix has coefficients of 1e100 and Noerlund-summation can define an equivalent which means to reduce them to 1e-120 or so... –  Gottfried Helms Sep 23 '11 at 20:38
Whew! Power[] was a bit tough to handle... anyway, for those who want to play around with Carleman matrices in Mathematica: CarlemanMatrix[f_, {x_, x0_}, {m_Integer, n_Integer}] := Prepend[Table[If[k == 0, Function[x, f][x0]^j, BellY[Table[{FactorialPower[j, i] (Which[#2 == 0, 1, #1 == 0, 0, True, #1^#2]) &[Function[x, f][x0], j - i], Derivative[i][Function[x, f]][x0]}, {i, k}]]/k!], {j, m}, {k, 0, n}], UnitVector[n + 1, 1]] –  Guess who it is. Sep 24 '11 at 10:43

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.