Scipy-root

创建时间: 2019-08-16 14:04:20
摘要: scipy求解方程，主要是非线性方程及非线性方程组

scipy.optimize.root

scipy求根使用的是root函数，结合源码root函数的相关说明。

声明

函数声明

def root(fun, x0, args=(), method='hybr', jac=None, tol=None, callback=None, options=None):

参数

fun : 可调用的

目标函数

x0 : numpy数组

估计的初值。

args : 元组，可选的

递给目标函数及其Jacobian矩阵的额外参数。

method : 字符串，可选的

定义求解器的类型，应该从以下列表选取：

method	参考
`hybr`	`optimize.root-hybr`
`lm`	`optimize.root-lm`
`broyden1`	`optimize.root-broyden1`
`broyden2`	`optimize.root-broyden2`
`anderson`	`optimize.root-anderson`
`linearmixing`	`optimize.root-linearmixing`
`diagbroyden`	`optimize.root-diagbroyden`
`excitingmixing`	`optimize.root-excitingmixing`
`krylov`	`optimize.root-krylov`
`df-sane`	`optimize.root-dfsane`
`

jac : 逻辑值或可调用的对象，可选的参考：

如果jac是逻辑值True，那么就假定fun返参考：矩阵。如果为False，那么Jacobian会通过数值方法估算。 jac也可以是可调用的对象，这种情况下，就像fun一样可以接受参数。

tol : 浮点数，可选的

迭代终止的容差。有关详细控制，使用求解器指定的选项options。

callback : 函数，可选的

可选的回调函数，以callback(x,f)的形式在每轮迭代调用，其中x是当前解，f是相应的函数值(residual?)。除了hybr及lm方法，都能使用回调函数。

options : 字典，可选的

求解器选项字典。如xtol或者maxiter。参考show_options()，可以通过show_options()给出可用的选项。

返回值

sol : OptimizeResult对象

结果以OptimizeResult对象的形式表示。一些重要的属性：x解数组，success表示算法成功退出与否的逻辑变量，message描述终止原因。更详细的信息请参考OptimizeResult。

注意

这一节描述可以通过method参数选择的求解器。默认的方法是hybr[1]。 hybr通过修改在MINPACK中实现的的Powell混合方法。 lm通过修改MINPACK中实现的Levenberg-Marquardt算法，以最小二乘法求解非线性方程组[1]。 df-sane是一种不使用导数的谱方法[3]。 broyden1,broyden2,anderson,linearmixing,diagbroyden,excitingmixing,krylov是使用回溯或全线性搜索的不精确的Newton方法。每个方法都对应一个特定的Jacobian近似。详细信息参考nonlin[2]。

Method broyden1 uses Broyden’s first Jacobian approximation, it is known as Broyden’s good method. 使用一阶Jacobian近似。
Method broyden2 uses Broyden’s second Jacobian approximation, it is known as Broyden’s bad method. 使用二阶Jacobian近似。
Method anderson uses (extended) Anderson mixing. 使用(扩展)Anderson混合。
Method Krylov uses Krylov approximation for inverse Jacobian. It is suitable for large-scale problem. 使用Krylov近似的Jacobian逆矩阵。适用于大规模问题。
Method diagbroyden uses diagonal Broyden Jacobian approximation.
Method linearmixing uses a scalar Jacobian approximation.
Method excitingmixing uses a tuned diagonal Jacobian approximation.

警告： The algorithms implemented for methods diagbroyden, linearmixing and excitingmixing may be useful for specific problems, but whether they will work may depend strongly on the problem.

versionadded:: 0.11.0

参考文献

示例

求解如下非线性方程组

\[\begin{split} \begin{align} x_0 +0.5 (x_0-x_1)^3-1=0 \\ 0.5(x_1 - x_0)^3 +x_1 = 0 \end{align} \end{split}\]

可视化可以发现求根即求解两条曲线的交点。求解及可视化代码如下：

# 求解
def fun(x):
    return [x[0]  + 0.5 * (x[0] - x[1])**3 - 1.0,
            0.5 * (x[1] - x[0])**3 + x[1]]
def jac(x):
    return np.array([[1 + 1.5 * (x[0] - x[1])**2,
                     -1.5 * (x[0] - x[1])**2],
                    [-1.5 * (x[1] - x[0])**2,
                     1 + 1.5 * (x[1] - x[0])**2]])
from scipy import optimize
sol = optimize.root(fun, [0, 0], jac=False, method='hybr')
sol

# 结果
        fjac: array([[-0.89914291,  0.43765515],
           [-0.43765515, -0.89914291]])
         fun: array([-1.11022302e-16,  0.00000000e+00])
     message: 'The solution converged.'
        nfev: 12
         qtf: array([ 1.19565972e-11, -4.12770392e-12])
           r: array([-2.16690469,  1.03701789, -1.10605417])
      status: 1
     success: True
           x: array([0.8411639, 0.1588361])

# 可视化
import matplotlib.pyplot as plt
def fun1(x):
    f = lambda x0: x0  + 0.5 * (x0 - x)**3 - 1.0
    s = optimize.root(f,0)
    return s.x
def fun2(x):
    f = lambda x0:0.5 * (x - x0)**3 + x
    s = optimize.root(f,0)
    return s.x
x1 = np.linspace(-10,10,100)
yi1 = np.array([fun1(i) for i in x1]) 
yi2 = np.array([fun2(i) for i in x1]) 
plt.plot(x1,yi1,'-',x1,yi2,'-',sol.x[1],sol.x[0],'ro')
plt.xlabel('$x_1$')
plt.ylabel('$x_0$')

plot

scipy.optimize.OptimizeResult

优化或者求根的结果是以OptimizeResult类给出。结合源码给出OptimizeResult相关说明。

声明

# 声明
class OptimizeResult(dict)

属性

x : ndarray

优化的解。

success : bool

优化器是否成功退出。

status : int

优化器的终止状态。其值取决于底层求解器。详细信息参考message。

message : str

描述终止原因。

fun, jac, hess: ndarray

目标函数、其Jacobian矩阵、其Hessian矩阵(如果存在的话)的值。 Hessian矩阵可能是近似值，参考问题中的函数文档。

hess_inv : object

目标函数的Hessian矩阵的逆；可能是近似值。并不是对所有求解器都可用，其类型可能是np.ndarray也可能是scipy.sparse.linalg.LinearOperator。

nfev, njev, nhev : int

目标函数及其Jacobian和Hessian矩阵的评价数（number of evaluations）。

nit : int

优化器迭代次数。

maxcv : float

最大约束冲突。

注意

对于特定的求解器，可能有额外的性质没有被列出。由于该类继承自dict，可以使用keys()方法查看哪些属性是可用的。

# 查看可用属性
sol.keys()

示例

非线性方程：单变量超越方程

\[x+2\cos (x)=0\]

# 求解
import numpy as np
from scipy.optimize import root

def func(x):
    return x+2*np.cos(x)

sol = root(func,1) # 默认使用`hybr`算法，及method='hybr'
sol

# 结果

            fjac: array([[-1.]])
            fun: array([0.])
        message: 'The solution converged.'
            nfev: 11
            qtf: array([3.32778249e-12])
            r: array([-2.71446049])
        status: 1
        success: True
            x: array([-1.02986653])

# 可视化

%matplotlib inline
import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.subplots()

xi = np.linspace(-5,5,1000)
yi = func(xi)
ax.plot(xi,yi)
ax.plot(xi,np.zeros(len(xi)),'--')
ax.plot(sol.x,func(sol.x),'ro')

plot

非线性方程组

\[\begin{split} \begin{aligned} x_0\cos (x_1) =4 \\ x_0x_1-x_1 =5 \end{aligned} \end{split}\]

这里自定义目标函数，返回函数及其Jacobian式，同时设置jac参数为True。使用Levenberg-Marquardt求解器。

# 求解
def func2(x):
    f=[x[0]*np.cos(x[1])-4,
       x[0]*x[1]-x[1]-5]
    df=np.array([[np.cos(x[1]),-x[0]*np.sin(x[1])],
                 [x[1],x[0]-1]])
    return f,df

sol = root(func2,[1,1],jac=True,method='lm') # func2返回函数及其jacbian矩阵，也可以分别给定
sol

# 结果
       cov_x: array([[ 0.87470958, -0.02852752],
           [-0.02852752,  0.01859874]])
        fjac: array([[ 7.52318843, -0.73161761],
           [ 0.24535902, -1.06922242]])
         fun: array([0., 0.])
        ipvt: array([2, 1], dtype=int32)
     message: 'The relative error between two consecutive iterates is at most 0.000000'
        nfev: 8
        njev: 7
         qtf: array([9.53474074e-13, 1.20388645e-13])
      status: 2
     success: True
           x: array([6.50409711, 0.90841421])

对比不使用jac加速，使用数值估算jac：

# 求解
def func3(x):
    f=[x[0]*np.cos(x[1])-4,
       x[0]*x[1]-x[1]-5]
    return f
sol = root(func3,[1,1],method='lm')
sol

# 结果
       cov_x: array([[ 0.87470957, -0.02852752],
           [-0.02852752,  0.01859874]])
        fjac: array([[ 7.52318845, -0.7316176 ],
           [ 0.24535901, -1.06922243]])
         fun: array([8.8817842e-16, 0.0000000e+00])
        ipvt: array([2, 1], dtype=int32)
     message: 'The relative error between two consecutive iterates is at most 0.000000'
        nfev: 22
         qtf: array([9.56065790e-13, 1.27078649e-13])
      status: 2
     success: True
           x: array([6.50409711, 0.90841421])

相比之下，不直接给出Jacbian表达式而是用数值估算，nfev会更大一些。