求解线性系统答案

【问题标题】：Solving linear system求解线性系统
【发布时间】：2017-12-13 09:51:56
【问题描述】：

我有 16 个方程：

a, b, c, d, e, f, g = symbols('a, b, c, d, e, f, g', real=True)
z00, z01, z02, z03 = symbols('z00, z01, z02, z03', real=True)
z10, z11, z12, z13 = symbols('z10, z11, z12, z13', real=True)
z20, z21, z22, z23 = symbols('z20, z21, z22, z23', real=True)
z30, z31, z32, z33 = symbols('z30, z31, z32, z33', real=True)    

Eq(z00, g)
Eq(z03, b + d + f + g)
Eq(-z00 + 6*z01 - 12*z02 + 8*z03, 8*b + 4*d + 2*f + g)
Eq(-8*z00 + 36*z01 - 54*z02 + 27*z03, 27*b + 9*d + 3*f + g)
Eq(z30, a + c + e + g)
Eq(z33, a + b + c + d + e + f + g)
Eq(-z30 + 6*z31 - 12*z32 + 8*z33, a + 8*b + c + 4*d + e + 2*f + g)
Eq(-8*z30 + 36*z31 - 54*z32 + 27*z33, a + 27*b + c + 9*d + e + 3*f + g)
Eq(-z00 + 6*z10 - 12*z20 + 8*z30, 8*a + 4*c + 2*e + g)
Eq(-z03 + 6*z13 - 12*z23 + 8*z33, 8*a + b + 4*c + d + 2*e + f + g)
Eq(z00 - 6*z01 + 12*z02 - 8*z03 - 6*z10 + 36*z11 - 72*z12 + 48*z13 + 12*z20 - 72*z21 + 144*z22 - 96*z23 - 8*z30 + 48*z31 - 96*z32 + 64*z33, 8*a + 8*b + 4*c + 4*d + 2*e + 2*f + g)
Eq(8*z00 - 36*z01 + 54*z02 - 27*z03 - 48*z10 + 216*z11 - 324*z12 + 162*z13 + 96*z20 - 432*z21 + 648*z22 - 324*z23 - 64*z30 + 288*z31 - 432*z32 + 216*z33, 8*a + 27*b + 4*c + 9*d + 2*e + 3*f + g)
Eq(-8*z00 + 36*z10 - 54*z20 + 27*z30, 27*a + 9*c + 3*e + g)
Eq(-8*z03 + 36*z13 - 54*z23 + 27*z33, 27*a + b + 9*c + d + 3*e + f + g)
Eq(8*z00 - 48*z01 + 96*z02 - 64*z03 - 36*z10 + 216*z11 - 432*z12 + 288*z13 + 54*z20 - 324*z21 + 648*z22 - 432*z23 - 27*z30 + 162*z31 - 324*z32 + 216*z33, 27*a + 8*b + 9*c + 4*d + 3*e + 2*f + g)
Eq(64*z00 - 288*z01 + 432*z02 - 216*z03 - 288*z10 + 1296*z11 - 1944*z12 + 972*z13 + 432*z20 - 1944*z21 + 2916*z22 - 1458*z23 - 216*z30 + 972*z31 - 1458*z32 + 729*z33, 27*a + 27*b + 9*c + 9*d + 3*e + 3*f + g)

并希望以a、b、c、...的形式表达（解决）所有z..。

所以我所做的就是将方程的左侧和右侧转换为矩阵，如下所示：

AA = Matrix([[z00], [z03], [-z00 + 6*z01 - 12*z02 + 8*z03], [-8*z00 + 36*z01 - 54*z02 + 27*z03], [z30], [z33], [-z30 + 6*z31 - 12*z32 + 8*z33], [-8*z30 + 36*z31 - 54*z32 + 27*z33], [-z00 + 6*z10 - 12*z20 + 8*z30], [-z03 + 6*z13 - 12*z23 + 8*z33], [z00 - 6*z01 + 12*z02 - 8*z03 - 6*z10 + 36*z11 - 72*z12 + 48*z13 + 12*z20 - 72*z21 + 144*z22 - 96*z23 - 8*z30 + 48*z31 - 96*z32 + 64*z33], [8*z00 - 36*z01 + 54*z02 - 27*z03 - 48*z10 + 216*z11 - 324*z12 + 162*z13 + 96*z20 - 432*z21 + 648*z22 - 324*z23 - 64*z30 + 288*z31 - 432*z32 + 216*z33], [-8*z00 + 36*z10 - 54*z20 + 27*z30], [-8*z03 + 36*z13 - 54*z23 + 27*z33], [8*z00 - 48*z01 + 96*z02 - 64*z03 - 36*z10 + 216*z11 - 432*z12 + 288*z13 + 54*z20 - 324*z21 + 648*z22 - 432*z23 - 27*z30 + 162*z31 - 324*z32 + 216*z33], [64*z00 - 288*z01 + 432*z02 - 216*z03 - 288*z10 + 1296*z11 - 1944*z12 + 972*z13 + 432*z20 - 1944*z21 + 2916*z22 - 1458*z23 - 216*z30 + 972*z31 - 1458*z32 + 729*z33]])
bb = Matrix([[g], [b + d + f + g], [8*b + 4*d + 2*f + g], [27*b + 9*d + 3*f + g], [a + c + e + g], [a + b + c + d + e + f + g], [a + 8*b + c + 4*d + e + 2*f + g], [a + 27*b + c + 9*d + e + 3*f + g], [8*a + 4*c + 2*e + g], [8*a + b + 4*c + d + 2*e + f + g], [8*a + 8*b + 4*c + 4*d + 2*e + 2*f + g], [8*a + 27*b + 4*c + 9*d + 2*e + 3*f + g], [27*a + 9*c + 3*e + g], [27*a + b + 9*c + d + 3*e + f + g], [27*a + 8*b + 9*c + 4*d + 3*e + 2*f + g], [27*a + 27*b + 9*c + 9*d + 3*e + 3*f + g]])

然后我尝试了：

from sympy.solvers.solveset import linsolve
linsolve((AA, bb), [z00])
linsolve((AA, bb), [z00, z01, z02, z03, z10, z11, z12, z13, z20, z21, z22, z23, z30, z31, z32, z33])

但我得到了EmptySet()。我知道至少z00 应该等于g 所以我期待至少某种非空答案。

请问，有谁知道如何求解这些方程，即为什么我会得到 EmptySet()？

【问题讨论】：

标签： python sympy

【解决方案1】：

试试solve()

_Z = symbols('z0(0:4),z1(0:4),z2(0:4),z3(0:4)', real=True)
solve(AA - bb, _Z)

不知道为什么linsolve() 像下面这样返回EmptySet，但我猜这与表达式中符号的顺序有关。

linsolve(AA - bb, _Z)

但下面的代码有效

A = Matrix([[ 1,    0,   0,    0,    0,    0,     0,   0,   0,     0,    0,     0,    0,   0,     0,   0],
        [ 0,    0,   0,    1,    0,    0,     0,   0,   0,     0,    0,     0,    0,   0,     0,   0],
        [-1,    6, -12,    8,    0,    0,     0,   0,   0,     0,    0,     0,    0,   0,     0,   0],
        [-8,   36, -54,   27,    0,    0,     0,   0,   0,     0,    0,     0,    0,   0,     0,   0],
        [ 0,    0,   0,    0,    0,    0,     0,   0,   0,     0,    0,     0,    1,   0,     0,   0],
        [ 0,    0,   0,    0,    0,    0,     0,   0,   0,     0,    0,     0,    0,   0,     0,   1],
        [ 0,    0,   0,    0,    0,    0,     0,   0,   0,     0,    0,     0,   -1,   6,   -12,   8],
        [ 0,    0,   0,    0,    0,    0,     0,   0,   0,     0,    0,     0,   -8,  36,   -54,  27],
        [-1,    0,   0,    0,    6,    0,     0,   0, -12,     0,    0,     0,    8,   0,     0,   0],
        [ 0,    0,   0,   -1,    0,    0,     0,   6,   0,     0,    0,   -12,    0,   0,     0,   8],
        [ 1,   -6,  12,   -8,   -6,   36,   -72,  48,  12,   -72,  144,   -96,   -8,  48,   -96,  64],
        [ 8,  -36,  54,  -27,  -48,  216,  -324, 162,  96,  -432,  648,  -324,  -64, 288,  -432, 216],
        [-8,    0,   0,    0,   36,    0,     0,   0, -54,     0,    0,     0,   27,   0,     0,   0],
        [ 0,    0,   0,   -8,    0,    0,     0,  36,   0,     0,    0,   -54,    0,   0,     0,  27],
        [ 8,  -48,  96,  -64,  -36,  216,  -432, 288,  54,  -324,  648,  -432,  -27, 162,  -324, 216],
        [64, -288, 432, -216, -288, 1296, -1944, 972, 432, -1944, 2916, -1458, -216, 972, -1458, 729]])
assert A * Matrix(_Z) == AA
linsolve((A, bb))

【讨论】：

非常感谢！最后一个linsolve((A, bb)) 应该是linsolve((A, bb), _Z)，我相信...
顺便说一句，你是怎么把Eq(64*z00 - 288*z01 + 432*z02...变成[64, -288, 432...的？ SymPy 中是否有执行此类转换的命令？
我是手动完成的 :) 起初我认为矩阵可能有些奇怪（单数？）。现在我想起来了，我可以通过循环遍历符号来手动提取系数。