在 Python 中运行求解方程函数

Question

我正在尝试运行此函数以求解 f(f(....(f(x))) = x 形式的方程，这是我的 python 代码：

from sympy import * 
x=symbols('x')

def give_list(x1):
    P,Q=x**2-1,x
    l=[]
    for k in range(x1):
        print(x1)
        U=P/Q
        P,Q=P.simplify(),Q.simplify()
        eqn = U - 1/U - x
        l+=(Poly(eqn.as_numer_denom()[0]).nroots(n=4))
        P,Q=P**2-Q**2,P*Q
        P,Q=P.simplify(),Q.simplify()
    return l
print(len(give_list(24)))

我收到了这个错误：

24
24
24
24
24
24
Traceback (most recent call last):
  File "C:\Users\MSI\AppData\Local\Programs\Python\Python39\lib\site-packages\sympy\polys\polytools.py", line 3649, in nroots
    roots = mpmath.polyroots(coeffs, maxsteps=maxsteps,
  File "C:\Users\MSI\AppData\Local\Programs\Python\Python39\lib\site-packages\mpmath\calculus\polynomials.py", line 195, in polyroots
    raise ctx.NoConvergence("Didn't converge in maxsteps=%d steps." \
mpmath.libmp.libhyper.NoConvergence: Didn't converge in maxsteps=50 steps.

在处理上述异常的过程中，又发生了一个异常：

Traceback (most recent call last):
  File "C:\Users\MSI\AppData\Local\Programs\Python\Python39\PE729.py", line 66, in <module>
    print(len(give_list(24)))
  File "C:\Users\MSI\AppData\Local\Programs\Python\Python39\PE729.py", line 62, in give_list
    l+=(Poly(eqn.as_numer_denom()[0]).nroots(n=4))
  File "C:\Users\MSI\AppData\Local\Programs\Python\Python39\lib\site-packages\sympy\polys\polytools.py", line 3657, in nroots
    raise NoConvergence(
mpmath.libmp.libhyper.NoConvergence: convergence to root failed; try n < 4 or maxsteps > 50

注意：我不能再减少 n 了，是否有另一种解决方案可以使用 x1 = 24 运行我的函数 give_list

Answer 1

这解决了f(f(....(f(x))) = x 为f(x) = x-1/x 的数字，正如您在 cmets 中所说的那样，这就是您正在尝试做的事情。

版本 1：保证找到所有解决方案，但速度相对较慢。

from sympy import * 

x=symbols('x')

def repeat(func, n):
    def new_func(x):
        for _ in range(n):
            x = func(x)
        return x
    return new_func


def f(x):
    return x-1/x 

for n in range(2, 7):
    eqn = repeat(f,n)(x)-x
    print(n, Poly(eqn.as_numer_denom()[0]).nroots())

这对我有用并给出：

2 [-0.707106781186548, 0.707106781186548]
3 [-1.46190220008154, -0.777861913430206, -0.507713305942872, 0.507713305942872, 0.777861913430206, 1.46190220008154]
4 [-1.98904379073655, -1.48628965095479, -1.30656296487638, -0.813473286151600, -0.707106781186548, -0.541196100146197, -0.415823381635519, 0.415823381635519, 0.541196100146197, 0.707106781186548, 0.813473286151600, 1.30656296487638, 1.48628965095479, 1.98904379073655]
5 [-2.41545660696940, -2.00145620757714, -1.85510188820248, -1.50181999459672, -1.41914220567326, -1.31604793846555, -1.24172627706540, -0.835961234160827, -0.774521517521973, -0.714491187612963, -0.685106200482040, -0.556197198403633, -0.516598196232952, -0.436395812155433, -0.360266484465194, 0.360266484465194, 0.436395812155433, 0.516598196232952, 0.556197198403633, 0.685106200482040, 0.714491187612963, 0.774521517521973, 0.835961234160827, 1.24172627706540, 1.31604793846555, 1.41914220567326, 1.50181999459672, 1.85510188820248, 2.00145620757714, 2.41545660696940]
6 [-2.78239673626141, -2.42299435953793, -2.29516266137635, -2.01028188414005, -1.93777043407717, -1.85946369466332, -1.80193773580484, -1.51283920812063, -1.46190220008154, -1.42171343247665, -1.40177496925008, -1.32167422188684, -1.29342081368380, -1.24697960371747, -1.20467134431239, -0.851830427655270, -0.810422925121491, -0.777861913430206, -0.764263390941947, -0.718336804559465, -0.707106781186548, -0.688393704827180, -0.673767837236344, -0.565058118281228, -0.544186041350644, -0.520277232398841, -0.507713305942872, -0.445041867912629, -0.423500707838532, -0.374569420894612, -0.322112140647414, 0.322112140647414, 0.374569420894612, 0.423500707838532, 0.445041867912629, 0.507713305942872, 0.520277232398841, 0.544186041350644, 0.565058118281228, 0.673767837236344, 0.688393704827180, 0.707106781186548, 0.718336804559465, 0.764263390941947, 0.777861913430206, 0.810422925121491, 0.851830427655270, 1.20467134431239, 1.24697960371747, 1.29342081368380, 1.32167422188684, 1.40177496925008, 1.42171343247665, 1.46190220008154, 1.51283920812063, 1.80193773580484, 1.85946369466332, 1.93777043407717, 2.01028188414005, 2.29516266137635, 2.42299435953793, 2.78239673626141]

第 2 版：相对来说要快得多，但可能会错过解决方案。

方法是将搜索区间划分为 10^-6 的步长，并检查 f_n(x)-x 的符号是否在这些步长期间发生了变化。如果是这样，请将其提供给找到零的求解器。

这种方法可能出现的一个问题是f_n(x)-x 有一个零但没有符号变化。这可能会发生，例如f(x)=x^2 的零，但不太可能。另一个是它确实改变了迹象，但我们必须检查更小的间隔来检测它。我选择 (10^-6, 20) 作为搜索间隔，因为考虑到 f_1,...f_6 的解决方案，20 似乎是一个合理的上限，并且由于 f_n-x 很奇怪，我们不需要担心负面解决方案（它们只是减去正解）。

希望第二个问题不会发生的一个原因是，这种方法确实找到了通过其他方法获得的 f_6 的所有正解。

fn = repeat(f,25)

epsilon = 1e-6
test_values = np.arange(0+epsilon,20,epsilon)
value_positive = fn(test_values) > 0
sign_change = (~value_positive[:-1])&(value_positive[1:])

@np.vectorize
def find_root(a,b):
    return optimize.root_scalar(fn, bracket=(a,b)).root

find_root(test_values[:-1][sign_change],test_values[1:][sign_change])

给予：

array([0.14548, 0.15015535, 0.15314301, ..., 6.52984421, 6.65976928,
       6.87379705])