处理复数的非线性方程组优化系统

Question

我想解一个非线性方程组。问题是，根据Octave/Matlab，解是复杂的，只有很小的虚部。我正试图将此移到python上，但不幸的是，我不确定应该如何优雅地处理它。你知道吗

在倍频程中，我可以直接使用fsolve，然后通过“real”函数传递解，得到数字的实数部分。问题是，它很容易解决它而不返回任何错误

不幸的是，在python中使用numpy在尝试求解方程时会返回错误。下面是用Python编写的公式：

import numpy as np
from scipy.optimize import fsolve
import scipy.io as spio

params = dict()
params['cbeta'] = 0.96
params['cdelta'] = 0.1
params['calpha'] = 0.33
params['cgamma'] = 1.2
params['clambda']= 1.0
params['csigma'] = 0.8
params['etau'] = 0.0

def steady_s(vars0):
    # unpacking paramters
    cbeta = params['cbeta']
    cdelta = params['cdelta']
    calpha = params['calpha']
    cgamma = params['cgamma']
    clambda= params['clambda']
    csigma = params['csigma']

    # guesses for initial values
    c = vars0[0]
    y = vars0[1]
    k = vars0[2]
    g = vars0[3]
    r = vars0[4]

    # == functions to minimize to find steady states == #
    f = np.empty((5,))

    # HH Euler
    f[0] = (1.0/c)*cbeta*(r + 1.0 - cdelta) - (1.0+g)/c

    # Goods market clearing
    f[1] = y - c - k*(1.0 + g) + k*(1.0-cdelta)

    # Capital Market clearing
    f[2] = r - (k)**(calpha-1.0)*calpha**2.0

    # production function for final good
    f[3] = y - k**calpha

    # growth rate
    pi = (calpha - 1.0) * k**calpha #small pi, this isnt actual profits
    f[4] = g - (cgamma - 1.0) * clambda * (csigma*clambda*pi)**(csigma/(1.0-csigma))

    return f

# == Initial Guesses == #
vars0 = np.ones((5,))

# == Solving for Steady State == #

xss = fsolve(steady_s, vars0)

在八度音阶中实现同样的功能可以得到以下解决方案：

 Columns 1 through 3:

   0.7851388 + 0.0000000i   0.8520544 + 0.0000000i   0.6155938 + 0.0000000i

 Columns 4 and 5:

   0.0087008 - 0.0000000i   0.1507300 - 0.0000000i

我把这个解通过“实”函数，得到我想要的结果。你知道吗

特别是，python甚至很难一次解决这些方程。特别是如果我尝试在函数外部运行f[4]并定义所有参数，它将返回一个nan值。你知道吗

任何帮助都将不胜感激！你知道吗

对我遗漏的/格式错误的任何内容提前道歉。你知道吗

Answer 1

实际上，scipy与复数作斗争。然而，一个名为mpmath的项目可以解决您的问题。这里：http://mpmath.org/。它过去总是伴随着同情(sympy.org网站). 您可以找到文档here：这个解决方案适合我：

from mpmath import findroot
import numpy as np
import scipy.io as spio

params = dict()
params['cbeta'] = 0.96
params['cdelta'] = 0.1
params['calpha'] = 0.33
params['cgamma'] = 1.2
params['clambda']= 1.0
params['csigma'] = 0.8
params['etau'] = 0.0

def steady_s(c,y,k,g,r):
    # unpacking paramters
    cbeta = params['cbeta']
    cdelta = params['cdelta']
    calpha = params['calpha']
    cgamma = params['cgamma']
    clambda= params['clambda']
    csigma = params['csigma']

    # guesses for initial values
    #c = vars0[0]
    #y = vars0[1]
    #k = vars0[2]
    #g = vars0[3]
    #r = vars0[4]

    # == functions to minimize to find steady states == #
    f = [0,0,0,0,0]

    # HH Euler
    f[0] = (1.0/c)*cbeta*(r + 1.0 - cdelta) - (1.0+g)/c

    # Goods market clearing
    f[1] = y - c - k*(1.0 + g) + k*(1.0-cdelta)

    # Capital Market clearing
    f[2] = r - (k)**(calpha-1.0)*calpha**2.0

    # production function for final good
    f[3] = y - k**calpha

    # growth rate
    pi = (calpha - 1.0) * k**calpha #small pi, this isnt actual profits
    f[4] = g - (cgamma - 1.0) * clambda * (csigma*clambda*pi)**(csigma/(1.0-csigma))

    return f

# == Initial Guesses == #
vars0 = list(np.ones((5,)))

# == Solving for Steady State == #
xss = findroot(steady_s, vars0)