3.2.1 非线性方程组求解
fsolve( )可以对非线性方程组进行求解，它的基本调用形式为fsolve(func, x0)。其中func是计算方程组误差的函数，它的参数x是一个数组，其值为方程组的一组可能的解。func返回将x代入方程组之后得到的每个方程的误差，x0为未知数的一组初始值。假设要对下面的方程组进行求解：
f1(u1,u2,u3)=0, f2(u1,u2,u3)=0, f3(u1,u2,u3)=0
那么func函数可以定义如下：
def func(x):
u1,u2,u3 = x
return [f1(u1,u2,u3), f2(u1,u2,u3), f3(u1,u2,u3)]
下面我们看一个对下列方程组求解的例子：
from math import sin, cos
from scipy import optimize
def f(x): ❶
x0, x1, x2 = x.tolist( ) ❷
return [
5*
x1+3,
4*
x0*
x0 - 2*
sin(x1*
x2),
x1*
x2 - 1.5
]
# f计算方程组的误差，[1,1,1]是未知数的初始值
result = optimize.fsolve(f, [1,1,1]) ❸
print result
print f(result)
[-0.70622057 -0.6 -2.5 ]
[0.0, -9.126033262418787e-14, 5.329070518200751e-15]
❶f( )是计算方程组的误差的函数，x参数是一组可能的解。fsolve( )在调用f( )时，传递给f( )的参数是一个数组。❷先调用数组的tolist( )方法，将其转换为Python的标准浮点数列表，然后调用math模块中的函数进行运算。因为在进行单个数值的运算时，标准浮点类型比NumPy的浮点类型要快许多，所以把数值都转换成标准浮点数类型，能缩短一些计算时间。❸调用fsolve( )时，传递计算误差的函数f( )以及未知数的初始值。
在对方程组进行求解时，fsolve( )会自动计算方程组在某点对各个未知数变量的偏导数，这些偏导数组成一个二维数组，数学上称之为雅可比矩阵。如果方程组中的未知数很多，而与每个方程有关联的未知数较少，即雅可比矩阵比较稀疏时，将计算雅可比矩阵的函数作为参数传递给fsolve( )，这能大幅度提高运算速度。笔者在一个模拟计算的程序中需要求解有50个未知数的非线性方程组。每个方程平均与6个未知数相关，通过传递计算雅可比矩阵的函数使fsolve( )的计算速度提高了4倍。
雅可比矩阵
雅可比矩阵是一阶偏导数以一定方式排列的矩阵，它给出了可微分方程与给定点的最优线性逼近，因此类似于多元函数的导数。例如前面的函数f1、f2、f3和未知数u1、u2、u3的雅可比矩阵如下：
下面使用雅可比矩阵对方程组进行求解。❶计算雅可比矩阵的函数j( )和f( )一样，其x参数是未知数的一组值，它计算非线性方程组在x处的雅可比矩阵。❷通过fprime参数将j( )传递给fsolve( )。由于本例中的未知数很少，因此计算雅可比矩阵并不能显著地提高计算速度。
def j(x): ❶
x0, x1, x2 = x.tolist( )
return [
[0, 5, 0],
[8*
x0, -2*
x2*
cos(x1*
x2), -2*
x1*
cos(x1*
x2)],
[0, x2, x1]
]
result = optimize.fsolve(f, [1,1,1], fprime=j) ❷
print result
print f(result)
[-0.70622057 -0.6 -2.5 ]
[0.0, -9.126033262418787e-14, 5.329070518200751e-15]