使用SymPy求解含未定常数的四元方程组无结果问题求助
SymPy求解四元方程组返回空解的问题分析与解决
问题背景
尝试用SymPy求解含4个未知量(k1、k2、k3、k4)的四元线性方程组,方程组包含多个未定常数,但运行代码后控制台仅输出Solution: [],无法得到预期的、以其他参数表示的k1-k4的表达式。运行环境为Spyder,执行结果如下:
In[54]: runfile('E:/Spyder/eigen_functions.py', wdir='E:/Spyder')
Solution: []
原代码重现
from sympy.solvers import solve from sympy import symbols, Eq, exp k1, k2, k3, k4, x, a, b, d, h, l,t = symbols(" k1 k2 k3 k4 x a b d h l t ") equation_1 = Eq((k1*exp(-a*t*1j) + k2*exp(-b*t*1j)+ k3*exp(a*t*1j) + k4*exp(b*t*1j)) ,1) equation_2 = Eq(((1/(2*l*h*(x+1)))*(-k1*((x**2)-1-2*(a*x)+(d**2)+(l**2)+a**2)*exp(-a*t*1j) -k2*((x**2)-1-2*(b*x)+(d**2)+(l**2)+b**2)*exp(-b*t*1j) -k3*((x**2)-1+2*(a*x)+(d**2)+(l**2)+a**2)*exp(a*t*1j) -k4*((x**2)-1+2*(b*x)+(d**2)+(l**2)+b**2)*exp(b*t*1j))) , 0) equation_3 = Eq(((1/(2*d*h*l*h*(x+1)))*(-k1*((l**2)*(x+1+a)-(d**2)*(x+1-a)-(x**3)-(x**2)+x+1-a*((x**2)+1+2*x)-(a**2)*(x-1)-a**3)*exp(-a*t*1j) -k2*((l**2)*(x+1+b)-(d**2)*(x+1-b)-(x**3)-(x**2)+x+1-b*((x**2)+1+2*x)-(b**2)*(x-1)-b**3)*exp(-b*t*1j) -k3*((l**2)*(x+1-a)-(d**2)*(x+1+a)-(x**3)-(x**2)+x+1+a*((x**2)+1+2*x)-(a**2)*(x-1)+a**3)*exp(a*t*1j) -k4*((l**2)*(x+1-b)-(d**2)*(x+1+b)-(x**3)-(x**2)+x+1+b*((x**2)+1+2*x)-(b**2)*(x-1)+b**3)*exp(b*t*1j))) ,0) equation_4 = Eq(((1/(2*d*h*(x+1)))*(-k1*(((d**2)+(l**2)-(a**2)-(x**2)+1+2*a*x))*exp(-a*t*1j) -k2*(((d**2)+(l**2)-(b**2)-(x**2)+1+2*b*x))*exp(-b*t*1j) -k3*(((d**2)+(l**2)-(a**2)-(x**2)+1-2*a*x))*exp(a*t*1j) -k4*(((d**2)+(l**2)-(b**2)-(x**2)+1-2*b*x))*exp(b*t*1j))) , 0) solution = solve((equation_1, equation_2, equation_3, equation_4), (k1, k2, k3, k4)) print("Solution:", solution)
问题原因及解决步骤
1. 消除冗余分母,简化方程组
方程2、3、4的右侧均为0,只要分母不为0(即x≠-1、l≠0、h≠0、d≠0),可以直接去掉分母,仅保留分子等于0,大幅简化计算量,避免SymPy因处理复杂分式而无法识别线性结构。
2. 利用指数项的线性独立性拆分方程
在a≠b、a≠0、b≠0的前提下,exp(-a*t*1j)、exp(-b*t*1j)、exp(a*t*1j)、exp(b*t*1j)是线性无关的项。因此原方程组可以拆分为每个指数项的系数分别等于对应右侧的系数:
- 方程1中,各指数项的系数对应k1、k2、k3、k4,右侧常数项为1;
- 方程2、3、4中,每个指数项的系数均为0。
3. 使用线性方程组求解器linsolve
原方程组是关于k1-k4的线性方程组,使用SymPy的linsolve比通用的solve更高效,专门处理线性系统,更容易得到解析解。
修改后的代码
from sympy import symbols, Eq, exp, linsolve, simplify # 定义符号,可添加约束(如a、b非零且不等) k1, k2, k3, k4, x, a, b, d, h, l, t = symbols("k1 k2 k3 k4 x a b d h l t", nonzero=True) a, b = symbols('a b', real=True, distinct=True) # -------------------------- # 步骤1:简化方程,去掉非零分母 # -------------------------- # 方程1保持不变 eq1 = Eq(k1*exp(-a*t*1j) + k2*exp(-b*t*1j) + k3*exp(a*t*1j) + k4*exp(b*t*1j), 1) # 方程2:提取分子并化简 coeff_k1_eq2 = simplify(-(x**2 - 1 - 2*a*x + d**2 + l**2 + a**2)) coeff_k2_eq2 = simplify(-(x**2 - 1 - 2*b*x + d**2 + l**2 + b**2)) coeff_k3_eq2 = simplify(-(x**2 - 1 + 2*a*x + d**2 + l**2 + a**2)) coeff_k4_eq2 = simplify(-(x**2 - 1 + 2*b*x + d**2 + l**2 + b**2)) eq2 = Eq(coeff_k1_eq2*k1*exp(-a*t*1j) + coeff_k2_eq2*k2*exp(-b*t*1j) + coeff_k3_eq2*k3*exp(a*t*1j) + coeff_k4_eq2*k4*exp(b*t*1j), 0) # 方程3:提取分子并化简 coeff_k1_eq3 = simplify(-(l**2*(x+1+a) - d**2*(x+1-a) - x**3 - x**2 + x + 1 - a*(x**2 + 1 + 2*x) - a**2*(x-1) - a**3)) coeff_k2_eq3 = simplify(-(l**2*(x+1+b) - d**2*(x+1-b) - x**3 - x**2 + x + 1 - b*(x**2 + 1 + 2*x) - b**2*(x-1) - b**3)) coeff_k3_eq3 = simplify(-(l**2*(x+1-a) - d**2*(x+1+a) - x**3 - x**2 + x + 1 + a*(x**2 + 1 + 2*x) - a**2*(x-1) + a**3)) coeff_k4_eq3 = simplify(-(l**2*(x+1-b) - d**2*(x+1+b) - x**3 - x**2 + x + 1 + b*(x**2 + 1 + 2*x) - b**2*(x-1) + b**3)) eq3 = Eq(coeff_k1_eq3*k1*exp(-a*t*1j) + coeff_k2_eq3*k2*exp(-b*t*1j) + coeff_k3_eq3*k3*exp(a*t*1j) + coeff_k4_eq3*k4*exp(b*t*1j), 0) # 方程4:提取分子并化简 coeff_k1_eq4 = simplify(-(d**2 + l**2 - a**2 - x**2 + 1 + 2*a*x)) coeff_k2_eq4 = simplify(-(d**2 + l**2 - b**2 - x**2 + 1 + 2*b*x)) coeff_k3_eq4 = simplify(-(d**2 + l**2 - a**2 - x**2 + 1 - 2*a*x)) coeff_k4_eq4 = simplify(-(d**2 + l**2 - b**2 - x**2 + 1 - 2*b*x)) eq4 = Eq(coeff_k1_eq4*k1*exp(-a*t*1j) + coeff_k2_eq4*k2*exp(-b*t*1j) + coeff_k3_eq4*k3*exp(a*t*1j) + coeff_k4_eq4*k4*exp(b*t*1j), 0) # -------------------------- # 步骤2:用linsolve求解线性方程组 # -------------------------- solution = linsolve((eq1, eq2, eq3, eq4), (k1, k2, k3, k4)) # 格式化输出解 print("Solution:") for idx, sol in enumerate(solution): print(f"解{idx+1}:") print(f"k1 = {simplify(sol[0])}") print(f"k2 = {simplify(sol[1])}") print(f"k3 = {simplify(sol[2])}") print(f"k4 = {simplify(sol[3])}")
额外优化建议
- 对系数表达式使用
sympy.simplify()处理,进一步减少计算复杂度; - 定义符号时添加约束条件(如
a、b为实数且不相等),帮助SymPy更高效地推导解析解; - 若解的表达式过于冗长,可使用
sympy.factor()或sympy.collect()对结果进行整理。
内容的提问来源于stack exchange,提问作者Lavoisier
相关产品推荐
相关产品推荐

