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

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

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

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

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 例子

多元多次方程组例如下面这种,三元二次方程组:

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


第二个方程组实在比较复杂,因此需要借助python。

二、python求解工具包

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

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

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

三、scipy方法

3.1 使用scipy的fsolve求解

关于scipy:下面这篇博文给的非常详细,https://blog.csdn.net/pipisorry/article/details/51106570

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

用fsolve相对初级,也相对简单易操作,代码较为简单,只用将方程的表达式写出运行即可。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,
    ]

 求解函数三个公式都为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 无法求解

部分较难情况无法求解

#!/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

4.1 二元一次方程组

较为简单,

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可以很好的进行求解

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

输出结果:

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

4.3 复数解

复数解也可以很好解出:

    # 复数解
    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 非线性求解

比如三角函数:

 程序均能较好解出

    # 非线性解
    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 较为复杂的二元二次方程

 此题较难,无论人来算,很难算出,用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))]

 复杂的问题终于解出,有四组实数解!

五、全部代码

#!/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确诊人数2月1日更新

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

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

概率相关实际问题汇总及解析

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

c++回溯法编程汇总

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

深度学习概览及主流模型演进

 

 

Logo

旨在为数千万中国开发者提供一个无缝且高效的云端环境,以支持学习、使用和贡献开源项目。

更多推荐