from scipy.optimize import fsolve, root
import numpy as np
# 定义方程内容
def f(x, *arg):
return arg[0] * 2 ** (1 - x) / (1 - x) + (1 - arg[0]) * 1.6 ** (1 - x) / (1 - x) - arg[0] * 3.85 ** (1 - x) / (
1 - x) - (1 - arg[0]) * 0.1 ** (1 - x) / (1 - x)
# 参数p为一个超参数
results = [[p * 0.1, fsolve(f, x0=0, args=(p * 0.1))[0]] for p in range(2, 10)]
print(np.array(results), '\n')
results = [[p * 0.1, root(f, x0=0, args=(p * 0.1))['x'][0]] for p in range(2, 10)]
print(np.array(results))
# output results
[[ 0.2 -0.94683705]
[ 0.3 -0.48657472]
[ 0.4 -0.14263228]
[ 0.5 0.14636333]
[ 0.6 0.4114561 ]
[ 0.7 0.67618002]
[ 0.8 0.9705809 ]
[ 0.9 1.3683912 ]]
[[ 0.2 -0.94683705]
[ 0.3 -0.48657472]
[ 0.4 -0.14263228]
[ 0.5 0.14636333]
[ 0.6 0.4114561 ]
[ 0.7 0.67618002]
[ 0.8 0.9705809 ]
[ 0.9 1.3683912 ]]
对此非线性方程组，根据方程图像（如下图）可知，应有4组不同的解。以下尝试用不同的方法求出该方程组的解。
"""利用gekko求解非线性方程组"""
from gekko import GEKKO
m = GEKKO()
x = m.Var(value=0) # 给定初值为0
y = m.Var(value=0) # 给定初值为0
m.Equations([x ** 2 / 4 + y ** 2 == 1,
(x - 0.2) ** 2 - y == 3])
m.solve(disp=False)
x, y = x.value, y.value
print(x, y)
输出结果：
[-1.2961338938] [-0.7615833719]
换不同初值计算得到的结果如下：
[-1.6818042485] [0.54118722964] # 给定初值为(-2,0)
[1.9760411678] [0.15432222765] # 给定初值为(2,0)
[1.8018969861] [-0.43392604594] # 给定初值为(2,-2)
[1.9760412095] [0.15432236862] # 给定初值为(10,10)
[1.801896954] [-0.4339261545] # 给定初值为(10,-10)
可知，用这种方法并不能得到方程组的全部解，并且最终得到的解为其解集中与给定的初值“距离”较近的一个。
optimize库中的fsolve函数可以用来对非线性方程组进行求解。
from scipy.optimize import fsolve
def f(X):
x = X[0]
y = X[1]
return [x ** 2 / 4 + y ** 2 - 1,
(x - 0.2) ** 2 - y - 3]
X0 = [0, 0]
result = fsolve(f, X0)
print(result)
输出结果：
[-1.29613389 -0.76158337]
换不同初值计算得到的结果如下：
[-1.68180425 0.54118723] # 给定初值为(-2,0)
[1.97604116 0.15432219] # 给定初值为(2,0)
[ 1.80189699 -0.43392605] # 给定初值为(2,-2)
[1.97604116 0.15432219] # 给定初值为(10,10)
[ 1.80189699 -0.43392605] # 给定初值为(10,-10)
可知，用这种方法也不能得到方程组的全部解，并且最终得到的解与给定初值有关。
from scipy.optimize import fsolve, root
def f(X):
x = X[0]
y = X[1]
return [x ** 2 / 4 + y ** 2 - 1,
(x - 0.2) ** 2 - y - 3]
X0 = [10, 10]
result1 = fsolve(f, X0)
result2 = root(f, X0)
print(result2)
输出结果：
fjac: array([[-0.2547064 , -0.96701843],
[ 0.96701843, -0.2547064 ]])
fun: array([-3.34943184e-12, 2.75734990e-12])
message: 'The solution converged.'
nfev: 22
qtf: array([-1.65320424e-10, -2.73193431e-10])
r: array([-3.70991104, 0.8956477 , 0.56891317])
status: 1
success: True
x: array([1.97604116, 0.15432219])
结果与fsolve函数得到的结果相同。
from scipy.optimize import leastsq
def f(X):
x = X[0]
y = X[1]
return [x ** 2 / 4 + y ** 2 - 1,
(x - 0.2) ** 2 - y - 3]
X0 = [10, 10]
h = leastsq(f, X0)
print(h)
输出结果：
(array([1.97604116, 0.15432219]), 2)
5.1 利用solve求解所有精确解
from sympy import symbols, Eq, solve, nsolve
x, y = symbols('x y')
eqs = [Eq(x ** 2 / 4 + y ** 2, 1),
Eq((x - 0.2) ** 2 - y, 3)]
print(solve(eqs, [x, y]))
输出结果：
[[-1.68180424847377 + 1.56760579250585e-32*I
0.541187229573922 - 3.01196919624356e-31*I]
[-1.29613389377477 + 1.95607066863502e-32*I
-0.761583371898353 + 3.93313832308616e-31*I]
[1.80189698634479 - 1.95607066863926e-32*I
-0.433926045139482 - 8.10475677027422e-31*I]
[1.97604115590375 - 1.56760579250161e-32*I
0.154322187463913 + 7.18358764343162e-31*I]]
可以看出，用这种方法能够得到方程组的全部解，并且为精确解，缺点是求解时间较长。
5.2 利用nsolve求解数值解
from sympy import symbols, Eq, nsolve
x, y = symbols('x y')
eqs = [Eq(x ** 2 / 4 + y ** 2, 1),
Eq((x - 0.2) ** 2 - y, 3)]
X0 = [3, 4]
print(nsolve(eqs, [x, y], X0))
输出结果：
Matrix([[1.97604115590375], [0.154322187463913]])
nsolve为数值求解，需要指定一个初始值，初始值会影响最终得到哪一个解（如果有多解的话），而且初始值设的不好，则可能找不到解。
scipy.optimize.root求解速度快，但只能得到靠近初始值的一个解。对形式简单、有求根公式的方程，sympy.solve能够得到所有严格解，但当方程组变量较多时，它求起来会很慢。而且对于不存在求根公式的复杂方程，sympy.solve无法求解。
以上为个人经验，希望能给大家一个参考，也希望大家多多支持。
