>>> from env_helper import info; info()
页面更新时间: 2023-06-24 12:31:10
运行环境:
    Linux发行版本: Debian GNU/Linux 12 (bookworm)
    操作系统内核: Linux-6.1.0-9-amd64-x86_64-with-glibc2.36
    Python版本: 3.11.2

2.3. 非线性方程组求解

optimize库中的 fsolve() 函数可以用来对非线性方程组进行求解。它的基本调用形式如下:

fsolve(func, x0)

func(x) 是计算方程组误差的函数,它的参数 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)]

下面是一个实际的例子,求解如下方程组的解:

• 5*x1 + 3 = 0
• 4*x0*x0 - 2*sin(x1*x2) = 0
• x1*x2 - 1.5 = 0

程序如下:

>>> from scipy.optimize import fsolve
>>> from math import sin,cos
>>> def f(x):
>>>     x0 = float(x[0])
>>>     x1 = float(x[1])
>>>     x2 = float(x[2])
>>>     return [
>>>         5*x1+3,
>>>         4*x0*x0 - 2*sin(x1*x2),
>>>         x1*x2 - 1.5
>>>     ]
>>>
>>> result = fsolve(f, [1,1,1])
>>>
>>> print(result)
>>> print(f(result))
[-0.70622057 -0.6        -2.5       ]
[0.0, -9.126033262418787e-14, 5.329070518200751e-15]

由于fsolve函数在调用函数f时,传递的参数为数组,因此如果直接使用数组中的元素计算的话,计算 速度将会有所降低,因此这里先用float函数将数组中的元素转换为Python中的标准浮点数,然后调用 标准math库中的函数进行运算。

在对方程组进行求解时,fsolve会自动计算方程组的雅可比矩阵,如果方程组中的未知数很多,而与每 个方程有关的未知数较少时,即雅可比矩阵比较稀疏时,传递一个计算雅可比矩阵的函数将能大幅度 提高运算速度。笔者在一个模拟计算的程序中需要大量求解近有50个未知数的非线性方程组的解。 每个方程平均与6个未知数相关,通过传递雅可比矩阵的计算函数使计算速度提高了4倍。

2.3.1. 雅可比矩阵

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

  ∂f 1 ∂f 1 ∂f 1  ∂u1 ∂u2 ∂u3     ∂f 2 ∂f 2 ∂f 2     ∂u1 ∂u2 ∂u3      ∂f 3 ∂f 3 ∂f 3 ∂u1 ∂u2 ∂u3

使用雅可比矩阵的fsolve实例如下,计算雅可比矩阵的函数j通过fprime参数传递给fsolve,函数j和函 数f一样,有一个未知数的解矢量参数x,函数j计算非线性方程组在矢量x点上的雅可比矩阵。由于这个 例子中未知数很少,因此程序计算雅可比矩阵并不能带来计算速度的提升。

>>> from scipy.optimize import fsolve
>>> from math import sin,cos
>>> def f(x):
>>>     x0 = float(x[0])
>>>     x1 = float(x[1])
>>>     x2 = float(x[2])
>>>     return [
>>>     5*x1+3,
>>>     4*x0*x0 - 2*sin(x1*x2),
>>>     x1*x2 - 1.5
>>> ]
>>>
>>>
>>> def j(x):
>>>     x0 = float(x[0])
>>>     x1 = float(x[1])
>>>     x2 = float(x[2])
>>>     return [
>>>         [0, 5, 0],
>>>         [8*x0, -2*x2*cos(x1*x2), -2*x1*cos(x1*x2)],
>>>         [0, x2, x1]
>>>     ]
>>> result = 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]