Python中具有分数幂的二项式展开答案

【问题标题】：Binomial expansion with fractional powers in PythonPython中具有分数幂的二项式展开
【发布时间】：2023-03-22 14:15:02
【问题描述】：

在 Scypy/numpy 中有没有一种快速的方法来扩展和求解二项式的分数幂？

例如，我想解下列方程

y * (1 + x)^4.8 = x^4.5

其中 y 是已知的（例如 1.03）。

这需要 (1 + x)^4.8 的二项式展开。

我希望为数百万个 y 值执行此操作，因此我正在寻找一种快速且快速的方法来解决此问题。

我尝试了 sympy 扩展（和简化），但它似乎不喜欢小数指数。我也在为 scipy fsolve 模块苦苦挣扎。

任何正确方向的指针将不胜感激。

编辑：

到目前为止，我发现的最简单的解决方案是为 x（和已知 y）的假定值生成一个真值表 (https://en.wikipedia.org/wiki/Truth_table)。这允许快速插值“真实”x 值。

y_true = np.linspace(7,12, 1e6)
x = np.linspace(10,15, 1e6)

a = 4.5
b = 4.8
y = x**(a+b) / (1 + x)**b

x_true = np.interp(y_true, y, x)

【问题讨论】：

标签： python math numpy scipy polynomial-math

【解决方案1】：

编辑：将输出与 y=1.03 的 Woldfram alpha 的输出进行比较后，看起来 fsolve 不会返回复数根。 https://stackoverflow.com/a/15213699/3456127 是一个类似的问题，可能会有更多帮助。

重新排列你的方程式：y = x^4.5 / (1+x)^4.8。 Scipy.optimize.fsolve() 需要一个函数作为其第一个参数。

要么：

from scipy.optimize import fsolve
import math
def theFunction(x):   
    return math.pow(x, 4.5) / math.pow( (1+x) , 4.8)

for y in millions_of_values:
    fsolve(theFunction, y)

或者使用lambda（匿名函数构造）：

from scipy.optimize import fsolve
import math
for y in millions_of_values:
    fsolve((lambda x: (math.pow(x, 4.5) / math.pow((1+x), 4.8))), y)

【讨论】：

我认为你需要为函数提供y；否则，它将如何解决y 的各种值？
我定义的函数不能求解 y，因为 y 是一个已知值，我们希望发现 x 的值，这是 fsolve 应该做的，但由于这个方程有复根，所以不能.示例：给定1.03 = x^4.5 / (1+x)^4.8 求解 x。 Wolfram alpha 显示 x 有 4 个复数值。
是的，但正如所写，您的 theFunction 将始终返回相同的值，无论 y 是什么。 fsolve 的语法将使用 y 作为根的 x 值的初始猜测，而不是作为函数的附加参数。
不幸的是，我的微积分能力没有我希望的那么强。我看到你在下面的答案中有return y * (1 + x)**(4.8) / x**(4.5) - 1。我会做return y * (1+x)**4.8 - x**4.5，但我不认为我是正确的。
不，这实际上也是正确的，我在我的版本中尝试过，但没有成功。我认为这与数值稳定性有关，尽管我没有对此进行验证。如果可以的话，最好将两个大数相除而不是相减。

【解决方案2】：

我认为您不需要二项式展开。 Horner's method 用于评估多项式意味着多项式的因式分解比展开式更好。

一般来说，非线性方程求解可以从符号微分中受益，这对于您的方程而言并不难手动完成。为导数提供解析表达式可以使求解器不必对导数进行数值估计。您可以编写两个函数：一个返回函数的值，另一个返回导数（即这个简单一维函数的函数Jacobian），如the docs for scipy.optimize.fsolve() 中所述。一些采用这种方法的代码：

import timeit
import numpy as np
from scipy.optimize import fsolve

def the_function(x, y): 
    return y * (1 + x)**(4.8)   /   x**(4.5) - 1

def the_derivative(x, y):
    l_dh = x**(4.5) * (4.8 * y * (1 + x)**(3.8))
    h_dl = y * (1 + x)**(4.8) * 4.5 * x**3.5
    square_of_whats_below = x**9
    return (l_dh - h_dl)/square_of_whats_below

print fsolve(the_function, x0=1, args=(0.1,))
print '\n\n'
print fsolve(the_function, x0=1, args=(0.1,), fprime=the_derivative)

%timeit fsolve(the_function, x0=1, args=(0.1,))
%timeit fsolve(the_function, x0=1, args=(0.1,), fprime=the_derivative)

...给我这个输出：

[ 1.79308495]



[ 1.79308495]
10000 loops, best of 3: 105 µs per loop
10000 loops, best of 3: 136 µs per loop

这表明在这种特殊情况下，分析微分不会导致任何加速。我的猜测是，该函数的数值逼近涉及更易于计算的函数，例如乘法、平方和/或加法，而不是分数取幂之类的函数。

您可以通过记录方程的对数并绘制它来获得额外的简化。通过一点代数，您应该能够获得ln_y 的显式函数，即y 的自然对数。如果我正确地完成了代数：

def ln_y(x):
    return 4.5 * np.log(x/(1.+x)) - 0.3 * np.log(1.+x)

你可以绘制这个函数，我已经为 lin-lin 和 log-log 绘制了这个函数：

%matplotlib inline
import matplotlib.pyplot as plt

x_axis = np.linspace(1, 100, num=2000)

f, ax = plt.subplots(1, 2, figsize=(8, 4))
ln_y_axis = ln_y(x_axis)

ax[0].plot(x_axis, np.exp(ln_y_axis))  # plotting y vs. x
ax[1].plot(np.log(x_axis), ln_y_axis)  # plotting ln(y) vs. ln(x)

这表明只要y 低于临界值，每个y 就有两个x 值。 y 的最小奇异值出现在 x=ln(15) 和 y 值为：

np.exp(ln_y(15))
0.32556278053267873

因此，您的示例 y 的 1.03 值导致 x 没有（真正的）解决方案。

我们从情节中辨别出的这种行为被我们之前进行的scipy.optimize.fsolve() 调用所概括：

print fsolve(the_function, x0=1, args=(0.32556278053267873,), fprime=the_derivative)
[ 14.99999914]

这表明最初猜测x=1，当y 是0.32556278053267873 时，得到x=15 作为解决方案。尝试更大的y 值：

print fsolve(the_function, x0=15, args=(0.35,), fprime=the_derivative)

导致错误：

/Users/curt/anaconda/lib/python2.7/site-packages/IPython/kernel/__main__.py:5: RuntimeWarning: invalid value encountered in power

错误的原因是 Python（或 numpy）中的 power 函数默认情况下不接受小数指数的负基数。你可以通过提供复数的幂来解决这个问题，即写x**(4.5+0j)而不是x**4.5，但是你真的对可以解决你的方程的复杂x值感兴趣吗？

【讨论】：