numpy - Sympy 名称 Heaviside 未在 lambdifygenerated 中定义

Question

Heaviside 函数应该内置到 Sympy 和 Numpy 中，但是下面的代码给出了错误 Name Heaviside not defined。尝试在将使用它的数值计算(基于 Traceback)之前在代码中自己定义 Heaviside 函数什么也没做——我想它应该在 lambdifygenerated 中定义。有解决方法吗？

from sympy import *
from IPython.display import display
mux, s, Px, Py, Pxe, Pye = symbols("mu_X s P_X P_Y P_X^* P_Y^*", positive=True)
vx, vy, cx, cy = symbols("v_X v_Y c_X c_Y", real=True)
pix = (Px-cx)*( mux*integrate(integrate(1,(vx,Min(1,Max(0,Px+Max(0,vy-Pye-s))),1)),(vy,0,1))
      +(1-mux)*integrate(integrate(1,(vx,Min(1,Max(0,Max(Pxe+s,Px)+Max(0,vy-Pye))),1)),(vy,0,1))
     )
piy = (Py-cy)*( (1-mux)*integrate(integrate(1,(vy,Min(1,Max(0,Py+Max(0,vx-Pxe-s))),1)),(vx,0,1))
      +mux*integrate(integrate(1,(vy,Min(1,Max(0,Max(Pye+s,Py)+Max(0,vx-Pxe))),1)),(vx,0,1))
     )
focx =diff(pix,Px)
focy =diff(piy,Py)
focxeq=focx.subs(Px,Pxe)
focyeq=focy.subs(Py,Pye)

import numpy as np
focx_lambda = lambdify((Pxe,Pye), focxeq, modules=['numpy', 'sympy'])
focy_lambda = lambdify((Pxe,Pye), focyeq, modules=['numpy', 'sympy'])
nsolve([focxeq.subs({mux:0.4,s:0.05,cx:0,cy:0.1}).evalf(),focyeq.subs({mux:0.4,s:0.05,cx:0,cy:0.1}).evalf()],(Pxe,Pye),(0.3,0.4))

回溯如下:

--------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-10-b7bc7e96827d> in <module>
     26 focx_lambda = lambdify((Pxe,Pye), focxeq, modules=['numpy', 'sympy'])
     27 focy_lambda = lambdify((Pxe,Pye), focyeq, modules=['numpy', 'sympy'])
---> 28 nsolve([focxeq.subs({mux:0.4,s:0.05,cx:0,cy:0.1}).evalf(),focyeq.subs({mux:0.4,s:0.05,cx:0,cy:0.1}).evalf()],(Pxe,Pye),(0.3,0.4))
     29 mux=0.4
     30 s=0.05

~/anaconda3/lib/python3.6/site-packages/sympy/utilities/decorator.py in func_wrapper(*args, **kwargs)
     88         dps = mpmath.mp.dps
     89         try:
---> 90             return func(*args, **kwargs)
     91         finally:
     92             mpmath.mp.dps = dps

~/anaconda3/lib/python3.6/site-packages/sympy/solvers/solvers.py in nsolve(*args, **kwargs)
   3045     J = lambdify(fargs, J, modules)
   3046     # solve the system numerically
-> 3047     x = findroot(f, x0, J=J, **kwargs)
   3048     if as_dict:
   3049         return [dict(zip(fargs, [sympify(xi) for xi in x]))]

~/anaconda3/lib/python3.6/site-packages/mpmath/calculus/optimization.py in findroot(ctx, f, x0, solver, tol, verbose, verify, **kwargs)
    926         # detect multidimensional functions
    927         try:
--> 928             fx = f(*x0)
    929             multidimensional = isinstance(fx, (list, tuple, ctx.matrix))
    930         except TypeError:

<lambdifygenerated-23> in _lambdifygenerated(Dummy_4515, _Dummy_4514)
      1 def _lambdifygenerated(Dummy_4515, _Dummy_4514):
----> 2     return (ImmutableDenseMatrix([[Dummy_4515*(mpf((0, 3602879701896397, -53, 52))*((-(_Dummy_4514 - Dummy_4515 + mpf((0, 4728779608739021, -52, 53)))*Heaviside(1 - Dummy_4515)*Heaviside(1 - Max(_Dummy_4514 + mpf((0, 3602879701896397, -56, 52)), _Dummy_4514 - Dummy_4515 + mpf((0, 4728779608739021, -52, 53)))) + Heaviside(1 - Dummy_4515)*Heaviside(1 - Max(_Dummy_4514 + mpf((0, 3602879701896397, -56, 52)), _Dummy_4514 - Dummy_4515 + mpf((0, 4728779608739021, -52, 53))))*Min(mpf((0, 1, 0, 1)), Max(_Dummy_4514 + mpf((0, 3602879701896397, -56, 52)), _Dummy_4514 - Dummy_4515 + mpf((0, 4728779608739021, -52, 53)))) + Min(mpf((0, 1, 0, 1)), _Dummy_4514 + mpf((0, 3602879701896397, -56, 52))) - Min(mpf((0, 1, 0, 1)), Max(_Dummy_4514 + mpf((0, 3602879701896397, -56, 52)), _Dummy_4514 - Dummy_4515 + mpf((0, 4728779608739021, -52, 53))))) if (Dummy_4515 >= 1) else (-(_Dummy_4514 - Dummy_4515 + mpf((0, 4728779608739021, -52, 53)))*Heaviside(1 - Dummy_4515)*Heaviside(1 - Max(_Dummy_4514 + mpf((0, 3602879701896397, -56, 52)), _Dummy_4514 - Dummy_4515 + mpf((0, 4728779608739021, -52, 53)))) + Heaviside(1 - Dummy_4515)*Heaviside(1 - Max(_Dummy_4514 + mpf((0, 3602879701896397, -56, 52)), _Dummy_4514 - Dummy_4515 + mpf((0, 4728779608739021, -52, 53))))*Min(mpf((0, 1, 0, 1)), Max(_Dummy_4514 + mpf((0, 3602879701896397, -56, 52)), _Dummy_4514 - Dummy_4515 + mpf((0, 4728779608739021, -52, 53)))) - Min(mpf((0, 1, 0, 1)), Max(_Dummy_4514 + mpf((0, 3602879701896397, -56, 52)), _Dummy_4514 - Dummy_4515 + mpf((0, 4728779608739021, -52, 53)))))) + mpf((0, 5404319552844595, -53, 53))*((0) if (Dummy_4515 >= mpf((0, 4278419646001971, -52, 52))) else (-(_Dummy_4514 - Dummy_4515 + mpf((0, 4278419646001971, -52, 52)))*Heaviside(1 - Max(_Dummy_4514, _Dummy_4514 - Dummy_4515 + 1))*Heaviside(_Dummy_4514 - Dummy_4515 - Max(_Dummy_4514, _Dummy_4514 - Dummy_4515 + mpf((0, 4278419646001971, -52, 52))) + 1) + Heaviside(1 - Max(_Dummy_4514, _Dummy_4514 - Dummy_4515 + 1))*Heaviside(_Dummy_4514 - Dummy_4515 - Max(_Dummy_4514, _Dummy_4514 - Dummy_4515 + mpf((0, 4278419646001971, -52, 52))) + 1)*Min(mpf((0, 1, 0, 1)), Max(_Dummy_4514, _Dummy_4514 - Dummy_4515 + mpf((0, 1, 0, 1))))) if (Dummy_4515 >= 1) else (0))) + mpf((0, 3602879701896397, -53, 52))*((-(_Dummy_4514 - Dummy_4515 + mpf((0, 4728779608739021, -52, 53)))*Min(mpf((0, 1, 0, 1)), _Dummy_4514 + mpf((0, 3602879701896397, -56, 52))) + (_Dummy_4514 - Dummy_4515 + mpf((0, 4728779608739021, -52, 53)))*Min(mpf((0, 1, 0, 1)), Max(_Dummy_4514 + mpf((0, 3602879701896397, -56, 52)), _Dummy_4514 - Dummy_4515 + mpf((0, 4728779608739021, -52, 53)))) + mpf((0, 1, -1, 1))*Min(mpf((0, 1, 0, 1)), _Dummy_4514 + mpf((0, 3602879701896397, -56, 52)))**2 - mpf((0, 1, -1, 1))*Min(mpf((0, 1, 0, 1)), Max(_Dummy_4514 + mpf((0, 3602879701896397, -56, 52)), _Dummy_4514 - Dummy_4515 + mpf((0, 4728779608739021, -52, 53))))**2) if (Dummy_4515 >= 1) else ((mpf((0, 1, 0, 1)) - Dummy_4515)*Min(mpf((0, 1, 0, 1)), _Dummy_4514 + mpf((0, 3602879701896397, -56, 52))) - (_Dummy_4514 - Dummy_4515 + mpf((0, 4728779608739021, -52, 53)))*Min(mpf((0, 1, 0, 1)), _Dummy_4514 + mpf((0, 3602879701896397, -56, 52))) + (_Dummy_4514 - Dummy_4515 + mpf((0, 4728779608739021, -52, 53)))*Min(mpf((0, 1, 0, 1)), Max(_Dummy_4514 + mpf((0, 3602879701896397, -56, 52)), _Dummy_4514 - Dummy_4515 + mpf((0, 4728779608739021, -52, 53)))) + mpf((0, 1, -1, 1))*Min(mpf((0, 1, 0, 1)), _Dummy_4514 + mpf((0, 3602879701896397, -56, 52)))**2 - mpf((0, 1, -1, 1))*Min(mpf((0, 1, 0, 1)), Max(_Dummy_4514 + mpf((0, 3602879701896397, -56, 52)), _Dummy_4514 - Dummy_4515 + mpf((0, 4728779608739021, -52, 53))))**2)) + mpf((0, 5404319552844595, -53, 53))*((-(_Dummy_4514 - Dummy_4515 + mpf((0, 4278419646001971, -52, 52)))*Min(mpf((0, 1, 0, 1)), _Dummy_4514) + (_Dummy_4514 - Dummy_4515 + mpf((0, 4278419646001971, -52, 52)))*Min(mpf((0, 1, 0, 1)), Max(_Dummy_4514, _Dummy_4514 - Dummy_4515 + mpf((0, 4278419646001971, -52, 52)))) + mpf((0, 1, -1, 1))*Min(mpf((0, 1, 0, 1)), _Dummy_4514)**2 - mpf((0, 1, -1, 1))*Min(mpf((0, 1, 0, 1)), Max(_Dummy_4514, _Dummy_4514 - Dummy_4515 + mpf((0, 4278419646001971, -52, 52))))**2) if (Dummy_4515 >= mpf((0, 4278419646001971, -52, 52))) else (-(_Dummy_4514 - Dummy_4515 + mpf((0, 4278419646001971, -52, 52)))*Min(mpf((0, 1, 0, 1)), _Dummy_4514) + (_Dummy_4514 - Dummy_4515 + mpf((0, 4278419646001971, -52, 52)))*Min(mpf((0, 1, 0, 1)), Max(_Dummy_4514, _Dummy_4514 - Dummy_4515 + mpf((0, 1, 0, 1)))) + mpf((0, 1, -1, 1))*Min(mpf((0, 1, 0, 1)), _Dummy_4514)**2 - mpf((0, 1, -1, 1))*Min(mpf((0, 1, 0, 1)), Max(_Dummy_4514, _Dummy_4514 - Dummy_4515 + mpf((0, 1, 0, 1))))**2) if (Dummy_4515 >= 1) else ((mpf((0, 4278419646001971, -52, 52)) - Dummy_4515)*Min(mpf((0, 1, 0, 1)), _Dummy_4514) - (_Dummy_4514 - Dummy_4515 + mpf((0, 4278419646001971, -52, 52)))*Min(mpf((0, 1, 0, 1)), _Dummy_4514) + (_Dummy_4514 - Dummy_4515 + mpf((0, 4278419646001971, -52, 52)))*Min(mpf((0, 1, 0, 1)), Max(_Dummy_4514, _Dummy_4514 - Dummy_4515 + mpf((0, 4278419646001971, -52, 52)))) + mpf((0, 1, -1, 1))*Min(mpf((0, 1, 0, 1)), _Dummy_4514)**2 - mpf((0, 1, -1, 1))*Min(mpf((0, 1, 0, 1)), Max(_Dummy_4514, _Dummy_4514 - Dummy_4515 + mpf((0, 4278419646001971, -52, 52))))**2))], [(_Dummy_4514 + mpf((1, 3602879701896397, -55, 52)))*(mpf((0, 5404319552844595, -53, 53))*((-(-_Dummy_4514 + Dummy_4515 + mpf((0, 4728779608739021, -52, 53)))*Heaviside(1 - _Dummy_4514)*Heaviside(1 - Max(Dummy_4515 + mpf((0, 3602879701896397, -56, 52)), -_Dummy_4514 + Dummy_4515 + mpf((0, 4728779608739021, -52, 53)))) + Heaviside(1 - _Dummy_4514)*Heaviside(1 - Max(Dummy_4515 + mpf((0, 3602879701896397, -56, 52)), -_Dummy_4514 + Dummy_4515 + mpf((0, 4728779608739021, -52, 53))))*Min(mpf((0, 1, 0, 1)), Max(Dummy_4515 + mpf((0, 3602879701896397, -56, 52)), -_Dummy_4514 + Dummy_4515 + mpf((0, 4728779608739021, -52, 53)))) + Min(mpf((0, 1, 0, 1)), Dummy_4515 + mpf((0, 3602879701896397, -56, 52))) - Min(mpf((0, 1, 0, 1)), Max(Dummy_4515 + mpf((0, 3602879701896397, -56, 52)), -_Dummy_4514 + Dummy_4515 + mpf((0, 4728779608739021, -52, 53))))) if (_Dummy_4514 >= 1) else (-(-_Dummy_4514 + Dummy_4515 + mpf((0, 4728779608739021, -52, 53)))*Heaviside(1 - _Dummy_4514)*Heaviside(1 - Max(Dummy_4515 + mpf((0, 3602879701896397, -56, 52)), -_Dummy_4514 + Dummy_4515 + mpf((0, 4728779608739021, -52, 53)))) + Heaviside(1 - _Dummy_4514)*Heaviside(1 - Max(Dummy_4515 + mpf((0, 3602879701896397, -56, 52)), -_Dummy_4514 + Dummy_4515 + mpf((0, 4728779608739021, -52, 53))))*Min(mpf((0, 1, 0, 1)), Max(Dummy_4515 + mpf((0, 3602879701896397, -56, 52)), -_Dummy_4514 + Dummy_4515 + mpf((0, 4728779608739021, -52, 53)))) - Min(mpf((0, 1, 0, 1)), Max(Dummy_4515 + mpf((0, 3602879701896397, -56, 52)), -_Dummy_4514 + Dummy_4515 + mpf((0, 4728779608739021, -52, 53)))))) + mpf((0, 3602879701896397, -53, 52))*((0) if (_Dummy_4514 >= mpf((0, 4278419646001971, -52, 52))) else (-(-_Dummy_4514 + Dummy_4515 + mpf((0, 4278419646001971, -52, 52)))*Heaviside(1 - Max(Dummy_4515, -_Dummy_4514 + Dummy_4515 + 1))*Heaviside(-_Dummy_4514 + Dummy_4515 - Max(Dummy_4515, -_Dummy_4514 + Dummy_4515 + mpf((0, 4278419646001971, -52, 52))) + 1) + Heaviside(1 - Max(Dummy_4515, -_Dummy_4514 + Dummy_4515 + 1))*Heaviside(-_Dummy_4514 + Dummy_4515 - Max(Dummy_4515, -_Dummy_4514 + Dummy_4515 + mpf((0, 4278419646001971, -52, 52))) + 1)*Min(mpf((0, 1, 0, 1)), Max(Dummy_4515, -_Dummy_4514 + Dummy_4515 + mpf((0, 1, 0, 1))))) if (_Dummy_4514 >= 1) else (0))) + mpf((0, 5404319552844595, -53, 53))*((-(-_Dummy_4514 + Dummy_4515 + mpf((0, 4728779608739021, -52, 53)))*Min(mpf((0, 1, 0, 1)), Dummy_4515 + mpf((0, 3602879701896397, -56, 52))) + (-_Dummy_4514 + Dummy_4515 + mpf((0, 4728779608739021, -52, 53)))*Min(mpf((0, 1, 0, 1)), Max(Dummy_4515 + mpf((0, 3602879701896397, -56, 52)), -_Dummy_4514 + Dummy_4515 + mpf((0, 4728779608739021, -52, 53)))) + mpf((0, 1, -1, 1))*Min(mpf((0, 1, 0, 1)), Dummy_4515 + mpf((0, 3602879701896397, -56, 52)))**2 - mpf((0, 1, -1, 1))*Min(mpf((0, 1, 0, 1)), Max(Dummy_4515 + mpf((0, 3602879701896397, -56, 52)), -_Dummy_4514 + Dummy_4515 + mpf((0, 4728779608739021, -52, 53))))**2) if (_Dummy_4514 >= 1) else ((mpf((0, 1, 0, 1)) - _Dummy_4514)*Min(mpf((0, 1, 0, 1)), Dummy_4515 + mpf((0, 3602879701896397, -56, 52))) - (-_Dummy_4514 + Dummy_4515 + mpf((0, 4728779608739021, -52, 53)))*Min(mpf((0, 1, 0, 1)), Dummy_4515 + mpf((0, 3602879701896397, -56, 52))) + (-_Dummy_4514 + Dummy_4515 + mpf((0, 4728779608739021, -52, 53)))*Min(mpf((0, 1, 0, 1)), Max(Dummy_4515 + mpf((0, 3602879701896397, -56, 52)), -_Dummy_4514 + Dummy_4515 + mpf((0, 4728779608739021, -52, 53)))) + mpf((0, 1, -1, 1))*Min(mpf((0, 1, 0, 1)), Dummy_4515 + mpf((0, 3602879701896397, -56, 52)))**2 - mpf((0, 1, -1, 1))*Min(mpf((0, 1, 0, 1)), Max(Dummy_4515 + mpf((0, 3602879701896397, -56, 52)), -_Dummy_4514 + Dummy_4515 + mpf((0, 4728779608739021, -52, 53))))**2)) + mpf((0, 3602879701896397, -53, 52))*((-(-_Dummy_4514 + Dummy_4515 + mpf((0, 4278419646001971, -52, 52)))*Min(mpf((0, 1, 0, 1)), Dummy_4515) + (-_Dummy_4514 + Dummy_4515 + mpf((0, 4278419646001971, -52, 52)))*Min(mpf((0, 1, 0, 1)), Max(Dummy_4515, -_Dummy_4514 + Dummy_4515 + mpf((0, 4278419646001971, -52, 52)))) + mpf((0, 1, -1, 1))*Min(mpf((0, 1, 0, 1)), Dummy_4515)**2 - mpf((0, 1, -1, 1))*Min(mpf((0, 1, 0, 1)), Max(Dummy_4515, -_Dummy_4514 + Dummy_4515 + mpf((0, 4278419646001971, -52, 52))))**2) if (_Dummy_4514 >= mpf((0, 4278419646001971, -52, 52))) else (-(-_Dummy_4514 + Dummy_4515 + mpf((0, 4278419646001971, -52, 52)))*Min(mpf((0, 1, 0, 1)), Dummy_4515) + (-_Dummy_4514 + Dummy_4515 + mpf((0, 4278419646001971, -52, 52)))*Min(mpf((0, 1, 0, 1)), Max(Dummy_4515, -_Dummy_4514 + Dummy_4515 + mpf((0, 1, 0, 1)))) + mpf((0, 1, -1, 1))*Min(mpf((0, 1, 0, 1)), Dummy_4515)**2 - mpf((0, 1, -1, 1))*Min(mpf((0, 1, 0, 1)), Max(Dummy_4515, -_Dummy_4514 + Dummy_4515 + mpf((0, 1, 0, 1))))**2) if (_Dummy_4514 >= 1) else ((mpf((0, 4278419646001971, -52, 52)) - _Dummy_4514)*Min(mpf((0, 1, 0, 1)), Dummy_4515) - (-_Dummy_4514 + Dummy_4515 + mpf((0, 4278419646001971, -52, 52)))*Min(mpf((0, 1, 0, 1)), Dummy_4515) + (-_Dummy_4514 + Dummy_4515 + mpf((0, 4278419646001971, -52, 52)))*Min(mpf((0, 1, 0, 1)), Max(Dummy_4515, -_Dummy_4514 + Dummy_4515 + mpf((0, 4278419646001971, -52, 52)))) + mpf((0, 1, -1, 1))*Min(mpf((0, 1, 0, 1)), Dummy_4515)**2 - mpf((0, 1, -1, 1))*Min(mpf((0, 1, 0, 1)), Max(Dummy_4515, -_Dummy_4514 + Dummy_4515 + mpf((0, 4278419646001971, -52, 52))))**2))]]))

NameError: name 'Heaviside' is not defined

我根据对 (Some function) is not defined with SymPy Lambdify 的回答添加了 focx_lambda = lambdify((Pxe,Pye), focxeq, modules=['numpy', 'sympy']) 但这并没有改变任何东西。

我自己定义 Heaviside 的方式是

def Heaviside(x):
    if x<0:
        out=0
    else:
        out=1
    return out

我还尝试了 from numpy import * 以防万一。这并没有改变任何东西。

Answer 1

lambdify 的几个问题似乎同时发生。我想我可以让事情正常进行，但你应该检查它是否有意义，因为我不熟悉具体的方程式。

一般来说，调用 from sympy import * 和 from numpy import * 会造成很多困惑。两个库中的许多函数都有相同的名称，而且它们真的不喜欢与其他变量一起工作。

另一方面，lambdify 不能很好地与 Heaviside 配合使用。此外，numpy 中的函数是小写的 and 需要两个参数:一个 x 值和一个 x2 来决定 应该发生什么x==0。作为补救措施，下面的代码将“Heaviside”替换为 lambda x: np.heaviside(x, 1)。

我无法让 sympy 的 nsolve 使用这些函数，所以我尝试了 scipy 的 fsolve。 fsolve 还需要一些技巧才能处理函数元组。

创建focx_lambda 时，重要的是除了函数参数Pxe 和Pye 之外的所有变量都接收一个固定值。因此，我在执行 lambdify 时替换了它们。

from sympy import symbols, integrate, Min, Max, diff, lambdify
from IPython.display import display

mux, s, Px, Py, Pxe, Pye = symbols("mu_X s P_X P_Y P_X^* P_Y^*", positive=True)
vx, vy, cx, cy = symbols("v_X v_Y c_X c_Y", real=True)
pix = (Px - cx) * (mux * integrate(integrate(1, (vx, Min(1, Max(0, Px + Max(0, vy - Pye - s))), 1)), (vy, 0, 1))
                   + (1 - mux) * integrate(integrate(1, (vx, Min(1, Max(0, Max(Pxe + s, Px) + Max(0, vy - Pye))), 1)),
                                           (vy, 0, 1))
                   )
piy = (Py - cy) * ((1 - mux) * integrate(integrate(1, (vy, Min(1, Max(0, Py + Max(0, vx - Pxe - s))), 1)), (vx, 0, 1))
                   + mux * integrate(integrate(1, (vy, Min(1, Max(0, Max(Pye + s, Py) + Max(0, vx - Pxe))), 1)),
                                     (vx, 0, 1))
                   )
focx = diff(pix, Px)
focy = diff(piy, Py)
focxeq = focx.subs(Px, Pxe)
focyeq = focy.subs(Py, Pye)

import numpy as np
from scipy.optimize import fsolve

modules = [{'Heaviside': lambda x: np.heaviside(x, 1)}, 'numpy']
values_for_parameters = {mux: 0.4, s: 0.05, cx: 0, cy: 0.1}
focx_lambda = lambdify((Pxe, Pye), focxeq.subs(values_for_parameters), modules=modules)
focy_lambda = lambdify((Pxe, Pye), focyeq.subs(values_for_parameters), modules=modules)

print(focx_lambda(0.3, 0.4))  # we need to check that the lambdify works, so this should print a floating point number
print(focy_lambda(0.3, 0.4))

def equations(p):
    x, y = p
    return focx_lambda(x, y), focy_lambda(x, y)

sol = fsolve(equations, (0.3, 0.4))
print(sol)  # [0.64701372 0.61726372]