使用 Python for 循环用数值方法求解方程

Question

我正在学习 Python 和概率。

我有这样的表达：

c = n! /((n-1)!1!) + 2*n! /((n-2)*2!) + 3*n!/((n-3)*3!) + ...+ n*n!/((n-n)*n!

0! = 1 和 !表示“阶乘”，即 n! = 1*2*3*...*(n-1)*n。 等于：

(a + n^b)*2^(c*n + d).  (^ signifies exponent)

我的目标是使用“蛮力”确定参数 a、b、c、d。

我使用上面的公式 c 计算了 n=3 (12), n=4 (38), n=5 (80), n=6 (192), n=7 (448)。

然后我将参数表示为两个整数的比率：即 a = a1/a2, b=b1/b2, c=c1/c2, d=d1/d2。

最后我定义了以下函数：

def com():
    parms = []
    for a1 in range(-10, 10):
        for a2 in range(1,11):
            for b1 in range(-10,10):
                for b2 in range(1, 11):
                    for c1 in range(-10,10):
                        for c2 in range(1, 11):
                            for d1 in range(-10 , 10):
                                for d2 in range(1 , 11):
                                    a = a1/a2
                                    b = b1/b2
                                    c = c1/c2
                                    d = d1/d2
                                    cr1 = ( 12 == (a + 3**b)*2**(c*3+d) )
                                    cr2 = ( 38 == (a + 4**b)*2**(c*4+d) )
                                    cr3 = ( 80 == (a + 5**b)*2**(c*5+d) )
                                    cr4 = ( 192 == (a + 6**b)*2**(c*6+d) )
                                    cr5 = ( 448 == (a + 7**b)*2**(c*7+d) )
                                    criterion = cr1 & cr2 & cr3 & cr4 & cr5
                                    if criterion == 1 :
                                        parms = [a, b, c, d]
                                        break
    return parms

但是，我的函数返回一个空列表。

你能解释一下吗？您对如何实现我的目标有什么建议吗？

您的建议将不胜感激。

Answer 1

首先，让我们看一下代码中的一些错误：

criterion = cr1 & cr2 & cr3 & cr4 & cr5

这是对值执行bitwise 操作。你可能想要：

criterion = cr1 and cr2 and cr3 and cr4 and cr5

Python 还提供了一个all 函数，你可以检查all 是否为True：

criterion = all([cr1,cr2,cr3,cr4,cr5])

现在，让我们看看你的if 声明：

if criterion == 1 :

既然您想知道criterion 是True 还是False，您可以简单地使用if criterion:。

最后，这种方法不太可能是True。以这一行为例：

cr1 = ( 12 == (a + 3**b)*2**(c*3+d) )

这些数字加起来必须正好是 12，否则会是 False，然后你会得到一个空列表。

另外，computers can't do decimal maths accurately。

要解决实际值，您必须使用参数替换。这不是一个编程问题，但数学和编程很好地结合在一起，所以这里是一个初学者：

12 = (a + 3**b)*2**(c*3+d)
12/(a+3**b) = 2**(c*3+d)

... 等等。用b、c 和d 得到a，然后用你的cr2 替换你拥有的数字a，并用c 得到b 和d.

再重复几次，你就有了四个值的数字。

Answer 2

您可以使用scipy curvefit 并使用非线性最小二乘法将函数 f 拟合到数据并找到参数 (a,b,c,d) 值。您提供的示例数据更适合函数(a + n**b)+2**(c*n + d)。请注意，最适合 a,b,c,d 值的是浮点数，而不是整数。

import numpy as np
from scipy.optimize import curve_fit

xdata = np.array(range(3,8))
ydata = np.array([12,38,80,192,448])

def func(n, a, b, c, d):
  return (a + n**b)+2**(c*n + d)

popt, pcov = curve_fit(func, xdata, ydata, p0=(1,1,1,1))
a, b, c, d = popt
print a, b, c, d # learnt parameters
# -5.62374782967 1.79345876905 1.29232902743 -0.328778229316

import matplotlib.pyplot as plt
plt.scatter(xdata, ydata)
plt.plot(xdata, func(xdata, a, b, c, d), '-r', label='np.poly')
plt.show()

在上面的拟合曲线中，蓝点是数据点，红线是带有学习参数的拟合函数，注意该函数是连续的，并针对n的所有值定义。

Answer 3

用大脑代替 Python：

您识别出修改后的二项式和，Sum k.Cnk 而不是 Sum Cnk。

然后注意到

k.n!/(k!(n-k)!) = n!/((k-1)!(n-k)!) = n.(n-1)!/(k-1)!(n-1-k+1)!) = n.Cn-1,k+1.

所以你的总和是n.2^(n-1)。

Answer 4

在急于使用蛮力之前，尝试估计参数的狭窄范围可能更明智。

总和的第一个值是

1, 4, 12, 32, 80, 192, 448, 1024, 2304, 5120, 11264

采用成对比率，您会发现它们越来越接近 2。

由于比率 (a + (n+1)^b)/(a + n^b) 将趋于 1，而比率 2^(c(n+1)+d)/2^(cn+d) 将趋于 2^c，因此提示您 c 可能接近 1。您可以通过尝试越来越大的值来检查此估计值n 的比率（例如，n=1000 产生比率 2.002002... = 2^1.00144...）。

然后为了降低表达式的增长率使其更易于处理，查看Sum/2^n 的值可能会很有趣，我们得到

1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11

哼哼。

---

不幸的是，这种情况对于正确演示该方法来说太容易了。总体思路是通过观察渐近行为来粗略估计一些参数。而当你有了这样一个粗略的估计（在一个很小的范围内），那么你可以通过某种方式取消它的效果，以便更好地看到其他的效果。

对于详尽的搜索，了解范围很重要。