python求解多元多次方程组或非线性方程组

Myth丶恋晨 2023-07-08 02:35 10阅读 0赞

**背景：**如何使用python求解多元多次方程组或者非线性方程组。

原创内容，转载注明出处！请勿用于商业用途！

（上篇用python拟合2019nCov感染人数的文章被不少博主转载了，发的比较早，不少博主在文章基础上添加新内容也新发了新的更新后的预测，或者加入一些新的模块。博文链接如下：）

[python实现logistic增长模型拟合2019-nCov确诊人数 ][python_logistic_2019-nCov_]

[python实现logistic增长模型拟合2019-nCov确诊人数2月1日更新 ][python_logistic_2019-nCov_2_1_]

[博客文章总目录-邢翔瑞的技术博客 ][-_]

**目录**

一、多元多次方程

1.1 定义

1.2 例子

二、python求解工具包

三、scipy方法

3.1 使用scipy的fsolve求解

3.2 非完备解

3.3 非线性方程的解

3.4 无法求解

四、sympy工具包求解

4.1 二元一次方程组

4.2 多解

4.3 复数解

4.4 非线性求解

4.5 较为复杂的二元二次方程

五、全部代码

--------------------

# 一、多元多次方程 #

## 1.1 定义 ##

我们常见的方程组有一元一次方程组，比如x+3=5这种，很简单很好解。

*  二元一次方程组，即方程组中有两个未知数，未知数的最高次数为1.
 *  二元二次方程组：方程组中有两个未知数，未知数的最高次数为2.。此类方程组均有公式解法或者成形的解法。

但是面临多元多次方程组，解法错综复杂，是数学家们研究的内容。为了更好的解决此类问题，我们可以用python来实现。

## 1.2 例子 ##

多元多次方程组例如下面这种，三元二次方程组：

![20200224143343595.png][]

下面这种，二元二次方程组。

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl8zNjQ3NDgwOQ_size_16_color_FFFFFF_t_70][]  
第二个方程组实在比较复杂，因此需要借助python。

# 二、python求解工具包 #

python求解方程组的工具包较多。例如：

*  numpy：numpy.linalg.solve 可以直接求解线性方程组，numpy是python非常常用的包，解的方程也较为初级。
 *  scipy：from scipy.optimize import fsolve，可以求解非线性方程组，使用较为方便，但是解集并不完备，可能漏掉一下解（后文会给个例子）scipy可以用于数学、科学、工程领域的常用软件包,可以处理插值、积分、优化，相对较初级易用
 *  sympy：此工具包功能相对强大，支持符号计算、高精度计算、解方程、微积分、组合数学、离散数学、几何学、概率与统计、物理学等方面的功能。github地址：[https://github.com/sympy/sympy][https_github.com_sympy_sympy]
 *  sage，不支持位运算，z3约束求解器，等其他工具包，本文不详述，感兴趣的可以查找相应的内容。

本文详细讲述scipy以及sympy求解多次方程的方法。

# 三、scipy方法 #

## 3.1 使用scipy的fsolve求解 ##

关于scipy：下面这篇博文给的非常详细，[https://blog.csdn.net/pipisorry/article/details/51106570][https_blog.csdn.net_pipisorry_article_details_51106570]

我们只将求解方程的部分。

用fsolve相对初级，也相对简单易操作，代码较为简单，只用将方程的表达式写出运行即可。fsolve近似看作用最小二乘法求解。不够很强大，很多情况下解集不完备或者无法解出。

例如对于![20200224143343595.png][]，首先要定义相应的函数：

def solve_function(unsolved_value):
        x,y,z=unsolved_value[0],unsolved_value[1],unsolved_value[2]
        return [
            x**2+y**2-10,
            y**2+z**2-34,
            x**2+z**2-26,
        ]

求解函数三个公式都为0时候的解，中括号内为初值\[0, 0, 0\]

solved=fsolve(solve_function,[0, 0, 0])

全部代码：

#!/usr/bin/python
    # -*- coding: UTF-8 -*-
    """
    python解方程
    """
    
    from scipy.optimize import fsolve
    
    def solve_function(unsolved_value):
        x,y,z=unsolved_value[0],unsolved_value[1],unsolved_value[2]
        return [
            x**2+y**2-10,
            y**2+z**2-34,
            x**2+z**2-26,
        ]
    
    solved=fsolve(solve_function,[0, 0, 0])
    print(solved)
    
    
    print("Program done!")
    
    """
    运行结果：
    [-1.  3.  5.]
    Program done!
    """

看出运行结果来看，此结果并非完备解集。因为x，y，z都是可正可负。例如1或者-1，3或者-3，5或者-5，但是此工具包只能解出一个解。

## 3.2 非完备解 ##

显而易见，x\*\*2-9=0的解为3或者-3

def solve_function(unsolved_value):
        x=unsolved_value[0]
        return [
            x**2-9,
        ]
    
    solved=fsolve(solve_function,[0])

但是程序只能得出一个结果3，但是得不到-3

## 3.3 非线性方程的解 ##

最简单的sin(x)=0.5,则x可能为π/6或者 5π/6

#!/usr/bin/python
    # -*- coding: UTF-8 -*-
    """
    python解方程
    """
    
    from scipy.optimize import fsolve
    from math import sin,cos
    
    def solve_function(unsolved_value):
        x=unsolved_value[0]
        return [
            sin(x)-0.5
        ]
    
    solved=fsolve(solve_function,[3.14])
    print(solved)
    solved=fsolve(solve_function,[0])
    print(solved)
    
    print("Program done!")

运行结果为：

[2.61799388]
    [0.52359878]
    Program done!

可以解出π/6或者 5π/6，中括号内为初始迭代的值。

## 3.4 无法求解 ##

部分较难情况无法求解

![20200224165317385.png][]

#!/usr/bin/python
    # -*- coding: UTF-8 -*-
    """
    python解方程
    """
    
    from scipy.optimize import fsolve
    
    def solve_function(unsolved_value):
        x,y=unsolved_value[0],unsolved_value[1]
        return [
            x*x+2*x*y,
            2*x*y-2*y*y
        ]
    
    solved=fsolve(solve_function,[6, -3])
    print(solved)
    
    print("Program done!")

无法求解会给出报错，和用最小二乘法迭代得到明显错误的解。

[1.64526700e-115 1.33665018e-115]
    A:\python\python\lib\site-packages\scipy\optimize\minpack.py:162: RuntimeWarning: The number of calls to function has reached maxfev = 600.
    Program done!
      warnings.warn(msg, RuntimeWarning)

# 四、sympy工具包求解 #

没安装可以在teiminal中pip install sympy，此工具包涉及支持符号计算、高精度计算、模式匹配、绘图、解方程、微积分、组合数学、离散 数学、几何学、概率与统计、物理学等方面的功能。功能较为强大，解方程组时性能也较好。

官方地址：[https://github.com/sympy/sympy][https_github.com_sympy_sympy]

## 4.1 二元一次方程组 ##

![20200224171509557.png][]

较为简单，

from sympy import *
    # 二元一次方程
    x = Symbol('x')
    y = Symbol('y')
    solved_value=solve([2*x+y-1, x-2*y], [x, y])
    print(solved_value)

此方法较为简单，但是相应的自变量应当写成符号的形式，x=Symbol('x')

求解后有分数解：

{x: 2/5, y: 1/5}
    Program done!

## 4.2 多解 ##

多解情况与复数解

例如，多个解的情况，sympy可以很好的进行求解

![20200224173101122.png][]

x = Symbol('x')
        solved_value=solve([x**2-9], [x])
        print(solved_value)

输出结果：

[(-3,), (3,)]

## 4.3 复数解 ##

复数解也可以很好解出：

![20200224174006527.png][]

# 复数解
        solved_value = solve([x ** 2 + 9], [x])
        print(solved_value)
        solved_value = solve([x ** 4 - 9], [x])
        print(solved_value)
    """
    运行结果：
    [(-3*I,), (3*I,)]
    [(-sqrt(3),), (sqrt(3),), (-sqrt(3)*I,), (sqrt(3)*I,)]
    """

复数解也能较好解出

## 4.4 非线性求解 ##

比如三角函数：

![20200224174602295.png][]

程序均能较好解出

# 非线性解
        solved_value = solve([sin(x) - 0.5], [x])
        print(solved_value)
        solved_value = solve([sin(x) - 1], [x])
        print(solved_value)
    
    """
    [(0.523598775598299,), (2.61799387799149,)]
    [(pi/2,)]
    """

## 4.5 较为复杂的二元二次方程 ##

![20200224165317385.png][]

此题较难，无论人来算，很难算出，用scipy工具包也迭代不出解。但是sympy强大的功能可以很好的解出此方程。

# 二元二次方程组
        x = Symbol('x')
        y=  Symbol('y')
        solved_value=solve([x**2+2*x*y-6,2*x*y-2*y**2+3], [x,y])
        print(solved_value)

有四组实数解：

[(-(-3 + sqrt(13))*sqrt(sqrt(13)/2 + 2), -sqrt(sqrt(13)/2 + 2)),
     ((-3 + sqrt(13))*sqrt(sqrt(13)/2 + 2), sqrt(sqrt(13)/2 + 2)), 
    (-sqrt(2 - sqrt(13)/2)*(-sqrt(13) - 3), -sqrt(2 - sqrt(13)/2)), 
    (sqrt(2 - sqrt(13)/2)*(-sqrt(13) - 3), sqrt(2 - sqrt(13)/2))]

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl8zNjQ3NDgwOQ_size_16_color_FFFFFF_t_70 1][]

复杂的问题终于解出，有四组实数解！

# 五、全部代码 #

#!/usr/bin/python
    # -*- coding: UTF-8 -*-
    
    """
    python解方程
    created by xingxinagrui on 2020.2.24
    """
    
    from scipy.optimize import fsolve
    from math import sin,cos
    from sympy import *
    
    # 1-4 scipy
    # 5-7 sympy
    part=7
    
    if part==1:
        # 求解非线性方程组
        def solve_function(unsolved_value):
            x=unsolved_value[0]
            return [
                sin(x)-0.5
            ]
    
        solved=fsolve(solve_function,[3.14])
        print(solved)
        solved=fsolve(solve_function,[0])
        print(solved)
    
    if part==2:
        # 求解三元二次方程组
        def solve_function(unsolved_value):
            x, y, z = unsolved_value[0], unsolved_value[1], unsolved_value[2]
            return [
                x ** 2 + y ** 2 - 10,
                y ** 2 + z ** 2 - 34,
                x ** 2 + z ** 2 - 26,
            ]
    
        solved = fsolve(solve_function, [0, 0, 0])
        print(solved)
    
    
    if part==3:
        #解的非完备性
        def solve_function(unsolved_value):
            x = unsolved_value[0]
            return [
                x ** 2 - 9,
            ]
    
        solved = fsolve(solve_function, [0])
        print(solved)
    
    
    if part == 4:
        # 较难无法求解
        def solve_function(unsolved_value):
            x, y = unsolved_value[0], unsolved_value[1]
            return [
                x * x + 2 * x * y,
                2 * x * y - 2 * y * y
            ]
    
        solved = fsolve(solve_function, [6, -3])
        print(solved)
    
    if part == 5:
        # 二元一次方程
        x = Symbol('x')
        y = Symbol('y')
        solved_value=solve([2*x+y-1, x-2*y], [x, y])
        print(solved_value)
    
    
    if part == 6:
        # 多解情况
        x = Symbol('x')
        solved_value=solve([x**2-9], [x])
        print(solved_value)
    
        # 复数解
        solved_value = solve([x ** 2 + 9], [x])
        print(solved_value)
        solved_value = solve([x ** 4 - 9], [x])
        print(solved_value)
    
    
        # 非线性解
        solved_value = solve([sin(x) - 0.5], [x])
        print(solved_value)
        solved_value = solve([sin(x) - 1], [x])
        print(solved_value)
    
    if part == 7:
        # 二元二次方程组
        x = Symbol('x')
        y=  Symbol('y')
        solved_value=solve([x**2+2*x*y-6,2*x*y-2*y**2+3], [x,y])
        print(solved_value)
    
    print("Program done!")

博主其他文章：

[博客文章总目录-邢翔瑞的技术博客 ][-_]

[python实现logistic增长模型拟合2019-nCov确诊人数 ][python_logistic_2019-nCov_]

[python实现logistic增长模型拟合2019-nCov确诊人数2月1日更新 ][python_logistic_2019-nCov_2_1_]

[支持向量机(Support Vector Machine,SVM)算法复杂度详解 ][Support Vector Machine_SVM_]

[王者荣耀中的数学原理及游戏策略(一)防御篇(护甲|魔抗|伤害运算机制) ][Link 1]

[概率相关实际问题汇总及解析 ][Link 2]

[机器学习算法基础问题(三)集成学习|adaboost与XGboost| EM算法 ][adaboost_XGboost_ EM_]

[c++回溯法编程汇总 ][c_]

[c++动态规划类算法编程汇总(四)集合的子集|最长子序列(矩阵)的和(积) | 最大子矩阵 ][c_ _ _]

[深度学习概览及主流模型演进 ][Link 3]

[python_logistic_2019-nCov_]: https://blog.csdn.net/weixin_36474809/article/details/104101055
[python_logistic_2019-nCov_2_1_]: https://blog.csdn.net/weixin_36474809/article/details/104119494
[-_]: https://blog.csdn.net/weixin_36474809/article/details/88884137
[20200224143343595.png]: https://img-blog.csdnimg.cn/20200224143343595.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl8zNjQ3NDgwOQ_size_16_color_FFFFFF_t_70]: https://img-blog.csdnimg.cn/20200224143534871.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl8zNjQ3NDgwOQ==,size_16,color_FFFFFF,t_70
[https_github.com_sympy_sympy]: https://github.com/sympy/sympy
[https_blog.csdn.net_pipisorry_article_details_51106570]: https://blog.csdn.net/pipisorry/article/details/51106570
[20200224165317385.png]: https://img-blog.csdnimg.cn/20200224165317385.png
[20200224171509557.png]: https://img-blog.csdnimg.cn/20200224171509557.png
[20200224173101122.png]: https://img-blog.csdnimg.cn/20200224173101122.png
[20200224174006527.png]: https://img-blog.csdnimg.cn/20200224174006527.png
[20200224174602295.png]: https://img-blog.csdnimg.cn/20200224174602295.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl8zNjQ3NDgwOQ_size_16_color_FFFFFF_t_70 1]: https://img-blog.csdnimg.cn/20200224181217867.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl8zNjQ3NDgwOQ==,size_16,color_FFFFFF,t_70
[Support Vector Machine_SVM_]: https://blog.csdn.net/weixin_36474809/article/details/103986845
[Link 1]: https://blog.csdn.net/weixin_36474809/article/details/102087639
[Link 2]: https://blog.csdn.net/weixin_36474809/article/details/101448598
[adaboost_XGboost_ EM_]: https://blog.csdn.net/weixin_36474809/article/details/101274046
[c_]: https://blog.csdn.net/weixin_36474809/article/details/100713385
[c_ _ _]: https://blog.csdn.net/weixin_36474809/article/details/100581124
[Link 3]: https://blog.csdn.net/weixin_36474809/article/details/100534384