文章目录
0. 问题实例
{ 10 x − y − 2 z = 72 − x + 10 y − 2 z = 83 − x − y + 5 z = 42 \left \{ \begin{aligned} 10x-y-2z=72 \\ -x+10y-2z=83 \\ -x-y+5z=42 \end{aligned} \right. ⎩⎪⎨⎪⎧10x−y−2z=72−x+10y−2z=83−x−y+5z=42
1. 利用gekko的GEKKO求解
"""利用gekko求解线性方程组"""
from gekko import GEKKO
m = GEKKO() # 定义模型
x = m.Var() # 定义模型变量,初值为0
y = m.Var()
z = m.Var()
m.Equations([10 * x - y - 2 * z == 72,
-x + 10 * y - 2 * z == 83,
-x - y + 5 * z == 42, ]) # 方程组
m.solve(disp=False) # 求解
x, y, z = x.value, y.value, z.value
print(x,y,z) # 打印结果
输出结果:
[11.0] [12.0] [13.0]
2. 利用scipy的linalg求解
from scipy import linalg
import numpy as np
A = np.array([[10, -1, -2], [-1, 10, -2], [-1, -1, 5]]) # A代表系数矩阵
b = np.array([72, 83, 42]) # b代表常数列
x = linalg.solve(A, b)
print(x)
输出结果:
[11. 12. 13.]
3. 利用scipy.optimize的root或fsolve求解
from scipy.optimize import root, fsolve
def f(X):
x = X[0]
y = X[1]
z = X[2] # 切分变量
return [10 * x - y - 2 * z - 72,
-x + 10 * y - 2 * z - 83,
-x - y + 5 * z - 42]
X0 = [1, 2, 3] # 设定变量初值
m1 = root(f, X0).x # 利用root求解并给出结果
m2 = fsolve(f, X0) # 利用fsolve求解并给出结果
print(m1)
print(m2)
输出结果:
[11. 12. 13.]
[11. 12. 13.]
4. 利用Numpy的linalg求解
import numpy as np
A = np.array([[10, -1, -2], [-1, 10, -2], [-1, -1, 5]]) # A为系数矩阵
b = np.array([72, 83, 42]) # b为常数列
inv_A = np.linalg.inv(A) # A的逆矩阵
x = inv_A.dot(b) # A的逆矩阵与b做点积运算
x = np.linalg.solve(A, b) # 5,6两行也可以用本行替代
print(x)
输出结果:
[11. 12. 13.]
import numpy as np
# A = np.mat([[10, -1, -2], [-1, 10, -2], [-1, -1, 5]]) # A为系数矩阵
# b = np.mat([[72], [83], [42]]) # b为常数列
A = np.mat("10, -1, -2; -1, 10, -2; -1, -1, 5") # A为系数矩阵
b = np.mat("72;83;42") # b为常数列
inv_A = np.linalg.inv(A) # A的逆矩阵
inv_A = A.I # A的逆矩阵
# x = inv_A.dot(b) # A的逆矩阵与b做点积运算
x = np.linalg.solve(A, b)
print(x)
输出结果:
[11. 12. 13.]
5. 利用sympy的solve和nsolve求解
5.1 利用solve求解所有精确解
from sympy import symbols, Eq, solve
x, y, z = symbols('x y z')
eqs = [Eq(10 * x - y - 2 * z, 72),
Eq(-x + 10 * y - 2 * z, 83),
Eq(-x - y + 5 * z, 42)]
print(solve(eqs, [x, y, z]))
输出结果:
{x: 11, y: 12, z: 13}
5.1 利用nsolve求解数值解
from sympy import symbols, Eq, nsolve
x, y, z = symbols('x y z')
eqs = [Eq(10 * x - y - 2 * z, 72),
Eq(-x + 10 * y - 2 * z, 83),
Eq(-x - y + 5 * z, 42)]
initialValue = [1, 2, 3]
print(nsolve(eqs, [x, y, z], initialValue))
输出结果:
Matrix([[11.0000000000000], [12.0000000000000], [13.0000000000000]])