优化非线性方程组 - 处理复数
Optimizing System of Non Linear Equations - Handling Complex Numbers
我正在尝试求解非线性方程组。问题是解决方案很复杂,根据 Octave/Matlab,虚部非常小。我正在尝试将其移至 python,但不幸的是我不确定我应该如何优雅地处理它。
在Octave中,我可以直接使用fsolve,然后将解通过"real"函数传递,得到数字的实部。问题是,它很容易解决它而不会返回任何错误
不幸的是,在尝试求解方程式时,在 python returns 错误中使用 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)
在 Octave 中实现同样的事情给出了这个解决方案:
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
我通过 Octave 中的 "real" 函数传递此解决方案,以获得我想要的结果。
特别是python连一次方程式都解不出来。特别是如果我尝试 运行 f[4] 在定义了所有参数的函数之外,它 returns 一个 nan 值。
如有任何帮助,我们将不胜感激!
提前向我道歉 missed/formatted-badly。
的确,scipy 很难处理复数。但是,一个名为 mpmath 的项目可以解决您的问题。此处:http://mpmath.org/. It used to come with sympy (sympy.org). You can find documentation 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)
我正在尝试求解非线性方程组。问题是解决方案很复杂,根据 Octave/Matlab,虚部非常小。我正在尝试将其移至 python,但不幸的是我不确定我应该如何优雅地处理它。
在Octave中,我可以直接使用fsolve,然后将解通过"real"函数传递,得到数字的实部。问题是,它很容易解决它而不会返回任何错误
不幸的是,在尝试求解方程式时,在 python returns 错误中使用 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)
在 Octave 中实现同样的事情给出了这个解决方案:
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
我通过 Octave 中的 "real" 函数传递此解决方案,以获得我想要的结果。
特别是python连一次方程式都解不出来。特别是如果我尝试 运行 f[4] 在定义了所有参数的函数之外,它 returns 一个 nan 值。
如有任何帮助,我们将不胜感激!
提前向我道歉 missed/formatted-badly。
的确,scipy 很难处理复数。但是,一个名为 mpmath 的项目可以解决您的问题。此处:http://mpmath.org/. It used to come with sympy (sympy.org). You can find documentation 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)