命令 linsolve 或 solve from sympy 无法求解一维线性方程
The commands linsolve or solve from sympy cannot solve a 1D linear equation
为了引入问题，我把我的代码写在下面：
import sympy as sp
sp.init_printing()
import numpy as np
import math
n=2
u,v=sp.symbols('u v')
a,b,c,d,h=sp.symbols('a b c d h')
k=sp.symbols('k')
diffmatrix=sp.zeros(n)
for i in range(n):
diffmatrix[i,i] = sp.symbols('D'+str(i+1))
globals()['D'+str(i+1)]=sp.symbols('D'+str(i+1))
f=[]
var=[]
var.append('u')
var.append('v')
# f.append(u+v)
# f.append(u-v**2)
f.append(a-c*u+d*v+u**2*v)
f.append(b+h*u-d*v-u**2*v)
f=sp.Matrix(f)
var=sp.Matrix(var)
jacobianmat=f.jacobian(var)
phi=[]
alpha=[]
beta=[]
psi=[]
gamma=[]
for i in range(n):
phi.append(sp.symbols('phi' + str(i)))
alpha.append(sp.symbols('alpha' + str(i)))
beta.append(sp.symbols('beta'+str(i)))
psi.append(sp.symbols('psi' + str(i)))
gamma.append(sp.symbols('gamma' + str(i)))
globals()['phi' + str(i)] = sp.symbols('phi' + str(i))
globals()['alpha' + str(i)] = sp.symbols('alpha' + str(i))
globals()['beta' + str(i)] = sp.symbols('beta' + str(i))
globals()['psi' + str(i)] = sp.symbols('psi' + str(i))
globals()['gamma' + str(i)] = sp.symbols('gamma' + str(i))
auxiliarterm=sp.linsolve(sp.Matrix([np.dot((jacobianmat-(k**2)*diffmatrix)[0:n-1],phi[0:n-1])])+sp.Matrix(jacobianmat-(k**2)*diffmatrix)[0:n-1,n-1],phi[0:n-1])
phi[0:n-1]=sp.Matrix(list(auxiliarterm))
phi[n-1]=1
doublesummationsecondorder=[]
for functionnumber in range(n):
tempsum=0
for i in range(n):
for j in range(n):
tempsum=sp.Add(tempsum,phi[i]*phi[j]*sp.diff(f[functionnumber],var[i],var[j]))
tempsum=-tempsum/4
doublesummationsecondorder.append(tempsum)
doublesummationsecondorder=sp.Matrix(doublesummationsecondorder)
alpha=sp.linsolve(sp.Matrix(np.dot(jacobianmat,alpha))-doublesummationsecondorder,alpha)
alpha=sp.Matrix(sp.Transpose(sp.Matrix(list(alpha))))
beta=sp.linsolve(sp.Matrix(np.dot(jacobianmat-4*k**2*diffmatrix,beta))-doublesummationsecondorder,beta)
beta=sp.Matrix(sp.Transpose(sp.Matrix(list(beta))))
doublesummationthirdorder1=[]
doublesummationthirdorder2=[]
triplesummationthirdorder=[]
for functionnumber in range(n):
tempsum1=0
tempsum2=0
tempsum3=0
for i in range(n):
for j in range(n):
tempsum1=sp.Add(tempsum1,phi[i]*(alpha[j]+beta[j]/2)*sp.diff(f[functionnumber],var[i],var[j]))
tempsum2=sp.Add(tempsum2,phi[i]*beta[j]*sp.diff(f[functionnumber],var[i],var[j]))
for l in range(n):
tempsum3=sp.Add(tempsum3,phi[i]*phi[j]*phi[l]*sp.diff(f[functionnumber],var[i],var[j],var[l]))
doublesummationthirdorder1.append(tempsum1)
doublesummationthirdorder2.append(tempsum2)
triplesummationthirdorder.append(tempsum3)
doublesummationthirdorder1=sp.Matrix(doublesummationthirdorder1)
doublesummationthirdorder2=sp.Matrix(doublesummationthirdorder2)
triplesummationthirdorder=sp.Matrix(triplesummationthirdorder)
auxiliarterm=sp.linsolve(sp.Matrix([np.dot((sp.matrices.Transpose(jacobianmat)-(k**2)*diffmatrix)[0:n-1],psi[0:n-1])])+sp.Matrix(sp.matrices.Transpose(jacobianmat)-(k**2)*diffmatrix)[0:n-1,n-1],psi[0:n-1])
psi[0:n-1]=sp.Matrix(list(auxiliarterm))
psi[n-1]=1
print('Almost ready')
auxiliarterm=sp.solve(sp.Matrix([np.dot((jacobianmat-(k**2)*diffmatrix)[0:n-1],gamma[0:n-1])])+sp.Matrix(doublesummationthirdorder1[0:n-1])+sp.Matrix(3/math.factorial(4)*sp.Matrix(triplesummationthirdorder[0:n-1])),gamma[0:n-1])
如您所见，当 n=2 时，最后求解的方程是一维（线性）方程。问题是，由于某种原因，我无法得到它的解决方案。在任何情况下，如果我尝试使用 f=[u+v,u-v^2] 给出的另一个向量场，那么一切正常。我想知道这是内存的问题，符号计算的效率还是我的代码有问题。
万事如意。
问题是求解器会尝试扩展您拥有的表达式，并且这些表达式在扩展时会变大，所以速度很慢。对于符号抵消的线性方程组，通常需要扩展。在这种情况下，因为它只是一个方程式，所以如果您确定存在唯一解，则可以跳过该步骤：
In [41]: a = Wild('a')
In [42]: b = Wild('b')
In [43]: match = eq[0].match(a*s[0]+b)
In [44]: sol = -match[b]/match[a]
In [45]: sol
Out[45]:
⎛ 2 2
⎛ 2⎞ ⎜ 2⋅D₁⋅k ⋅u + 2⋅c⋅u + d⋅v - 3⋅u ⋅v
2⋅u⋅⎝d + u ⎠⋅⎜───────────────────────────────────────────────────────────── + ──────────────────────────────────────────────────────────────────────────────
⎜ 2 4 2 2 2 2 2 ⎛ 3 8 3 6 3 6 2 2 6 2 6
⎝2⋅D₁ ⋅k + 4⋅D₁⋅c⋅k - 8⋅D₁⋅k ⋅u⋅v + 2⋅c - 8⋅c⋅u⋅v + 8⋅u ⋅v 2⋅⎝32⋅D₁ ⋅D₂⋅k + 8⋅D₁ ⋅d⋅k + 8⋅D₁ ⋅k ⋅u + 72⋅D₁ ⋅D₂⋅c⋅k - 144⋅D₁ ⋅D₂⋅k ⋅u⋅
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
2 4 2 4 3 2
8⋅D₁ ⋅d⋅k ⋅u + 8⋅D₁ ⋅k ⋅u + 10⋅D₁⋅c⋅d⋅k ⋅u +
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
2 4 2 4 2 2 4 2 4 2 4 2 2 4 3 2 4 4 4 2
v + 18⋅D₁ ⋅c⋅d⋅k + 18⋅D₁ ⋅c⋅k ⋅u - 2⋅D₁ ⋅d⋅h⋅k - 32⋅D₁ ⋅d⋅k ⋅u⋅v - 2⋅D₁ ⋅h⋅k ⋅u - 32⋅D₁ ⋅k ⋅u ⋅v + 48⋅D₁⋅D₂⋅c ⋅k - 192⋅D₁⋅D₂⋅c⋅k ⋅u⋅v + 192⋅D₁⋅D₂⋅k ⋅u
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
2 3 2 2 2 2 2 2 3 2 4 2 2 3 2 2
10⋅D₁⋅c⋅k ⋅u + 4⋅D₁⋅d ⋅k ⋅v - 2⋅D₁⋅d⋅h⋅k ⋅u - 8⋅D₁⋅d⋅k ⋅u ⋅v - 2⋅D₁⋅h⋅k ⋅u - 12⋅D₁⋅k ⋅u ⋅v + 2⋅c ⋅d⋅u + 2⋅c ⋅u + c⋅d ⋅v - 2⋅c⋅d⋅h⋅u - 2⋅c⋅d⋅u ⋅v - 2⋅c⋅h
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
2 2 2 2 2 2 2 2 2 2 2 3 2 2 2 2 2
⋅v + 12⋅D₁⋅c ⋅d⋅k + 12⋅D₁⋅c ⋅k ⋅u - 4⋅D₁⋅c⋅d⋅h⋅k - 40⋅D₁⋅c⋅d⋅k ⋅u⋅v - 4⋅D₁⋅c⋅h⋅k ⋅u - 40⋅D₁⋅c⋅k ⋅u ⋅v + 8⋅D₁⋅d⋅h⋅k ⋅u⋅v + 32⋅D₁⋅d⋅k ⋅u ⋅v + 8⋅D₁⋅h⋅k ⋅
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
2
D₁⋅k + c - 2⋅u⋅v
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
3 4 2 2 4
⋅u - 3⋅c⋅u ⋅v - d ⋅h⋅v + 2⋅d⋅h⋅u ⋅v + 3⋅h⋅u ⋅v
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
3 2 4 2 3 2 2 2 2 2 2 2 3 3 3 3 2 2 2 2 2 2 3
u ⋅v + 32⋅D₁⋅k ⋅u ⋅v + 8⋅D₂⋅c ⋅k - 48⋅D₂⋅c ⋅k ⋅u⋅v + 96⋅D₂⋅c⋅k ⋅u ⋅v - 64⋅D₂⋅k ⋅u ⋅v + 2⋅c ⋅d + 2⋅c ⋅u - 2⋅c ⋅d⋅h - 8⋅c ⋅d⋅u⋅v - 2⋅c ⋅h⋅u - 8⋅c ⋅u ⋅v
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
⎞
⎟
───────────────────────────────────────────────────────────────────────────────⎟
2 2 3 4 2 2 2 4 2⎞⎟
+ 8⋅c⋅d⋅h⋅u⋅v + 8⋅c⋅d⋅u ⋅v + 8⋅c⋅h⋅u ⋅v + 8⋅c⋅u ⋅v - 8⋅d⋅h⋅u ⋅v - 8⋅h⋅u ⋅v ⎠⎠
──────────────────────────────────────────────────────────────────────────────── - ─────────────────────────────────────────────────────────────────────────
3 8 3 6 3 6 2 2 6 2 6
16⋅D₁ ⋅D₂⋅k + 4⋅D₁ ⋅d⋅k + 4⋅D₁ ⋅k ⋅u + 36⋅D₁ ⋅D₂⋅c⋅k - 72⋅D₁ ⋅D₂⋅k ⋅u
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
2 4 2 4 2 2 4 2 4 2 4 2 2 4 3 2 4 4 4 2 2
⋅v + 9⋅D₁ ⋅c⋅d⋅k + 9⋅D₁ ⋅c⋅k ⋅u - D₁ ⋅d⋅h⋅k - 16⋅D₁ ⋅d⋅k ⋅u⋅v - D₁ ⋅h⋅k ⋅u - 16⋅D₁ ⋅k ⋅u ⋅v + 24⋅D₁⋅D₂⋅c ⋅k - 96⋅D₁⋅D₂⋅c⋅k ⋅u⋅v + 96⋅D₁⋅D₂⋅k ⋅u ⋅v + 6
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
⎛ 4 4 3 2 2 3 2 2 2 2 2 4 ⎞
u⋅⎝4⋅D₁⋅D₂⋅d⋅k ⋅u + 4⋅D₁⋅D₂⋅k ⋅u + 4⋅D₂⋅c⋅d⋅k ⋅u + 4⋅D₂⋅c⋅k ⋅u + 2⋅D₂⋅d ⋅k ⋅v - 4⋅D₂⋅d⋅k ⋅u ⋅v - 6⋅D₂⋅k ⋅u ⋅v⎠
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
2 2 2 2 2 2 2 2 2 2 3 2 2 2 2 2 3
⋅D₁⋅c ⋅d⋅k + 6⋅D₁⋅c ⋅k ⋅u - 2⋅D₁⋅c⋅d⋅h⋅k - 20⋅D₁⋅c⋅d⋅k ⋅u⋅v - 2⋅D₁⋅c⋅h⋅k ⋅u - 20⋅D₁⋅c⋅k ⋅u ⋅v + 4⋅D₁⋅d⋅h⋅k ⋅u⋅v + 16⋅D₁⋅d⋅k ⋅u ⋅v + 4⋅D₁⋅h⋅k ⋅u ⋅v + 16
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
2
- D₁⋅k - c + 2⋅u⋅v
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
2 4 2 3 2 2 2 2 2 2 2 3 3 3 3 2 2 2 2 2 2 3
⋅D₁⋅k ⋅u ⋅v + 4⋅D₂⋅c ⋅k - 24⋅D₂⋅c ⋅k ⋅u⋅v + 48⋅D₂⋅c⋅k ⋅u ⋅v - 32⋅D₂⋅k ⋅u ⋅v + c ⋅d + c ⋅u - c ⋅d⋅h - 4⋅c ⋅d⋅u⋅v - c ⋅h⋅u - 4⋅c ⋅u ⋅v + 4⋅c⋅d⋅h⋅u⋅v + 4
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
───────────────────────────────────────────────────────────── - ────────────────────────────────────────────────────────────────────────────────────────────
2 2 3 4 2 2 2 4 2 ⎛ 2 ⎞ ⎛ 3 8 3 6 3 6 2 2 6 2 6
⋅c⋅d⋅u ⋅v + 4⋅c⋅h⋅u ⋅v + 4⋅c⋅u ⋅v - 4⋅d⋅h⋅u ⋅v - 4⋅h⋅u ⋅v ⎝D₁⋅k + c - 2⋅u⋅v⎠⋅⎝16⋅D₁ ⋅D₂⋅k + 4⋅D₁ ⋅d⋅k + 4⋅D₁ ⋅k ⋅u + 36⋅D₁ ⋅D₂⋅c⋅k - 72⋅D₁ ⋅D₂⋅k
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
2 4 2 4 2 2 4 2 4 2 4 2 2 4 3 2 4 4 4 2 2
⋅u⋅v + 9⋅D₁ ⋅c⋅d⋅k + 9⋅D₁ ⋅c⋅k ⋅u - D₁ ⋅d⋅h⋅k - 16⋅D₁ ⋅d⋅k ⋅u⋅v - D₁ ⋅h⋅k ⋅u - 16⋅D₁ ⋅k ⋅u ⋅v + 24⋅D₁⋅D₂⋅c ⋅k - 96⋅D₁⋅D₂⋅c⋅k ⋅u⋅v + 96⋅D₁⋅D₂⋅k ⋅u ⋅v +
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
⎛ 2⎞ ⎛ 4 4 3 2 2 3 2 2 2 2 2 4 ⎞
v⋅⎝d + u ⎠⋅⎝4⋅D₁⋅D₂⋅d⋅k ⋅u + 4⋅D₁⋅D₂⋅k ⋅u + 4⋅D₂⋅c⋅d⋅k ⋅u + 4⋅D₂⋅c⋅k ⋅u + 2⋅D₂⋅d ⋅k ⋅v - 4⋅D₂⋅d⋅k ⋅u ⋅v - 6⋅D₂⋅k ⋅u ⋅v⎠
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
2 2 2 2 2 2 2 2 2 2 3 2 2 2 2 2 3
6⋅D₁⋅c ⋅d⋅k + 6⋅D₁⋅c ⋅k ⋅u - 2⋅D₁⋅c⋅d⋅h⋅k - 20⋅D₁⋅c⋅d⋅k ⋅u⋅v - 2⋅D₁⋅c⋅h⋅k ⋅u - 20⋅D₁⋅c⋅k ⋅u ⋅v + 4⋅D₁⋅d⋅h⋅k ⋅u⋅v + 16⋅D₁⋅d⋅k ⋅u ⋅v + 4⋅D₁⋅h⋅k ⋅u ⋅v +
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
2 4 2 3 2 2 2 2 2 2 2 3 3 3 3 2 2 2 2 2 2 3
16⋅D₁⋅k ⋅u ⋅v + 4⋅D₂⋅c ⋅k - 24⋅D₂⋅c ⋅k ⋅u⋅v + 48⋅D₂⋅c⋅k ⋅u ⋅v - 32⋅D₂⋅k ⋅u ⋅v + c ⋅d + c ⋅u - c ⋅d⋅h - 4⋅c ⋅d⋅u⋅v - c ⋅h⋅u - 4⋅c ⋅u ⋅v + 4⋅c⋅d⋅h⋅u⋅v +
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
2
⎛ 2⎞
0.75⋅⎝d + u ⎠
──────────────────────────────────────────────────────────────── - ────────────────────
2 2 3 4 2 2 2 4 2⎞ 2
4⋅c⋅d⋅u ⋅v + 4⋅c⋅h⋅u ⋅v + 4⋅c⋅u ⋅v - 4⋅d⋅h⋅u ⋅v - 4⋅h⋅u ⋅v ⎠ ⎛ 2 ⎞
⎝D₁⋅k + c - 2⋅u⋅v⎠
───────────────────────────────────────────────────────────────────────────────────────
另请注意，其中有一个 Float (0.75)。我会确保这里只使用有理数，因为浮点系数的多项式运算没有很好的定义。
正如 Oscar 指出的那样，表达式很大，求解过程在处理所有符号时陷入困境。以下例程将简化表达式，让您使用一小组符号：
def keepx(eq, x):
"""collapse non-x dependent Add and Mul into Duumy symbols
and return e, r where e is the new expresion and r is
the dictionary that can restore e: eq == e.xreplace(r)
Examples
========
>>> from sympy import solve, symbols
>>> from sympy.abc import x
>>> y = symbols('y', positive=True)
>>> eq = 1/(3*x*y + 1) + 4*x*y + 5
Dummy symbols with no name are used to keep track of the
constants. The retain unique identity, however:
>>> e, r = keepx(eq, x); e
_ + _*x + 1/(_ + _*x)
>>> assert e.xreplace(r) == eq
>>> sol = solve(eq, x)
>>> _sol = [i.xreplace(r).factor() for i in solve(e, x)]
>>> assert sol == _sol
"""
reps = {}
def store(i):
d = Dummy('')
reps[d] = i
return d
def do(e):
if e.is_Add or e.is_Mul:
i, d = e.as_independent(x)
if i is not e.identity:
i = store(i)
d = do(d)
return e.func(i, d)
if not e.args:
return e
return e.func(*[do(i) for i in e.args])
return do(eq), reps
要在您的案例中使用它来获得解决方案，只需执行以下操作：
>>> eq, x = (sp.Matrix([np.dot((jacobianmat-(k**2)*diffmatrix)[
0:n-1],gamma[0:n-1])])+sp.Matrix(doublesummationthirdorder1[0:n-
1])+sp.Matrix(3/math.factorial(4)*sp.Matrix(triplesummationthirdorder[0:n-
1])),gamma[0:n-1])
>>> e, r = keepx(eq, x)
>>> auxiliarterm = solve(e, x)[0].xreplace(r)
它仍然是一个笨拙的表达式，但现在您可以快速解决问题（以及您可以在需要屏蔽“噪音”以便能够使用表达式的核心时使用的例程） .
为了引入问题，我把我的代码写在下面：
import sympy as sp
sp.init_printing()
import numpy as np
import math
n=2
u,v=sp.symbols('u v')
a,b,c,d,h=sp.symbols('a b c d h')
k=sp.symbols('k')
diffmatrix=sp.zeros(n)
for i in range(n):
diffmatrix[i,i] = sp.symbols('D'+str(i+1))
globals()['D'+str(i+1)]=sp.symbols('D'+str(i+1))
f=[]
var=[]
var.append('u')
var.append('v')
# f.append(u+v)
# f.append(u-v**2)
f.append(a-c*u+d*v+u**2*v)
f.append(b+h*u-d*v-u**2*v)
f=sp.Matrix(f)
var=sp.Matrix(var)
jacobianmat=f.jacobian(var)
phi=[]
alpha=[]
beta=[]
psi=[]
gamma=[]
for i in range(n):
phi.append(sp.symbols('phi' + str(i)))
alpha.append(sp.symbols('alpha' + str(i)))
beta.append(sp.symbols('beta'+str(i)))
psi.append(sp.symbols('psi' + str(i)))
gamma.append(sp.symbols('gamma' + str(i)))
globals()['phi' + str(i)] = sp.symbols('phi' + str(i))
globals()['alpha' + str(i)] = sp.symbols('alpha' + str(i))
globals()['beta' + str(i)] = sp.symbols('beta' + str(i))
globals()['psi' + str(i)] = sp.symbols('psi' + str(i))
globals()['gamma' + str(i)] = sp.symbols('gamma' + str(i))
auxiliarterm=sp.linsolve(sp.Matrix([np.dot((jacobianmat-(k**2)*diffmatrix)[0:n-1],phi[0:n-1])])+sp.Matrix(jacobianmat-(k**2)*diffmatrix)[0:n-1,n-1],phi[0:n-1])
phi[0:n-1]=sp.Matrix(list(auxiliarterm))
phi[n-1]=1
doublesummationsecondorder=[]
for functionnumber in range(n):
tempsum=0
for i in range(n):
for j in range(n):
tempsum=sp.Add(tempsum,phi[i]*phi[j]*sp.diff(f[functionnumber],var[i],var[j]))
tempsum=-tempsum/4
doublesummationsecondorder.append(tempsum)
doublesummationsecondorder=sp.Matrix(doublesummationsecondorder)
alpha=sp.linsolve(sp.Matrix(np.dot(jacobianmat,alpha))-doublesummationsecondorder,alpha)
alpha=sp.Matrix(sp.Transpose(sp.Matrix(list(alpha))))
beta=sp.linsolve(sp.Matrix(np.dot(jacobianmat-4*k**2*diffmatrix,beta))-doublesummationsecondorder,beta)
beta=sp.Matrix(sp.Transpose(sp.Matrix(list(beta))))
doublesummationthirdorder1=[]
doublesummationthirdorder2=[]
triplesummationthirdorder=[]
for functionnumber in range(n):
tempsum1=0
tempsum2=0
tempsum3=0
for i in range(n):
for j in range(n):
tempsum1=sp.Add(tempsum1,phi[i]*(alpha[j]+beta[j]/2)*sp.diff(f[functionnumber],var[i],var[j]))
tempsum2=sp.Add(tempsum2,phi[i]*beta[j]*sp.diff(f[functionnumber],var[i],var[j]))
for l in range(n):
tempsum3=sp.Add(tempsum3,phi[i]*phi[j]*phi[l]*sp.diff(f[functionnumber],var[i],var[j],var[l]))
doublesummationthirdorder1.append(tempsum1)
doublesummationthirdorder2.append(tempsum2)
triplesummationthirdorder.append(tempsum3)
doublesummationthirdorder1=sp.Matrix(doublesummationthirdorder1)
doublesummationthirdorder2=sp.Matrix(doublesummationthirdorder2)
triplesummationthirdorder=sp.Matrix(triplesummationthirdorder)
auxiliarterm=sp.linsolve(sp.Matrix([np.dot((sp.matrices.Transpose(jacobianmat)-(k**2)*diffmatrix)[0:n-1],psi[0:n-1])])+sp.Matrix(sp.matrices.Transpose(jacobianmat)-(k**2)*diffmatrix)[0:n-1,n-1],psi[0:n-1])
psi[0:n-1]=sp.Matrix(list(auxiliarterm))
psi[n-1]=1
print('Almost ready')
auxiliarterm=sp.solve(sp.Matrix([np.dot((jacobianmat-(k**2)*diffmatrix)[0:n-1],gamma[0:n-1])])+sp.Matrix(doublesummationthirdorder1[0:n-1])+sp.Matrix(3/math.factorial(4)*sp.Matrix(triplesummationthirdorder[0:n-1])),gamma[0:n-1])
如您所见，当 n=2 时，最后求解的方程是一维（线性）方程。问题是，由于某种原因，我无法得到它的解决方案。在任何情况下，如果我尝试使用 f=[u+v,u-v^2] 给出的另一个向量场，那么一切正常。我想知道这是内存的问题，符号计算的效率还是我的代码有问题。
万事如意。
问题是求解器会尝试扩展您拥有的表达式，并且这些表达式在扩展时会变大，所以速度很慢。对于符号抵消的线性方程组，通常需要扩展。在这种情况下，因为它只是一个方程式，所以如果您确定存在唯一解，则可以跳过该步骤：
In [41]: a = Wild('a')
In [42]: b = Wild('b')
In [43]: match = eq[0].match(a*s[0]+b)
In [44]: sol = -match[b]/match[a]
In [45]: sol
Out[45]:
⎛ 2 2
⎛ 2⎞ ⎜ 2⋅D₁⋅k ⋅u + 2⋅c⋅u + d⋅v - 3⋅u ⋅v
2⋅u⋅⎝d + u ⎠⋅⎜───────────────────────────────────────────────────────────── + ──────────────────────────────────────────────────────────────────────────────
⎜ 2 4 2 2 2 2 2 ⎛ 3 8 3 6 3 6 2 2 6 2 6
⎝2⋅D₁ ⋅k + 4⋅D₁⋅c⋅k - 8⋅D₁⋅k ⋅u⋅v + 2⋅c - 8⋅c⋅u⋅v + 8⋅u ⋅v 2⋅⎝32⋅D₁ ⋅D₂⋅k + 8⋅D₁ ⋅d⋅k + 8⋅D₁ ⋅k ⋅u + 72⋅D₁ ⋅D₂⋅c⋅k - 144⋅D₁ ⋅D₂⋅k ⋅u⋅
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
2 4 2 4 3 2
8⋅D₁ ⋅d⋅k ⋅u + 8⋅D₁ ⋅k ⋅u + 10⋅D₁⋅c⋅d⋅k ⋅u +
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
2 4 2 4 2 2 4 2 4 2 4 2 2 4 3 2 4 4 4 2
v + 18⋅D₁ ⋅c⋅d⋅k + 18⋅D₁ ⋅c⋅k ⋅u - 2⋅D₁ ⋅d⋅h⋅k - 32⋅D₁ ⋅d⋅k ⋅u⋅v - 2⋅D₁ ⋅h⋅k ⋅u - 32⋅D₁ ⋅k ⋅u ⋅v + 48⋅D₁⋅D₂⋅c ⋅k - 192⋅D₁⋅D₂⋅c⋅k ⋅u⋅v + 192⋅D₁⋅D₂⋅k ⋅u
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
2 3 2 2 2 2 2 2 3 2 4 2 2 3 2 2
10⋅D₁⋅c⋅k ⋅u + 4⋅D₁⋅d ⋅k ⋅v - 2⋅D₁⋅d⋅h⋅k ⋅u - 8⋅D₁⋅d⋅k ⋅u ⋅v - 2⋅D₁⋅h⋅k ⋅u - 12⋅D₁⋅k ⋅u ⋅v + 2⋅c ⋅d⋅u + 2⋅c ⋅u + c⋅d ⋅v - 2⋅c⋅d⋅h⋅u - 2⋅c⋅d⋅u ⋅v - 2⋅c⋅h
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
2 2 2 2 2 2 2 2 2 2 2 3 2 2 2 2 2
⋅v + 12⋅D₁⋅c ⋅d⋅k + 12⋅D₁⋅c ⋅k ⋅u - 4⋅D₁⋅c⋅d⋅h⋅k - 40⋅D₁⋅c⋅d⋅k ⋅u⋅v - 4⋅D₁⋅c⋅h⋅k ⋅u - 40⋅D₁⋅c⋅k ⋅u ⋅v + 8⋅D₁⋅d⋅h⋅k ⋅u⋅v + 32⋅D₁⋅d⋅k ⋅u ⋅v + 8⋅D₁⋅h⋅k ⋅
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
2
D₁⋅k + c - 2⋅u⋅v
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
3 4 2 2 4
⋅u - 3⋅c⋅u ⋅v - d ⋅h⋅v + 2⋅d⋅h⋅u ⋅v + 3⋅h⋅u ⋅v
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
3 2 4 2 3 2 2 2 2 2 2 2 3 3 3 3 2 2 2 2 2 2 3
u ⋅v + 32⋅D₁⋅k ⋅u ⋅v + 8⋅D₂⋅c ⋅k - 48⋅D₂⋅c ⋅k ⋅u⋅v + 96⋅D₂⋅c⋅k ⋅u ⋅v - 64⋅D₂⋅k ⋅u ⋅v + 2⋅c ⋅d + 2⋅c ⋅u - 2⋅c ⋅d⋅h - 8⋅c ⋅d⋅u⋅v - 2⋅c ⋅h⋅u - 8⋅c ⋅u ⋅v
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
⎞
⎟
───────────────────────────────────────────────────────────────────────────────⎟
2 2 3 4 2 2 2 4 2⎞⎟
+ 8⋅c⋅d⋅h⋅u⋅v + 8⋅c⋅d⋅u ⋅v + 8⋅c⋅h⋅u ⋅v + 8⋅c⋅u ⋅v - 8⋅d⋅h⋅u ⋅v - 8⋅h⋅u ⋅v ⎠⎠
──────────────────────────────────────────────────────────────────────────────── - ─────────────────────────────────────────────────────────────────────────
3 8 3 6 3 6 2 2 6 2 6
16⋅D₁ ⋅D₂⋅k + 4⋅D₁ ⋅d⋅k + 4⋅D₁ ⋅k ⋅u + 36⋅D₁ ⋅D₂⋅c⋅k - 72⋅D₁ ⋅D₂⋅k ⋅u
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
2 4 2 4 2 2 4 2 4 2 4 2 2 4 3 2 4 4 4 2 2
⋅v + 9⋅D₁ ⋅c⋅d⋅k + 9⋅D₁ ⋅c⋅k ⋅u - D₁ ⋅d⋅h⋅k - 16⋅D₁ ⋅d⋅k ⋅u⋅v - D₁ ⋅h⋅k ⋅u - 16⋅D₁ ⋅k ⋅u ⋅v + 24⋅D₁⋅D₂⋅c ⋅k - 96⋅D₁⋅D₂⋅c⋅k ⋅u⋅v + 96⋅D₁⋅D₂⋅k ⋅u ⋅v + 6
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
⎛ 4 4 3 2 2 3 2 2 2 2 2 4 ⎞
u⋅⎝4⋅D₁⋅D₂⋅d⋅k ⋅u + 4⋅D₁⋅D₂⋅k ⋅u + 4⋅D₂⋅c⋅d⋅k ⋅u + 4⋅D₂⋅c⋅k ⋅u + 2⋅D₂⋅d ⋅k ⋅v - 4⋅D₂⋅d⋅k ⋅u ⋅v - 6⋅D₂⋅k ⋅u ⋅v⎠
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
2 2 2 2 2 2 2 2 2 2 3 2 2 2 2 2 3
⋅D₁⋅c ⋅d⋅k + 6⋅D₁⋅c ⋅k ⋅u - 2⋅D₁⋅c⋅d⋅h⋅k - 20⋅D₁⋅c⋅d⋅k ⋅u⋅v - 2⋅D₁⋅c⋅h⋅k ⋅u - 20⋅D₁⋅c⋅k ⋅u ⋅v + 4⋅D₁⋅d⋅h⋅k ⋅u⋅v + 16⋅D₁⋅d⋅k ⋅u ⋅v + 4⋅D₁⋅h⋅k ⋅u ⋅v + 16
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
2
- D₁⋅k - c + 2⋅u⋅v
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
2 4 2 3 2 2 2 2 2 2 2 3 3 3 3 2 2 2 2 2 2 3
⋅D₁⋅k ⋅u ⋅v + 4⋅D₂⋅c ⋅k - 24⋅D₂⋅c ⋅k ⋅u⋅v + 48⋅D₂⋅c⋅k ⋅u ⋅v - 32⋅D₂⋅k ⋅u ⋅v + c ⋅d + c ⋅u - c ⋅d⋅h - 4⋅c ⋅d⋅u⋅v - c ⋅h⋅u - 4⋅c ⋅u ⋅v + 4⋅c⋅d⋅h⋅u⋅v + 4
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
───────────────────────────────────────────────────────────── - ────────────────────────────────────────────────────────────────────────────────────────────
2 2 3 4 2 2 2 4 2 ⎛ 2 ⎞ ⎛ 3 8 3 6 3 6 2 2 6 2 6
⋅c⋅d⋅u ⋅v + 4⋅c⋅h⋅u ⋅v + 4⋅c⋅u ⋅v - 4⋅d⋅h⋅u ⋅v - 4⋅h⋅u ⋅v ⎝D₁⋅k + c - 2⋅u⋅v⎠⋅⎝16⋅D₁ ⋅D₂⋅k + 4⋅D₁ ⋅d⋅k + 4⋅D₁ ⋅k ⋅u + 36⋅D₁ ⋅D₂⋅c⋅k - 72⋅D₁ ⋅D₂⋅k
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
2 4 2 4 2 2 4 2 4 2 4 2 2 4 3 2 4 4 4 2 2
⋅u⋅v + 9⋅D₁ ⋅c⋅d⋅k + 9⋅D₁ ⋅c⋅k ⋅u - D₁ ⋅d⋅h⋅k - 16⋅D₁ ⋅d⋅k ⋅u⋅v - D₁ ⋅h⋅k ⋅u - 16⋅D₁ ⋅k ⋅u ⋅v + 24⋅D₁⋅D₂⋅c ⋅k - 96⋅D₁⋅D₂⋅c⋅k ⋅u⋅v + 96⋅D₁⋅D₂⋅k ⋅u ⋅v +
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
⎛ 2⎞ ⎛ 4 4 3 2 2 3 2 2 2 2 2 4 ⎞
v⋅⎝d + u ⎠⋅⎝4⋅D₁⋅D₂⋅d⋅k ⋅u + 4⋅D₁⋅D₂⋅k ⋅u + 4⋅D₂⋅c⋅d⋅k ⋅u + 4⋅D₂⋅c⋅k ⋅u + 2⋅D₂⋅d ⋅k ⋅v - 4⋅D₂⋅d⋅k ⋅u ⋅v - 6⋅D₂⋅k ⋅u ⋅v⎠
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
2 2 2 2 2 2 2 2 2 2 3 2 2 2 2 2 3
6⋅D₁⋅c ⋅d⋅k + 6⋅D₁⋅c ⋅k ⋅u - 2⋅D₁⋅c⋅d⋅h⋅k - 20⋅D₁⋅c⋅d⋅k ⋅u⋅v - 2⋅D₁⋅c⋅h⋅k ⋅u - 20⋅D₁⋅c⋅k ⋅u ⋅v + 4⋅D₁⋅d⋅h⋅k ⋅u⋅v + 16⋅D₁⋅d⋅k ⋅u ⋅v + 4⋅D₁⋅h⋅k ⋅u ⋅v +
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
2 4 2 3 2 2 2 2 2 2 2 3 3 3 3 2 2 2 2 2 2 3
16⋅D₁⋅k ⋅u ⋅v + 4⋅D₂⋅c ⋅k - 24⋅D₂⋅c ⋅k ⋅u⋅v + 48⋅D₂⋅c⋅k ⋅u ⋅v - 32⋅D₂⋅k ⋅u ⋅v + c ⋅d + c ⋅u - c ⋅d⋅h - 4⋅c ⋅d⋅u⋅v - c ⋅h⋅u - 4⋅c ⋅u ⋅v + 4⋅c⋅d⋅h⋅u⋅v +
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
2
⎛ 2⎞
0.75⋅⎝d + u ⎠
──────────────────────────────────────────────────────────────── - ────────────────────
2 2 3 4 2 2 2 4 2⎞ 2
4⋅c⋅d⋅u ⋅v + 4⋅c⋅h⋅u ⋅v + 4⋅c⋅u ⋅v - 4⋅d⋅h⋅u ⋅v - 4⋅h⋅u ⋅v ⎠ ⎛ 2 ⎞
⎝D₁⋅k + c - 2⋅u⋅v⎠
───────────────────────────────────────────────────────────────────────────────────────
另请注意，其中有一个 Float (0.75)。我会确保这里只使用有理数，因为浮点系数的多项式运算没有很好的定义。
正如 Oscar 指出的那样，表达式很大，求解过程在处理所有符号时陷入困境。以下例程将简化表达式，让您使用一小组符号：
def keepx(eq, x):
"""collapse non-x dependent Add and Mul into Duumy symbols
and return e, r where e is the new expresion and r is
the dictionary that can restore e: eq == e.xreplace(r)
Examples
========
>>> from sympy import solve, symbols
>>> from sympy.abc import x
>>> y = symbols('y', positive=True)
>>> eq = 1/(3*x*y + 1) + 4*x*y + 5
Dummy symbols with no name are used to keep track of the
constants. The retain unique identity, however:
>>> e, r = keepx(eq, x); e
_ + _*x + 1/(_ + _*x)
>>> assert e.xreplace(r) == eq
>>> sol = solve(eq, x)
>>> _sol = [i.xreplace(r).factor() for i in solve(e, x)]
>>> assert sol == _sol
"""
reps = {}
def store(i):
d = Dummy('')
reps[d] = i
return d
def do(e):
if e.is_Add or e.is_Mul:
i, d = e.as_independent(x)
if i is not e.identity:
i = store(i)
d = do(d)
return e.func(i, d)
if not e.args:
return e
return e.func(*[do(i) for i in e.args])
return do(eq), reps
要在您的案例中使用它来获得解决方案，只需执行以下操作：
>>> eq, x = (sp.Matrix([np.dot((jacobianmat-(k**2)*diffmatrix)[
0:n-1],gamma[0:n-1])])+sp.Matrix(doublesummationthirdorder1[0:n-
1])+sp.Matrix(3/math.factorial(4)*sp.Matrix(triplesummationthirdorder[0:n-
1])),gamma[0:n-1])
>>> e, r = keepx(eq, x)
>>> auxiliarterm = solve(e, x)[0].xreplace(r)
它仍然是一个笨拙的表达式，但现在您可以快速解决问题（以及您可以在需要屏蔽“噪音”以便能够使用表达式的核心时使用的例程） .