# 将 j 中的 d 替换为 1 s = {i: j.subs(d, 1) for i, j in solution.items()} # 将 {d: 1} 放进字典 s.update({d: 1})
# ========================== # 如果不知道 d 是自由变量,也可以求出来
# 选一个存在于解中的变量 m = list(solution.keys())[0] # 求出自由变量 n = list(solution[m].free_symbols)[0] # 求出特解 s = {i: j.subs(n, 1) for i, j in solution.items()} s.update({n: 1})
当方程组的解不为整数时,我们也可以求出一组最小(互质)的整数解:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
# python >= 3.9 才可用 math.lcm,否则须手动实现 # 又吃到了 python 的语法糖~ from math import lcm from sympy import *
# 求解 ...
# 求出分母 fm = [] for i in solution.values(): # 将分母加入数组 fm.append(fraction(i)[1])
# 计算最小公倍数 x = lcm(*fm) solution_ = {i: int(j*x) for i, j in s.items()}
defparse(s: str) -> dict: """ 解析化学式,计算原子数 s: 化学式 返回 dict 型 """ e = {} # 储存每种原子的个数 i = 0
while i < (l := len(s)): if s[i] == "(": # 要找到和当前左括号对应的右括号 j = i + 1 k = 1# 括号栈长度 while j < l and k > 0: if s[j] == "(": k += 1 elif s[j] == ")": k -= 1 j += 1 # 递归处理括号内 kh = parse(s[i + 1 : j - 1])
# 括号后的数字 kh_num_ = "" while j < l and s[j].isdigit(): kh_num_ += s[j] j += 1 kh_num = int(kh_num_) if kh_num_ else1
# 合并到 e 字典中 for p, q in kh.items(): e[p] = e.get(p, 0) + q * kh_num # 跳过已处理的括号和数字 i = j
elif s[i].isupper(): # 原子符号 el = s[i] i += 1 while i < l and s[i].islower(): el += s[i] i += 1
# 原子数(方法同上) el_num_ = "" while i < l and s[i].isdigit(): el_num_ += s[i] i += 1 el_num = int(el_num_) if el_num_ else1
# 合并到 e 字典中 e[el] = e.get(el, 0) + el_num
else: i += 1
return e
回到主函数,
1 2 3 4 5 6 7 8 9 10 11 12 13
# 解析化学式,计算原子数 tot = {} for i inrange(lr): z = parse(r[i]) for k, v in z.items(): if k notin tot: tot[k] = [f"{v} * r{i}", ""] else: tot[k][0] += f" + {v} * r{i}" for i inrange(lp): z = parse(p[i]) for k, v in z.items(): tot[k][1] += f"{ ' + 'if tot[k][1] else'' }{v} * p{i}"
然后建立方程组,
1 2 3 4 5 6
# 建立方程组 eqs = [] for k, v in tot.items(): left_expr = eval(v[0], {"__builtins__": None}, sym_dict) right_expr = eval(v[1], {"__builtins__": None}, sym_dict) eqs.append(Eq(left_expr, right_expr))
from math import lcm from sympy import symbols, Eq, solve, fraction
defparse(s: str) -> dict: """ 解析化学式,计算原子数 s: 化学式 返回 dict 型 """ e = {} # 储存每种原子的个数 i = 0
while i < (l := len(s)): if s[i] == "(": # 要找到和当前左括号对应的右括号 j = i + 1 k = 1# 括号栈长度 while j < l and k > 0: if s[j] == "(": k += 1 elif s[j] == ")": k -= 1 j += 1 # 递归处理括号内 kh = parse(s[i + 1 : j - 1])
# 括号后的数字 kh_num_ = "" while j < l and s[j].isdigit(): kh_num_ += s[j] j += 1 kh_num = int(kh_num_) if kh_num_ else1
# 合并到 e 字典中 for p, q in kh.items(): e[p] = e.get(p, 0) + q * kh_num # 跳过已处理的括号和数字 i = j
elif s[i].isupper(): # 原子符号 el = s[i] i += 1 while i < l and s[i].islower(): el += s[i] i += 1
# 原子数(方法同上) el_num_ = "" while i < l and s[i].isdigit(): el_num_ += s[i] i += 1 el_num = int(el_num_) if el_num_ else1
# 合并到 e 字典中 e[el] = e.get(el, 0) + el_num
else: i += 1
return e
defmain(): # 输入并分割化学式 input_ = lambda x: [j for i ininput(f"{x} >>> ").split("+") if (j := i.strip())] r = input_("反应物") p = input_("生成物")
# 生成未知数 lr = len(r) lp = len(p) rExp = "" pExp = "" for i inrange(lr): rExp += f"r{i} " for i inrange(lp): pExp += f"p{i} "
# 创建符号变量字典,供 eval 使用 sym_dict = {str(sym): sym for sym in (*rs, *ps)}
# 解析化学式,计算原子数 tot = {} for i inrange(lr): z = parse(r[i]) for k, v in z.items(): if k notin tot: tot[k] = [f"{v} * r{i}", ""] else: tot[k][0] += f" + {v} * r{i}" for i inrange(lp): z = parse(p[i]) for k, v in z.items(): tot[k][1] += f"{ ' + 'if tot[k][1] else'' }{v} * p{i}"
# 建立方程组 eqs = [] for k, v in tot.items(): left_expr = eval(v[0], {"__builtins__": None}, sym_dict) right_expr = eval(v[1], {"__builtins__": None}, sym_dict) eqs.append(Eq(left_expr, right_expr))
# 解方程组,求出通解 solution = solve(eqs, (*rs, *ps))
# 选一个存在于解中的变量 m = list(solution.keys())[0]
# 求出自由变量 n = list(solution[m].free_symbols)[0]
# 求出特解 s = {i: j.subs(n, 1) for i, j in solution.items()} s.update({n: 1})
# 求出分母 fm = [fraction(i)[1] for i in solution.values()]
# 计算最小公倍数 x = lcm(*fm) solution_ = {i: int(j * x) for i, j in s.items()}
# 输出结果 r_f = " + ".join( [ f"{ x_ if (x_ := solution_[sym]) != 1else'' }{r[i]}" for i, sym inenumerate(rs) ] ) p_f = " + ".join( [ f"{ x_ if (x_ := solution_[sym]) != 1else'' }{p[i]}" for i, sym inenumerate(ps) ] ) print(f"{r_f} == {p_f}")
# 创建变量 vars = sp.symbols(f'a0:{num_compounds}', positive=True, real=True)
# 建立方程组 equations = [] for elem in all_elements: left = sum(vars[i] * parsed[i].get(elem, 0) for i inrange(len(reactants))) right = sum(vars[len(reactants)+i] * parsed[len(reactants)+i].get(elem, 0) for i inrange(len(products))) equations.append(sp.Eq(left, right))
# 固定第一个系数为1 equations.append(sp.Eq(vars[0], 1))
# 解方程 solution = sp.solve(equations, vars)
# 处理可能缺失的变量 coeffs = [solution.get(v, 0) for v invars]
# 把分数转为整数 denominators = [sp.fraction(val)[1] for val in coeffs if val != 0] if denominators: lcm = sp.lcm(denominators) coeffs = [int(val * lcm) for val in coeffs]
# 解析元素和原子数 elements = {} # 匹配元素符号(可能是一个大写字母 followed by 小写字母)和可选的数字 pattern = r'([A-Z][a-z]*)(\d*)' matches = re.findall(pattern, compound)
for elem, count in matches: if count == '': count = 1 else: count = int(count) if elem in elements: elements[elem] += count else: elements[elem] = count
return elements
defget_all_elements(reactants, products): """获取所有涉及的元素""" elements = set() for comp in reactants + products: elements.update(comp.keys()) returnsorted(elements)
# 初始化矩阵 matrix = [[0for _ inrange(num_compounds)] for _ inrange(num_elements)]
# 填充反应物(正系数) for i, elem inenumerate(elements): for j, comp inenumerate(reactants): matrix[i][j] = comp.get(elem, 0)
# 填充产物(负系数,因为移到等号左边) for i, elem inenumerate(elements): for j, comp inenumerate(products): matrix[i][len(reactants) + j] = -comp.get(elem, 0)
return matrix
defbalance_equation(reactants_list, products_list): """配平化学方程式""" # 解析反应物和产物 reactants = [parse_compound(r) for r in reactants_list] products = [parse_compound(p) for p in products_list] elements = get_all_elements(reactants, products)