非线性方程组求解

Posted 2021-08-26 陨星落云

非线性方程组求解

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) # x0,x1,x2的值
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 倍。

雅可比矩阵
雅可比矩阵是一阶偏导数以一定方式排列的矩阵，它给出了可微分方程与给定点的最优线 性逼近，因此类似于多元函数的导数。例如前面的函数fl、f2、f3 和未知数u1、u2、u3 的雅可比矩阵如下：

​ 下面使用雅可比矩阵对方程组进行求解。计算雅可比矩阵的函数j()和f()— 样，其 x 参数是未知数的一组值，它计算非线性方程组在x 处的雅可比矩阵。通过fprime () 参数将 j()传递 给 fsolve() 。由于本例中的未知数很少，因此计算雅可比矩阵并不能显著地提高计算速度。

def f(x):
    x0,x1,x2 = x.tolist()
    return [5*x1+3,4*x0*x0-2*sin(x1*x2),x1*x2-1.5]

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]

摘自《python科学计算》