时间戳

首页 > 解决方案 > 来自 sympy 的命令 linsolve 或求解无法求解一维线性方程

python - 来自 sympy 的命令 linsolve 或求解无法求解一维线性方程

问题描述

我在下面编写代码以介绍问题：

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,uv^2] 给出的另一个向量场，那么一切正常。我想知道这是内存问题，符号计算的效率还是我的代码有问题。

一切都将不胜感激。

标签： pythonmatrixsympylinear-algebra

解决方案

正如 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)

它仍然是一个笨拙的表达式，但现在您有一个快速解决问题的方法（以及当您必须屏蔽“噪音”以便能够使用表达式的核心时可以使用的例程）。

推荐阅读