0. 问题实例

$$\{\begin{array}{r}10x-y-2z=72\\ -x+10y-2z=83\\ -x-y+5z=42\end{array}$$

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]])