scipy 如何在Python中构建扩展数学表达式?

l5tcr1uw  于 2023-02-22  发布在  Python
关注(0)|答案(4)|浏览(195)

我有一个问题,我需要求解一个次数随着每次迭代(循环)而增加的多项式,表达式为
(1/(1+x)^1)+(1/(1+x)^2)+(1/(1+x)^3)+(1/(1+x)^4)+(1/(1+x)^5)
正如你所看到的,指数随着每次迭代而增长,这个表达式将被用在函数“fsolve”中。我不想计算这个表达式,而是构建它以便在函数“fsolve”中使用它。(我正在使用scipy.optimize导入fsolve)提前感谢
编辑

def func(x):
    return["Here i will have my mathematical expression in the question"]
print(fsolve(func, [0.05]))

其中0.05是当等于零时表达式的解的起始猜测。
编辑2
整个问题是:

C1 - (1-1/(1+x)^n - C2*sum[1/(1+x) + 1/(1+x)^2)+ ... + 1/(1+x)^n]

其中C1C2是两个常数,n是迭代次数The boundaries of the sum-expression is from 1 to n

eiee3dmh

eiee3dmh1#

我想用一个理解:

exp = " + ".join([f"(1/(1+x)^{i})" for i in range(1, 6)])
print(exp)

为您提供:

(1/(1+x)^1) + (1/(1+x)^2) + (1/(1+x)^3) + (1/(1+x)^4) + (1/(1+x)^5)
ykejflvf

ykejflvf2#

Scipy提供了期望数值函数作为参数的数值求解器。没有必要编写一个返回函数表示的字符串的函数,那将不起作用。
实现您所要求的一个简单方法是使用factory(这里我们将使用decorated function):

import numpy as np
from scipy import optimize
    
def factory(order=1):
    @np.vectorize
    def wrapped(x):
        return np.sum([1/np.power(1 + x, i + 1) for i in range(order)])
    return wrapped

关键思想是将order参数传递给装饰器(factory),装饰器将返回带有fsolve兼容签名的装饰函数(wrapped)。

  • wrapped函数是您的目标函数,它返回所需数量wrt orderx
  • factory函数是一种常用的order参数化方法,每次您想要研究特定的order时,都可以避免重写整个过程,此函数返回wrapped修饰函数。

@np.vectorize decorator是一个numpy的矢量化工具,它可以对数组元素进行操作,并返回每个值的结果。

x = np.linspace(-20, 20, 2000)
fig, axe = plt.subplots()
for order in range(1, 6):
    y = factory(order)(x) # This call will return 1 value if not vectorized, returns 2000 values as expected when vectorized
    axe.plot(x, y, label="Order %d" % order)

然后简单地迭代您的订单,以获得您想要求解的函数:

for order in range(1, 11):
    function = factory(order=order)
    solution = optimize.fsolve(function, x0=-2.1)
    print(order, solution)

# 1 [-1.93626047e+83]
# 2 [-2.]
# 3 [-1.64256077e+83]
# 4 [-2.]
# 5 [-1.47348643e+83]
# 6 [-2.]
# 7 [-1.42261052e+83]
# 8 [-2.]
# 9 [-1.43409773e+83]
# 10 [-2.]

实际上,由于fsolve正在寻找根,而在您的设置中有没有真实的根的订单。Analyzing the function of order 5将显示there is no real roots,但当x变为无穷大(不是根)时,我们仍然在实域中。对于偶数订单,根很容易通过上述过程找到。

因此,您可能会有float overflow错误,有时高数字是毫无意义的w.r.t. a root finding problem为奇数订单。

OP更新

基于您在文章中详细介绍的新函数,您可以向目标函数添加额外的参数,如下所示:

def factory(order=1):
    @np.vectorize
    def wrapped(x, C1, C2):
        return C1 - (1 - 1/np.power(1 + x, order) - C2*np.sum([1/np.power(1 + x, i + 1) for i in range(order)]))
    return wrapped

然后我们像往常一样求解几个订单,我选择了(C1, C2) = (10,10)

for order in range(1, 11):
    function = factory(order=order)
    solution = optimize.fsolve(function, x0=-2.1, args=(10, 10))
    print(order, solution)

# 1 [-2.22222222]
# 2 [-3.20176167]
# 3 [-2.10582151]
# 4 [-2.73423813]
# 5 [-2.06950043]
# 6 [-2.54187609]
# 7 [-2.0517491]
# 8 [-2.43340793]
# 9 [-2.04122279]
# 10 [-2.36505166]
zyfwsgd6

zyfwsgd63#

如果你可以放松对fsolve的使用,你也可以简单地通过使用numpy.polynomials和做一些代数运算来解决复数的这个问题。
首先认识到根的职能是根的分子是一个多项式:

def numerator(order=1):
    term = Polynomial([1, 1])  # m(x) = x + 1
    return sum([term**i for i in range(order)])

然后只需对您要分析的每个订单应用roots方法:

for i in range(1, 11):
    print(i, numerator(order=i).roots())

# 1 []
# 2 [-2.]
# 3 [-1.5-0.8660254j -1.5+0.8660254j]
# 4 [-2.+0.j -1.-1.j -1.+1.j]
# 5 [-1.80901699-0.58778525j -1.80901699+0.58778525j -0.69098301-0.95105652j
#    -0.69098301+0.95105652j]
# 6 [-2. +0.j        -1.5-0.8660254j -1.5+0.8660254j -0.5-0.8660254j
#    -0.5+0.8660254j]
# 7 [-1.90096887-0.43388374j -1.90096887+0.43388374j -1.22252093-0.97492791j
#     -1.22252093+0.97492791j -0.3765102 -0.78183148j -0.3765102 +0.78183148j]
# 8 [-2.        +0.j         -1.70710678-0.70710678j -1.70710678+0.70710678j
#    -1.        -1.j         -1.        +1.j         -0.29289322-0.70710678j
#    -0.29289322+0.70710678j]
# 9 [-1.93969262-0.34202014j -1.93969262+0.34202014j -1.5       -0.8660254j
#    -1.5       +0.8660254j  -0.82635182-0.98480775j -0.82635182+0.98480775j
#    -0.23395556-0.64278761j -0.23395556+0.64278761j]
# 10 [-2.        +0.j         -1.80901699-0.58778525j    -1.80901699+0.58778525j
#     -1.30901699-0.95105652j -1.30901699+0.95105652j -0.69098301-0.95105652j
#     -0.69098301+0.95105652j -0.19098301-0.58778525j -0.19098301+0.58778525j]
dtcbnfnu

dtcbnfnu4#

你确定你只需要五项吗?从分析上讲,如果有无穷多项,你的整个表达式的计算结果是1 + 1/x,这是平凡可逆的。

相关问题