The Wolfram Emerging Leaders Program: High School (WELP:HS) is a remote, semester-long program designed for exceptional graduates from the Wolfram High School Summer Research Program to continue to develop their skills. Students work in a small group over the course of a semester to do either a research-based project or a product development–based project under the guidance of expert mentors. Projects are centered around a topic at the intersection of the group’s interests.

Authors

◼
  • Bertie Bennett
    • ◼
  • Christopher Gilbert
    • ◼
  • Gregory Roudenko

    • Abstract

    This project explores non-integer iteration of functions of one or more variables. We present algebraic solutions to roots of some generalized functions, and detail various numerical and algebraic root-finding algorithms.

    Introduction

    Basic Definitions

    Composition, denoted with the symbol ◦, involves inputting one function’s output into another’s input. For example,
    (f◦g)(x)
    is equivalent to
    f(g(x))
    . When the same function is composed with itself multiple times, (e.g. f(f(f(x)))) it is said to be nested.
    Let
    f◦f◦f...
    n
    where
    n∈
    (i.e. nesting) be denoted
    n
    f
    , where ◦ denotes composition.
    n
    f
    may be called the nth functional exponent of f.

    Implementation
    

    Properties

    Composition of the same nested function is additive, so
    n
    f
    ◦
    m
    f
    =
    n+m
    f
    .
    Nesting a nested function is multiplicative, so
    m
    (
    n
    f
    )
    =
    nm
    f
    .
    In these ways, functional exponentiation and composition have similar algebraic properties to exponentiation, such as the rules that
    a
    e
    b
    e
    =
    a+b
    e
    and
    b
    (
    a
    e
    )
    =
    a×b
    e
    .

    Generalized Nesting

    Integers

    Although nesting so far has only been defined for the natural numbers, it can be extended based on the properties above. If we define
    0
    f
    (x)=x
    and
    -n
    f
    =
    n
    (
    -1
    f
    )
    (where
    -1
    f
    is the inverse function of f), then the above properties of composition and nesting now hold for integer nestings.
    E.g. for some function
    f
    with an inverse
    -1
    f
    :
    ◼
  • By the definition of the inverse function,
    (
    1
    f
    ◦
    -1
    f
    )(x)=x
    • ◼
  • Equivalently, with our definitions and properties,
    (
    1
    f
    ◦
    -1
    f
    )(x)=
    1-1
    f
    (x)=
    0
    f
    (x)=x

    • Inverses and Rationals

    In the same way that composition being additive leads to a generalization to the integers, nesting being multiplicative leads to a generalization to the rationals.
    Let the m-th functional root be denoted
    1/m
    f
    . This function has the property that nesting it m times will result in the original function
    f
    , because
    m
    (
    1/m
    f
    )
    =
    m/m
    f
    =
    1
    f
    E.g.
    1/2
    f
    ◦
    1/2
    f
    =
    1/2+1/2
    f
    =
    1
    f
    In other words, f’s square root g must satisfy g(g(x)) = f(x).
    By the same logic, let the functional exponent
    n/m
    f
    be defined as the m-th root of
    n
    f
    , or, equivalently, the nth nesting of the m-th root of
    f
    .

    Reals

    For some
    w∈
    , define
    w
    f
    (x)=
    
    (i->∞)
    
    i
    2
    w
    
    i
    2
    
    f
    (x)
    . This is simply the limit of a rational approximation of w.

    Complexes

    Occasionally, it may be even possible to generalize functional iteration to complex. That is, the i-th exponent of f,
    
    f
    , must satisfy the property that
    
    (
    
    f
    )
    =
    2
    
    f
    =
    -1
    f
    .

    Goal of the Project

    In general, the goal of the project will be finding and exploring solutions or approximations of non-Integer functional exponents.
    Basic Examples of Explicit Solutions
    Beyond trivial examples where f(f(x)) = f(x) (e.g. f(x) = a, f(x) = x, f(x) = Ramp(x)), here we will explore algebraic solutions to simple functional exponents.

    Helper Functions
    

    Affine Functions

    An affine function is a function which can be written
    f(x)=a+bx
    . Using Solve[], one can find a simple, affine solution to the functional m-th root of f.
    This method only works if the structure of
    1/2
    f
    is known. That is, the affine nature of
    1/2
    f
    can be guessed, and because an exact solution exists, it can be seen that this guess is true. Such a structure becomes more difficult to find as a function becomes more complex, and is the main challenge in finding exact solutions to most functions.
    Set up affine function form (and use 1/7th as the iterate):
    In[]:=
    m=7;​​f[x_]:=a+bx​​fRoot[x_]:=r+sx
    Find
    m
    (
    1/m
    f
    )
    (x)
    . Then, solve for the coefficients in
    1/m
    f
    :
    In[]:=
    fRootNested=(fRoot@*m)[x]//PowerExpand//Expand
    Out[]=
    r+rs+r
    2
    s
    +r
    3
    s
    +r
    4
    s
    +r
    5
    s
    +r
    6
    s
    +
    7
    s
    x
    In[]:=
    fCoefs=coefList[f[x],{0,1}];​​fRootCoefs=coefList[fRoot[x],{0,1}];​​fRootNestedCoefs=coefList[fRootNested,{0,1}];
    Solve for the coefficients of
    (
    1/m
    f
    )
    :
    In[]:=
    (fullSolution=Solve[​​And@@Thread[#1==#2&[Values[fCoefs],Values[fRootNestedCoefs]]],​​Values[fRootCoefs[[;;-2]]]​​])//Column
    Out[]=
    r
    -a+a
    1/7
    b
    -1+b
    ,s
    1/7
    b
    
    r
    -a-
    1/7
    (-1)
    a
    1/7
    b
    -1+b
    ,s-
    1/7
    (-1)
    1/7
    b
    
    r
    -a+
    2/7
    (-1)
    a
    1/7
    b
    -1+b
    ,s
    2/7
    (-1)
    1/7
    b
    
    r
    -a-
    3/7
    (-1)
    a
    1/7
    b
    -1+b
    ,s-
    3/7
    (-1)
    1/7
    b
    
    r
    -a+
    4/7
    (-1)
    a
    1/7
    b
    -1+b
    ,s
    4/7
    (-1)
    1/7
    b
    
    r
    -a-
    5/7
    (-1)
    a
    1/7
    b
    -1+b
    ,s-
    5/7
    (-1)
    1/7
    b
    
    r
    -a+
    6/7
    (-1)
    a
    1/7
    b
    -1+b
    ,s
    6/7
    (-1)
    1/7
    b
    
    Return the only real solution:
    In[]:=
    sol=Normal[Association[Join[{r->0,s->0},First[fullSolution]]]]
    Out[]=
    r
    -a+a
    1/7
    b
    -1+b
    ,s
    1/7
    b
    
    Based on the above computation, a general m-th root formula is given by:
    1/m
    f
    (x)=a
    1/m
    b
    -1
    b-1
    +
    1/m
    b
    x
    ​This can be extended to a general exponentiation formula for rationals, reals, and complexes:
    z
    f
    (x)=a
    z
    b
    -1
    b-1
    +
    z
    b
    x
    Plot any affine function with a z-th iterate:
    Out[]=
    ​
    z
    a
    b
    -10
    -5
    5
    10
    -30
    -20
    -10
    10
    20
    30
    f(x)=a+bx
    z
    f
    ​(x)
    Plot of affine f and
    z
    f
    Plot any affine function for any real z:
    Out[]=
    ​
    a
    b
    Plot of
    z
    f
    (x) where f(x)=a+bx

    Power Functions / Monomials

    Generalized monomial:
    Plot of x^2 and its iterates:
    Plot of complex iterates of x^2:
    Inverse Function Composition
    Although easiest to identify in the case of affine and monomial, this can be generalized to all functions, and to all function-inverse function pairs:
    And so we arrive at the above formula. As this method of proof can be generalized to all rational exponents, it can be extended to the reals and complexes through a similar method as given above for the extending of functional nesting into functional exponentiation.

    Formulas found by this technique

    Ln(Sin(Exp(x)))

    And indeed, if this root is iterated with itself, it gives an approximation for f:
    Numerical Methods
    As many functional roots are very difficult to find, and in many cases (such as polynomials on the full complex plane), are proven to have no solutions, this project now examines methods for finding functional iterates (but mainly functional square roots) using numerical methods.

    Method 1: Recursive Function Values

    Input

    Input a function and a finite range over which the solution will be found. The step size indicates the distance between the points where error is calculated.
    Here is a plot of the inputted function over the range:

    Explanation

    We are trying to solve the equation g(g(x)) = f(x) for all x.
    ◼
  • Given an initial x, xi, we can make a guess as to the value of g:
    ​
    • ◼
  • g(xi) = w
    • ◼
  • Because we know that g(g(x)) = f(x), then g(g(xi)) = f(xi), and since g(xi) = w, then:
    • ◼
  • g(w) = f(xi)
    • ◼
  • This gives a second value of g. For a third, we can repeat the same thing. g(g(x)) = f(x), so g(g(w)) = f(w), and since g(w) = f(xi):
    • ◼
  • g(f(xi)) = f(w)
    • If you continue this, you get as many points as you want that will all lie on the g, assuming that g(xi) = w. By searching over all xi and all w, we can find the solution which minimizes the squared difference between a linear spline of the approximating points and the original function.
    This is a function which more efficiently implements the process described above:
    FInally, here are more helper functions relating to it, which are fairly self-explanatory.

    Search

    For every xii and w, the Search[] function checks the error of it, setting minP to be the values of xii and w which minimize error, and returning that minimum error. By then subdividing the region around minP, a better and better solution can be found. As the first iteration is the biggest search (in order to minimize the effect of local minima), it takes much longer than all other iterations.

    Output

    Based on the search, here is the best function found so far:
    And here are the xii and w values, and the error which were achieved:

    Analysis

    Although this method does find some points on the function, and the general “shape” of it for positive numbers (especially small positive numbers), it has a final total error (i.e. the area under the curve of the squared difference between the function and it’s approximation) of:

    Method 2: Simultaneous function value optimization

    Explanation

    Instead of being smart about nesting functions with a single initial guess, this method tries to determine the whole function at once in order to move the values closer to what is needed. It is extremely susceptible to local minima, but it can often do a much better job of optimizing simpler functions, such as sine or some polynomials.
    First, take a guess at a functional root g. We used the g from method 1 here. One possible route for further exploration is to randomize the initial function g, then search from there. In this way, some local minima for the optimization function may be avoided.
    Next, find the differences between that and the original function:
    delta(x) = f(x) - g(g(x))
    Finally, simply add the delta to update g:
    new_g(x) = g(x) + temperature * delta(x)
    For some points, g(x) will not be affected much, but g(g(x)) will be moved closer to f. In this case, g will be positively affected. For other points, it may shoot off to infinity or become discontinuous, so there is no real guarantee that this will work.

    Implementation

    First, this code represents both the original function and the approximation as a set of points, with xVals being the x values of the points, fVals being the corresponding y values of the original function, and aVals being the approximated y values. Using a linear spline, it approximates the corresponding functions from these lists of points.
    Here are plots of the interpolated values for the original function and the function from method 1:

    Optimization

    This repeatedly applies the method described to optimize the function.
    Plot the progress made so far:
    You can see that in the case of the example function, this method does eventually diverge. However, before then, it does find a quite close solution.

    Output

    Here is the output of the optimization:
    And the nested approximation compared to the original:

    Analysis

    This approximation already looks much more promising, and indeed, when the loss is found, it comes out to be much smaller, though with some improvement still possible.

    Method 3: Genetic Algorithm of Fourier and Taylor Series Approximations

    Integral Loss Function and Minimization

    Another method of finding a numerical solution to functional equations is to evolve functions from an initial guess using a genetic algorithm. The genetical algorithm minimizes an integral loss function defined here:
    We can either directly vary the values of a function at each point (by decomposing it into a list of x values and y values and mutating those) or instead represent our iterate approximant as a function belonging to a family (e.g. Taylor Series/Fourier Series) and mutate the coefficients.

    Built-In Optimization of Fourier Series

    For some functions like Sin and Cos, fourier series approximations of their roots converge very quickly, and so a good approximation can be found with relatively little computational power. For example, a 9th degree summation of sines can be minimized to a squared error of less than 10^-7:
    However, this method is quite slow and scales up exponentially with more terms. Using a different Method parameter for NMinimize may reduce this.

    Rust-Based Approximations

    Because of the speed concerns, and also for greater flexibility in the types of series-based approximations that we wanted to test, we wrote a genetic algorithm in the low-level programming language Rust, which can approximate functions by means of Fourier series, and Taylor series, and Linear splines.

    Results for Fourier Approximation

    For an example of the output of this code, here is one output for a fourier approximation of the square root of sin:
    Converting this into Wolfram Language, we get the following approximation:
    This approximation achieves an incredibly small loss, the best of any method tested so far:
    And it is indeed an almost perfect match for sin:

    Results for Taylor Series

    While not as perfect as for sin, this process still yields a small error value, especially when considering that this function was often not handled well by other approximation algorithms:

    Generalizability

    Although one can fine-tune the optimizations for better results, or choose between Taylor series,Ffourier series, and linear spline approximations, here are the loss values for a variety of functions all trained with the same parameters over the range [-1, 1]:
    ◼
  • f(x) = x^2 + 1 | 0.000111133
    • ◼
  • f(x) = 2x + 3 | 2.94388*10^-6
    • ◼
  • f(x) = cos(x) | 0.0189814
    • ◼
  • f(x) = e^x | 0.000010337
    • ◼
  • f(x) = ln(2 + x) | 0.0000279523
    • ◼
  • f(x) = {x < 0: x + 1, x^2} | 0.0768464
    • Although some functions achieve a much better result than others, this shows that there does generally exist a functional square root for many functions, and the existence of solutions (especially one over a fixed range), is not limited to a select few functions.
    For comparison with the previous two examples, we should ideally test with the same range of [-2 Pi, 2 Pi]. However, every method of approximation I use with the genetic algorithm seems to diverge, getting quite poor results (worse than the previous method by a couple orders of magnitude). Because of this, I’ve included both values:
    ◼
  • f(x) = Sin(x) + Sigmoid(x) where -1 < x < 1 | 8.05862*10^-7
    • ◼
  • f(x) = Sin(x) + SIgmoid(x) where -2 Pi < x < 2 Pi | 0.433118
    • Conclusion

    Results

    Throughout this paper, we discussed a formal definition for non-integer functional iteration, and presented both formulaic solutions for specific functions as well as numerical techniques for finding fairly accurate approximations of functional roots. Most of our formulas were able to be generalized to complex iteration of complex functions, which implies that it is often possible to extend iteration to complex functions and complex exponents.

    Future Work

    For formulaic solutions, we mostly used solutions in terms of arithmetic operations. However, using infinite sums or integrals may allow one to apply our techniques for solving functional iteration to a wider variety of function types. In particular, solving for a half iterate by means of a fourier series may be possible for many more functions, especially periodic functions. Promisingly, Newton Series can often be used to solve this type of problem, as is used in https://mathoverflow.net/questions/17605/how-to-solve-ffx-cosx/44727#44727.
    As for numerical approximations, although we found a variety of methods, the ones which yielded the best results relied heavily on random search and generalized function approximation (instead of iterate-specific improvements). This implies that more efficient algorithms such as gradient descent, simulated annealing, etc. may converge much more quickly on solutions. A more complete and general method of approximation may use one of these algorithms in conjunction with other optimization methods to achieve a maximally accurate numerical approximation.

    Citations

    1
    .
    Curtright, T. Functional roots of the sine function. Notes. (2011).
    2
    .
    Erdős, P., and Jabotinsky, E. On analytic iteration. J. Analyse Math. 8 (1960/61), 361–376.
    3
    .
    Szekeres, G. Fractional iteration of exponentially growing functions. J. Austral. Math. Soc. 2 (1961/62), 301–320.
    4
    .
    Curtright, Thomas, Xiang Jin, and Cosmas Zachos. “Approximate Solutions of Functional Equations.” Journal of physics. A, Mathematical and theoretical 44.40 (2011): 405205–12. Web.
    5
    .
    Chen, T., & Scheinerman, E. (2022). Finding a Compositional Square Root of Sine. The American Mathematical Monthly, 129(9), 816–830. https://doi.org/10.1080/00029890.2022.2105054

    Acknowledgements

    Lastly, we’d like to thank a couple people who made our WELP project possible. Many thanks to Eryn Gillam and Rory Foulger, who make this program possible. We’d also like to thank our mentor Daniele Ceravolo for everything he helped us with along the way. Finally, thanks to every YouTube Video, textbook, and teacher who gave us the math and the computer skills necessary to get here.

    CITE THIS NOTEBOOK

    On functional iteration and roots​
    by Wolfram Education Programs​
    Wolfram Community, STAFF PICKS, February 4, 2025
    ​https://community.wolfram.com/groups/-/m/t/3379973