sympy 解多元方程组

#这里想用sympy解一个多元方程，但是运行以后进入循环了？（很久都没有结果）代码如下：

from sympy import *

V, W, X, Y, Z = 4, 4, 8, 5, 2

alpha = 1  # Volumetric shape factor of the particles
C_in = 249.971
Q_in = 0.0165
rho_NaCl = 2170  # NaCl density [kg m^-3]
rho_water = 1000  # Water density [kg m^-3]
T = 10 + 273.15  # System's temperature [K]
C_eq = 111.09  # Equilibrium concentration []
Volume = 50 + V  # Crystallizer volume [dm^3]

# Nucleation
A_n = (1 + X/10) * 10**8  # [m^-3 h^-1]
B_n = 30 + (Y/10)        # [K]

# Growth
A_g = (2 + X/10) * 10**-3  # [m h^-1]
B_g = 50 + (Y/10)         # [K]

# Setting Symbol for functions
S = Symbol('S')           # Supersaturation
M_V = Symbol('M_V')       # Mass of suspended crystals per volume
J = Symbol('J')           # Nucleation rate
#A_n = Symbol('A_n')
#B_n = Symbol('B_n')
#T = Symbol('T')
G = Symbol('G')           # Growth rate
#A_g = Symbol('A_g')
#B_g = Symbol('B_g')
#alpha = Symbol('alpha')
#C_in = Symbol('C_in')
C_out = Symbol('C_out')
#Q_in = Symbol('Q_in')
#C_eq = Symbol('C_eq')
#rho_NaCl = Symbol('rho_NaCl')
#rho_water = Symbol('rho_water')
#Volume = Symbol('Volume')
n_0 = Symbol('n_0')
tau = Symbol('tau')      # Residence time

solved_value = solve([S - C_out / C_eq,
                      M_V - (C_in - C_out),
                      M_V - (6 * rho_NaCl * alpha * n_0 * (G * tau)**4),
                      n_0 * G - J,
                      G - (A_g * (S - 1) * exp(-1*(B_g/T*(S-1)))),
                      J - (A_n * S * exp(-1 * (B_n/(T * (log(C_out / C_eq))**2)))),
                      tau * Q_in - Volume],
                     [S, M_V, J, G, C_out, n_0, tau])
print(solved_value)

• 帮你找了个相似的问题, 你可以看下: https://ask.csdn.net/questions/7767876
• 这篇博客你也可以参考下：sympy 求导，并带入值进行计算
• 除此之外, 这篇博客: Sympy代数符号运算库中的 方法 部分也许能够解决你的问题, 你可以仔细阅读以下内容或跳转源博客中阅读:
  

• evalf() 函数可以用求出代数式的浮点数。

  >>> (1/x).evalf(subs={x: 3.0}, n=21)
  0.333333333333333333333
  rather than
  >>> (1/x).subs({x: 3.0}).evalf(21)
  0.333333333333333314830
  >>> E.evalf(10)
  2.718281828
  >>> pi.evalf(15)
  3.14159265358979
  >>> sqrt(2).evalf()  #查看实数的有效数字，缺省15
  1.41421 35623 7310
  >>> sqrt(3).evalf()
  1.73205080756888
  

• 以下回答由chatgpt基于相关博客总结生成:

  要使用Sympy求解多元方程组，需要先将方程组转化成需要的形式，然后用Sympy提供的函数进行求解。具体步骤如下：

  1. 定义方程组

  首先需要将方程组写成符号表达式的形式，每个变量都用Symbol()函数定义，并将它们组成方程式的元组：

  import sympy as sp
  
  # 定义符号变量
  x1, x2, x3 = sp.symbols('x1 x2 x3')
  
  # 定义方程组
  eq1 = sp.Eq(3*x1 + 2*x2 - x3, 1)
  eq2 = sp.Eq(2*x1 - 2*x2 + 4*x3, -2)
  eq3 = sp.Eq(-x1 + (1/2)*x2 - x3, 0)
  
  # 将方程组写成元组的形式
  eqs = [eq1, eq2, eq3]
  

  1. 解方程组

  将方程组传入 solve() 函数，并指定需要求解的未知数，可以得到方程组的解：

  # 解多元方程组
  solutions = sp.solve(eqs, (x1, x2, x3))
  
  # 打印解
  print(solutions)
  

  指定需要求解的未知数是可选项，如果不指定，solve()函数返回的解将包含所有未知数的值。

  注意：如果方程组的解不唯一，solve()函数将返回字典。字典的key是每个未知数，value是该未知数的解。如果方程组无解，则返回空列表。

  1. 特殊情况处理

  在解多元方程组时，有一些特殊的情况需要考虑。比如，如果一个未知数出现在方程组的所有等式中，那么它的值可以直接从其中一条等式中得到。在Sympy中，可以使用subs()函数实现这个功能：

  # 解f(x)=g(x)中的未知数x
  x = sp.Symbol('x')
  f = x**2 - 4
  g = x - 2
  sol = sp.solve(sp.Eq(f, g), x)
  
  # 如果方程组的解不唯一，处理特殊情况
  if isinstance(sol, list):
      # 如果f(x)=g(x)中的未知数x出现在f和g中
      if x in f.free_symbols and x in g.free_symbols:
          # 用g(x)中的x替换f(x)中的x，得到一个一元方程式
          eq = f.subs(x, g)
          # 解一元方程式得到x的值
          sol = sp.solve(eq, x)
  

  如果方程组包含超越函数（如sin(x)、exp(x)等），则需要使用nsolve()函数。nsolve()函数需要指定超越函数的初始值，初始值可以通过实验得到：

  # 解 sin(x) = x
  x = sp.Symbol('x')
  equation = sp.sin(x) - x
  sol = sp.nsolve(equation, 0.5)
  

  这里指定初始值为0.5。

该回答引用GPT
你的代码中出现了一个循环，可能是由于符号解决方案的复杂性导致的。你可以尝试使用更具体的求解器，例如linsolve或nonlinsolve，或者尝试简化方程组。以下是一个简化的示例代码：

from sympy import *

# 设置符号
S, M_V, J, G, C_out, n_0, tau = symbols('S M_V J G C_out n_0 tau')

# 定义方程
eq1 = Eq(S - C_out / 111.09)
eq2 = Eq(M_V - (249.971 - C_out))
eq3 = Eq(M_V - (6 * 2170 * 1 * n_0 * (G * tau)**4))
eq4 = Eq(n_0 * G - J)
eq5 = Eq(G - ((2 + 0.1 * 8) * 10**-3 * (S - 1) * exp(-1*((50 + 0.1 * 5)/T*(S-1)))))
eq6 = Eq(J - ((1 + 0.1 * 4) * 10**8 * S * exp(-1 * ((30 + 0.1 * 5)/(T * (log(C_out / 111.09))**2)))))
eq7 = Eq(tau * 0.0165 - (50 + 4))

# 解方程
sol = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7], [S, M_V, J, G, C_out, n_0, tau])

# 输出解
print(sol)

您的方程组可能很复杂，这可能导致SymPy在求解时花费很长时间。尝试对方程进行简化或者使用数值方法求解。可以使用SymPy的nsolve函数来尝试数值求解。首先，需要为您的变量提供一个初始猜测值。以下是如何使用nsolve求解方程组的示例：

from sympy import symbols, exp, log, nsolve

V, W, X, Y, Z = 4, 4, 8, 5, 2

alpha = 1  # Volumetric shape factor of the particles
C_in = 249.971
Q_in = 0.0165
rho_NaCl = 2170  # NaCl density [kg m^-3]
rho_water = 1000  # Water density [kg m^-3]
T = 10 + 273.15  # System's temperature [K]
C_eq = 111.09  # Equilibrium concentration []
Volume = 50 + V  # Crystallizer volume [dm^3]

# Nucleation
A_n = (1 + X/10) * 10**8  # [m^-3 h^-1]
B_n = 30 + (Y/10)        # [K]

# Growth
A_g = (2 + X/10) * 10**-3  # [m h^-1]
B_g = 50 + (Y/10)         # [K]

# Setting Symbol for functions
S, M_V, J, G, C_out, n_0, tau = symbols('S M_V J G C_out n_0 tau')

# Define the equations
equations = [
    S - C_out / C_eq,
    M_V - (C_in - C_out),
    M_V - (6 * rho_NaCl * alpha * n_0 * (G * tau)**4),
    n_0 * G - J,
    G - (A_g * (S - 1) * exp(-1*(B_g/T*(S-1)))),
    J - (A_n * S * exp(-1 * (B_n/(T * (log(C_out / C_eq))**2)))),
    tau * Q_in - Volume
]

# Provide an initial guess
initial_guess = [1, 1, 1, 1, 1, 1, 1]

# Solve the equations numerically
solutions = nsolve(equations, (S, M_V, J, G, C_out, n_0, tau), initial_guess)

print(solutions)

这个例子定义了您的方程组，提供了一个初始猜测值，并使用nsolve函数对方程组进行了数值求解。根据这个解决方案，您应该能够获得更快的求解结果。注意，您可能需要调整初始猜测值以获得更好的收敛。