>>> 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]