迭代动力系统

来自整数值迭代函数的有向图

在 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 次迭代后达到。在这段代码中我们展示了对于所有小于 10,000 的整数（在 3N 中），最多需要 13 次迭代。

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

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

由此产生的大型有向图对于测试网络软件很有用。

一般问题

给定数字 n，幂 p 和基数 b，将 F(n; p, b) 定义为 n 的数字（在基数 b 下）的 p 次幂的总和。上面的例子对应于 f(n)=F(n; 3,10)，下面将 F(n; p, b) 实现为函数 powersum(n,p,b)。上面由映射 n:->f(n) 定义的迭代动力系统（在 3N 上）收敛到一个不动点；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_digraph。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)}")