【问题标题】:Python scipy fsolve works incorrectlyPython scipy fsolve 工作不正确
【发布时间】:2015-12-01 21:09:11
【问题描述】:

我有方程式:

import numpy as np
from scipy import optimize

def wealth_evolution(price, wealth=10, rate=0.01, q=1, realEstate=0.1, prev_price=56):
    sum_wantedEstate = 100
    for delta in range(1,4):
        z = rate - ((price-prev_price) / (price + q / rate))
        k = delta * np.divide(1.0, float(np.maximum(0.0, z)))
        wantedEstate = (wealth / (price + q / rate)) * np.minimum(k, 1) - realEstate
        sum_wantedEstate += wantedEstate
    return sum_wantedEstate

所以我找到了这个方程的解:

sol = optimize.fsolve(wealth_evolution, 200)

但如果我将sol 代入等式,我将不会得到0 (welth_evolution(sol))。为什么会发生? fsolve 找到 f(x)=0 的根。

更新: full_output 给出:

(array([ 2585200.]), {'qtf': array([-99.70002298]), 'nfev': 14, 'fjac': array([[-1.]]), 'r': array([  3.45456519e-11]), 'fvec': array([ 99.7000116])}, 5, 'The iteration is not making good progress, as measured by the \n  improvement from the last ten iterations.')

【问题讨论】:

  • 你可能想显示minimal reproducible example
  • 如果您发布由full_info=True 生成的调试输出也会很有帮助。
  • @cel 实际上,它非常接近。运行代码,用wealth_evolution(sol)检查结果。
  • Roma,听听@ali_m 的建议(除了要使用的参数是full_output=True)。特别是查看ier 和错误消息。 (还要弄清楚为什么会收到关于被零除的警告。)当使用full_output=True 时,函数返回ier=5,这是一个错误条件。我不明白为什么该函数不会引发异常。
  • 其实第一次运行没有full_output=True的时候,代码会生成一个RuntimeWarning关于没有取得好的进展。这表明您应该使用full_output=True 重新运行代码。 (有点烦人的是,如果您以相同的初始猜测再次运行它,您不会收到该警告。这是一个示例,python 的默认警告行为可能具有欺骗性,scipy 可能应该覆盖它并始终生成警告.)

标签: python scipy equation


【解决方案1】:

你尝试过绘制你的函数吗?

import numpy as np
from scipy import optimize
from matplotlib import pyplot as plt
small = 1e-30
def wealth_evolution(price, wealth=10, rate=0.01, q=1, realEstate=0.1, prev_price=56):
    sum_wantedEstate = 100
    for delta in range(1,4):
        z = rate - ((price-prev_price) / (price + q / rate))
        k = delta * np.divide(1.0, float(np.maximum(small, z)))
        wantedEstate = (wealth / (price + q / rate)) * np.minimum(k, 1) - realEstate
        sum_wantedEstate += wantedEstate
    return sum_wantedEstate




price_range = np.linspace(0,10000,10000)
we = [wealth_evolution(p) for p in price_range]

plt.plot(price_range,we)
plt.xlabel('price')
plt.ylabel('wealth_evolution(price)')
plt.show()

至少对于您指定的参数,它没有根,这是fsolve 试图找到的。如果你想最小化一个函数,你可以试试fmin。对于这个函数,这将无济于事,因为它似乎只是渐近衰减到 99.7 左右。所以最小化会导致无限的价格

所以要么你必须忍受这个想出一个不同的功能来优化限制你的搜索范围(在这种情况下你不必搜索,因为它只是最大值......)。

【讨论】:

  • 人们总是让我感到惊讶,人们使用 Python 之类的不错的交互式工具,然后期望它神奇地找到他们“知道”存在的东西,而无需检查(通过绘图)它是否确实存在。你应该总是绘制你的函数!这应该为您提供关于使用哪个求解器(或最小化器)的非常好的信息,以及关于从哪里开始搜索的好主意。与在卡上运行 Fortran 程序的日子相比(是的,我这样做了),现在很容易做到。
猜你喜欢
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2013-12-03
  • 1970-01-01
  • 1970-01-01
  • 2015-08-18
  • 1970-01-01
  • 2018-03-25
相关资源
最近更新 更多