Python中分数次幂的二项式展开

Question

在Scypy/numpy中，有没有一种快速的方法可以将二元数展开并求解成分数次幂？在

例如，我想解下面的方程

y*（1+x）^4.8=x^4.5

其中y已知（例如1.03）。在

这需要（1+x）^4.8的二项式展开。在

我希望对数百万个y值做这个，所以我在寻找一个好的和快速的方法来解决这个问题。在

我尝试过sympy展开（和简化），但它似乎不喜欢分数指数。我也在和scipy fsolve模块做斗争。在

任何指向正确方向的建议都将不胜感激。在

编辑：

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

Answer 1

我认为你不需要二项式展开。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)

…给我这个输出：

^{pr2}$

这表明，在这种特殊情况下，分析微分并没有导致任何加速。我的猜测是，对函数的数值逼近涉及到更容易计算的函数，如乘法、平方和/或加法，而不是像分数幂这样的函数。在

你可以通过记录你的方程并绘制它来获得额外的简化。用一点代数，你应该能够得到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低于临界值，每个x都有两个值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是{}，则得到{}作为解。尝试更大的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值感兴趣吗？在

Answer 2

编辑：将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（匿名函数构造）：

^{pr2}$