Scipy-非线性方程组求解

 原文地址:http://old.sebug.net/paper/books/scipydoc/scipy_intro.html,转载请注明出处!

ptimize库中的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倍。

雅可比矩阵

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

\begin{bmatrix} \dfrac{\partial f1}{\partial u1} & \dfrac{\partial f1}{\partial u2} & \dfrac{\partial f1}{\partial u3} \\[9pt] \dfrac{\partial f2}{\partial u1} & \dfrac{\partial f2}{\partial u2} & \dfrac{\partial f2}{\partial u3} \\[9pt] \dfrac{\partial f3}{\partial u1} & \dfrac{\partial f3}{\partial u2} & \dfrac{\partial f3}{\partial u3} \\ \end{bmatrix}

使用雅可比矩阵的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]
    ]

 

 

 

 

 

 

 

  • 0
    点赞
  • 12
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值