Conference rusure::math

1570.0. "Repeated composition." by CADSYS::COOPER (Topher Cooper) Fri Feb 21 1992 17:26

From: dwilson@cvbnet.prime.com (David Wilson x5694 4-1600)
Newsgroups: sci.math
Subject: f o f = exp
Message-ID: <2768@cvbnetPrime.COM>
Date: 19 Feb 92 19:37:10 GMT
Organization: Computervision, A Prime Computer Company, Bedford, MA, USA
Lines: 226



    If you are interested in iterated functions from an experimental
    standpoint, here is a quick primer.  I do not keep junk around to
    repost, so save it while its hot.

        ---------------------------------------------------------

    Let f be a function from R to R.  Let [f^n] denote f composed on
    itself n times.  For n in Z+, [f^n] can be recursively defined as
    follows:

        [f^1] = f
        f o [f^a] = [f^(a+1)]

    From P1 and P2, we can conclude the following for n in Z+:

        [f^a] o [f^b] = [f^(a+b)]
        [[f^a]^b] = [f^(ab)]

    The natural extension of [f^n] to n in Z is:

        [f^0] = I               (I = identity function)
        [f^(-n)] = [f^n]^-1     (f^-1 = inverse of f)

    Extending to n in Q:

        [[f^(a/b)]^b] = [f^a]

    The natural extension to n in R would be:

        [f^n](x) = lim [f^m](x)
                   m->n

    All this, of course, is a gross oversimplification that ignores
    extremely messy existence, uniqueness, and domain issues.

        ---------------------------------------------------------

    For some functions [f^n] for n in R seems to have a "natural"
    definition.  This definition can be had by finding, by inspection,
    a formula for [f^n] with n in Z+ and extending the formula
    "naturally" to n in R.  For instance:

        f(x)                    [f^n](x)

        x + a                   x + na
        ax + b (a != 1)         (a^n)x + ((a^n-1)/(a-1))b
        bx^a (a != 1)           b^((a^n-1)/(a-1))x^(a^n)

    We will call such a formula the "natural formula" for [f^n].

        ---------------------------------------------------------

    Suppose now, that we have an f for which we cannot find a natural
    formula for [f^n].  Nevertheless, suppose we can find g such that

         lim  [g^-k][f^k]
        k->inf

    exists (on some subdomain of R).  It then turns out that

        [f^n] = lim  [f^-k][g^n][f^k].
               k->inf

    modulo domain difficulties.

        ---------------------------------------------------------

    For example, suppose we wish to find [f^n], where f(x) = x^2+1.
    Iterating f does not seem to yield a natural formula:

        [f^1](x) = x^2 + 1
        [f^2](x) = x^4 + 2x^2 + 2
        [f^3](x) = x^8 + 4x^6 + 8x^4 + 8x^2 + 5
        ...

    However, g(x) = x^2 has the natural formula [g^n](x) = x^(2^n).
    Also,

         lim  [g^-k][f^k]
        k->inf

    exists on the domain R (exercise for the reader).  We therefore
    have

        [f^n] = lim  [f^-k][g^n][f^k]
               k->inf

    Still letting f(x) = x^2+1, we wish to find [f^1/2](0) to about
    20 digits ([f^1/2] composed on itself yields f):

	1.  Choose integer k large enough that (f o [f^k])(0) =
	    (g o [g^k])(0) to within relative accuracy 10e-20.
	    In this case, k = 8.

	2.  Compute [f^k](0), i.e., iterate f 8 times on 0, giving
	    about 4.4127887745906175987e22.

	3.  Perform [g^1/2] on the result.  The natural definition
	    of [g^1/2](x) is x^sqrt(2).  This gives about
	    1.0579385394282632787e32.

	4.  Perform [f^-k] on the result, i.e., perform [f^-1] k
	    times on the result, where [f^-1](x) = sqrt(x-1).  This
	    gives 6.4209450439082838148e-1.

    Hence, [f^1/2](0) = 6.4209450439082838148e-1 to about 20 digits.
    
    It turns out that [f^n](x), where f(x) = x^2+1, is very fast to
    compute for reasonable n and x.  The following is a K&R C routine
    to do the job.  (It will, however, give sqrt domain errors when
    [f^n](x) does not exist.)

        ---------------------------------------------------------

#include <math.h>

/* g(x) = x^2 */
double g(x) double x; { return x*x; }

/* g(x) = x^2 iterated n times by the natural formula [g^n] = x^(2^n) */
double g_iterated(n, x) double n, x;
{ return exp(log(x)*exp(log(2.0)*n)); }

/* f(x) = x^2+1 */
double f(x) double x; { return x*x+1; }

/* Inverse of f(x) = x^2+1 */
double f_inverse(x) double x; { return sqrt(x-1); }

/* f(x) = x^2+1 iterated n times by limit formula, to machine accuracy */
double f_iterated(n, x) double n, x;
{
    int k;

    /* Iterate f on x, counting k iterations, */
    /* until result is so large that f(x) and g(x) are equal */
    /* to within machine accuracy */
    for (k = 0; f(x) > g(x); k++)
	x = f(x);

    /* Apply [g^n] to result */
    x = g_iterated(n, x);

    /* Iterate inverse of f k times on f */
    for (; k > 0; k--)
	x = f_inverse(x);

    /* And we are done */
    return x;
}

        ---------------------------------------------------------

    Finally, to attack the iteration problem for exp, the exponential
    function.  There seems to be no g with a natural formula for [g^n]
    such that

         lim  [g^-k][exp^k]
        k->inf

    is defined for any subdomain of R.  We are therefore forced to use
    another approach to compute [exp^n] (in particular, h = [exp^1/2],
    which satisfies h o h = exp).

    The key is the Taylor series for exp

                  inf
	exp(x) =  sum  x^i/i!
		 i = O

    The partial sums of this series,

                   p 
	e_p(x) =  sum  x^i/i!
		 i = O

    form a sequence of polynomials which approximate exp as closely
    as desired on any desired range.  To compute [exp^n](x), for
    0 <= n <= 1, do the following,
    
	1.  Choose p so that f = e_p approximates exp to within the
	    desired accuracy near x.

	2.  Let g(x) = x^p/p!, the term of highest exponent of e_p(x).
	    It turns out that

	    1.	 lim  [g^-k][f^k],
		k->inf

	    exists on an infinite subdomain of R.  (This is true whenever
	    g(x), the term of highest exponent of polynomial f(x), has a
	    positive coefficient.  Exercise 2)

	    2.  [g^n] has a natural formula.

	3.  It therefore follows that we can compute [f^n](x) as we did
	    for f(x) = x^2+1, and for 0 <= n <= 1, we can reasonably
	    expect that [f^n] will approximate [exp^n] as closely as f
	    approximates exp on the same range (my unproved assertion).

	    The problem with this computation is that, for tiny
	    accuracies, k in 

	    [f^n] = lim  [f^-k][g^n][f^k]
		   k->inf

	    may have to be chosen very large in order to get [g^k](x)
	    to dominate [f^k](x) to within accuracy.  Also, [f^-1] will
	    have to be accomplished by an approximation method, such as
	    Newton's method.

	4.  It turns out, that for any real n >= 0, [exp^n] can be
	    extended to domain R.  I leave out the details as I am
	    unaware of them.

        ---------------------------------------------------------

    Have a nice day.

-- 
David W. Wilson (dwilson@cvbnet.prime.com)
J.H.Whitney (was Prime Computer (was Computervision Corp.)), Bedford, MA
Disclaimer: "Truth is just truth...You can't have opinions about truth."
- Peter Schikele, introduction to P.D.Q. Bach's oratorio "The Seasonings."