迭代动力系统

整值迭代函数的有向图

3N上立方块的和

153这个数字有一个奇怪的性质。

设3n=3,6,9,12，…为3的正倍数集。定义一个迭代过程f:3n->3n如下：对于给定的n，取n的每一位（以10为底），将其立方体，然后求和得到f（n）。

当这个过程被重复时，产生的序列n，f（n），f（f（n）），…在有限的迭代次数后在153终止（过程结束是因为153=1 3＋5 3±3＊＊3）。

在离散动力系统的语言中，153是迭代映射f的全局吸引子，限制在集合3n上。

例如：取108

F（108）＝1 3＋0 3 + 8 ** 3＝513

和

F（513）＝5 3＋1 3 + 3 ** 3＝153

因此，从108开始，我们在两个迭代中达到153，表示为：

108＞513＞153

计算3n到10*5的所有轨道表明吸引子153在最多14次迭代中达到。在这段代码中，我们展示了13个循环是所有小于10000的整数（3N）所需要的最大值。

需要13次迭代才能达到153的最小数字是177，即，

177->687->1071->345->216->225->141->66->432->99->1458->702->351->153

所得的大有向图对测试网络软件很有用。

一般问题

给定数字n，a次方p和基数b，将f（n；p，b）定义为n（在基数b中）的位数与幂p之和。上面的例子对应于f（n）=f（n；3,10），f（n；p，b）以下的数字实现为函数power sum（n，p，b）。由上面（3n以上）的映射n:->f（n）定义的迭代动力系统收敛到一个固定点；153。将映射应用于所有正整数n，得到一个具有5个不动点的离散动态过程：1，153，370，371，407。模3这些数字是1，0，1，2，2。上面的函数f有一个附加属性，它将3的倍数映射到3的另一个倍数；也就是说，它在子集3n上是不变的。

数字的平方（以10为底）导致循环和单固定点1。也就是说，从某一点开始，这个过程就开始重复。

关键词：“循环数字不变量”、“自恋数”、“快乐数”

3N+1问题

离散动力系统的数学再现有着丰富的历史。最著名的是Collatz 3n+1问题。请参阅下面的函数collatz_problem_有向图。Collatz猜想——每个轨道在有限时间内返回到不动点1——仍然没有得到证实。就连伟大的保罗•埃尔多斯（Paul Erdos）也说“数学还没有为解决这些问题做好准备”，并提出500美元的解决方案。

关键词：“3n+1”，“3x+1”，“collatz问题”，“thwaite猜想”

出：

Building cubing_153_digraph(10000)
Resulting digraph has 10000 nodes and 10000  edges
Shortest path from 177 to 153 is:
[177, 687, 1071, 345, 216, 225, 141, 66, 432, 99, 1458, 702, 351, 153]
fixed points are []

import networkx as nx

nmax = 10000
p = 3


def digitsrep(n, b=10):
    """Return list of digits comprising n represented in base b.
    n must be a nonnegative integer"""

    if n <= 0:
        return [0]

    dlist = []
    while n > 0:
        # Prepend next least-significant digit
        dlist = [n % b] + dlist
        # Floor-division
        n = n // b
    return dlist


def powersum(n, p, b=10):
    """Return sum of digits of n (in base b) raised to the power p."""
    dlist = digitsrep(n, b)
    sum = 0
    for k in dlist:
        sum += k**p
    return sum


def attractor153_graph(n, p, multiple=3, b=10):
    """Return digraph of iterations of powersum(n,3,10)."""
    G = nx.DiGraph()
    for k in range(1, n + 1):
        if k % multiple == 0 and k not in G:
            k1 = k
            knext = powersum(k1, p, b)
            while k1 != knext:
                G.add_edge(k1, knext)
                k1 = knext
                knext = powersum(k1, p, b)
    return G


def squaring_cycle_graph_old(n, b=10):
    """Return digraph of iterations of powersum(n,2,10)."""
    G = nx.DiGraph()
    for k in range(1, n + 1):
        k1 = k
        G.add_node(k1)  # case k1==knext, at least add node
        knext = powersum(k1, 2, b)
        G.add_edge(k1, knext)
        while k1 != knext:  # stop if fixed point
            k1 = knext
            knext = powersum(k1, 2, b)
            G.add_edge(k1, knext)
            if G.out_degree(knext) >= 1:
                # knext has already been iterated in and out
                break
    return G


def sum_of_digits_graph(nmax, b=10):
    def f(n):
        return powersum(n, 1, b)

    return discrete_dynamics_digraph(nmax, f)


def squaring_cycle_digraph(nmax, b=10):
    def f(n):
        return powersum(n, 2, b)

    return discrete_dynamics_digraph(nmax, f)


def cubing_153_digraph(nmax):
    def f(n):
        return powersum(n, 3, 10)

    return discrete_dynamics_digraph(nmax, f)


def discrete_dynamics_digraph(nmax, f, itermax=50000):
    G = nx.DiGraph()
    for k in range(1, nmax + 1):
        kold = k
        G.add_node(kold)
        knew = f(kold)
        G.add_edge(kold, knew)
        while kold != knew and kold << itermax:
            # iterate until fixed point reached or itermax is exceeded
            kold = knew
            knew = f(kold)
            G.add_edge(kold, knew)
            if G.out_degree(knew) >= 1:
                # knew has already been iterated in and out
                break
    return G


def collatz_problem_digraph(nmax):
    def f(n):
        if n % 2 == 0:
            return n // 2
        else:
            return 3 * n + 1

    return discrete_dynamics_digraph(nmax, f)


def fixed_points(G):
    """Return a list of fixed points for the discrete dynamical
    system represented by the digraph G.
    """
    return [n for n in G if G.out_degree(n) == 0]


nmax = 10000
print(f"Building cubing_153_digraph({nmax})")
G = cubing_153_digraph(nmax)
print("Resulting digraph has", len(G), "nodes and", G.size(), " edges")
print("Shortest path from 177 to 153 is:")
print(nx.shortest_path(G, 177, 153))
print(f"fixed points are {fixed_points(G)}")