本文较长,没心情看的请直接跳到完整代码处。

起源

众所周知,我消失了十八天。事实上,我在育本集中营浪费了十八天的美丽青春——知识似乎没学会多少,但是我成功见识到了高一化学氧化还原反应方程式配平的难度。
文章中配平方程式的思想,主要来自于化学教书匠所谓“邪修”。“邪修”是指将一种物质的计量数记为11,依据原子守恒逐一推出每种物质的计量数,有时需要列方程。
“邪修”虽好,十分通用,但有时候计算量相当大,如配平下面这个方程式,“邪修”就远不及算转移电子数效率高:

3 Fe3O4+28 HNO3=9 Fe(NO3)3+NO+14 H2O3\ \text{Fe}_3\text{O}_4+28\ \text{HNO}_{3} = 9\ \text{Fe(NO}_3\text{)}_3+\text{NO}\uparrow+14\ \text{H}_2\text{O}

人脑的计算力不够,但是电脑足够。那天晚上,我在宿舍萌生了写这篇文章的想法,一直拖到回家。

前置知识:
依据质量守恒定律,化学反应前后,原子种类、数目、质量均不发生改变,即原子守恒。

“邪修”稍加变动,可以理解为:对于每种元素的原子,都有左、右两边原子数相等。利用待定系数法,设出每一个计量数,就得到了一个方程组,求解并选出最小整数解即可。

举个例子,对于简单的(初中的)置换反应Fe+HClFeCl2+H2\text{Fe}+\text{HCl} \to \text{FeCl}_2+\text{H}_2\uparrow,我们可以按照这种方式进行配平:

设每一项系数:

a Fe+b HCl=c FeCl2+d H2a\ \text{Fe}+b\ \text{HCl} = c\ \text{FeCl}_2+d\ \text{H}_2\uparrow

然后对于每种原子列出方程:

{(Fe)  a=c(H)    b=2d(Cl)  b=2c\left\{\begin{matrix} \texttt{(Fe)}\ \ a=c\\ \texttt{(H)}\ \ \ \ b=2d\\ \texttt{(Cl)}\ \ b=2c \end{matrix}\right.

然后消元即可得到通解:

{a=db=2dc=d\left\{\begin{matrix} a=d\\ b=2d\\ c=d \end{matrix}\right.

d=1d=1,则

{a=1b=2c=1d=1\left\{\begin{matrix} a=1\\ b=2\\ c=1\\ d=1 \end{matrix}\right.

代回原方程式,再省略掉系数11,即成功配平的方程式:

Fe+2 HCl=FeCl2+H2\text{Fe}+2\ \text{HCl} = \text{FeCl}_2+\text{H}_2\uparrow


理论

通过网上冲浪,我了解到 sympy 是一个可以进行数学符号计算第三方 Python 模块,例如解方程组。

比如,我希望利用它求上面方程组{a=cb=2db=2c\left\{\begin{matrix}a=c\\b=2d\\b=2c\end{matrix}\right. 的通解。可以这样做:

1
2
3
4
5
6
7
8
9
10
11
12
13
from sympy import symbols, Eq, solve

# 四个未知数
a, b, c, d = symbols('a b c d')

# 三个方程
eq1 = Eq(a, c)
eq2 = Eq(b, 2*d)
eq3 = Eq(b, 2*c)

# 解方程组,求出通解
solution = solve((eq1, eq2, eq3), (a, b, c, d))
print(solution)

程序输出字典 {a: d, b: 2*d, c: d},即方程组通解为:{a=db=2dc=d\left\{\begin{matrix}a=d\\b=2d\\c=d\end{matrix}\right..

d=1d=1,则得到特解:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
# 将 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()}

稍加整合,我们就可以完成整个计算部分的程序了。


实践

承接例子

承接上面的例子,依然是对于铁和稀盐酸的置换反应,我们可以写出计算程序:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
from math import lcm
from sympy import symbols, Eq, solve, fraction

# 四个未知数
a, b, c, d = symbols('a b c d')

# 三个方程
eq1 = Eq(a, c)
eq2 = Eq(b, 2*d)
eq3 = Eq(b, 2*c)

# 解方程组,求出通解
solution = solve((eq1, eq2, eq3), (a, b, c, 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})

# 求出分母
fm = [fraction(i)[1] for i in solution.values()]

# 计算最小公倍数
x = lcm(*fm)
solution_ = {i: int(j*x) for i, j in s.items()}

print(solution_)

程序成功输出 {a: 1, b: 2, c: 1, d: 1}

拓展到一般情况

首先,从控制台读取并处理用户的输入:

1
2
3
4
# 输入并分割化学式
input_ = lambda x: [j for i in input(f"{x} >>> ").split("+") if (j := i.strip())]
r = input_("反应物")
p = input_("生成物")

生成对应的未知数,其中反应物的为r1,r2,...,rlr1r_1,r_2,...,r_{\text{lr}-1},生成物的为p1,p2,...,plp1p_1,p_2,...,p_{\text{lp}-1}

1
2
3
4
5
6
7
8
9
# 生成未知数
lr = len(r)
lp = len(p)
rExp = ""
pExp = ""
for i in range(lr):
rExp += f"r{i} "
for i in range(lp):
pExp += f"p{i} "

然后设出未知数:

1
2
3
4
5
6
7
8
9
10
# 设出未知数,确保返回的是元组形式
rs = symbols(rExp.strip())
ps = symbols(pExp.strip())
if not isinstance(rs, tuple):
rs = (rs,)
if not isinstance(ps, tuple):
ps = (ps,)

# 创建符号变量字典,供 eval 使用
sym_dict = {str(sym): sym for sym in (*rs, *ps)}

接下来,解析化学式,计算原子数。这里我设计了一个 parse 方法。由于我还不会正则,所以处理格式问题很麻烦,用正则会大幅减少工作量。

所以这一段写得我手痒痒的,感觉 python 这么写很啰嗦,不如 C++ 了

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
def parse(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_ else 1

# 合并到 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_ else 1

# 合并到 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 in range(lr):
z = parse(r[i])
for k, v in z.items():
if k not in tot:
tot[k] = [f"{v} * r{i}", ""]
else:
tot[k][0] += f" + {v} * r{i}"
for i in range(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))

这几步和例子的解法不太一样,是因为面对未知,程序需要自动地动态生成符号变量,用迭代的方法来处理解决未知问题。

最后是输出:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
# 输出结果
r_f = " + ".join(
[
f"{ x_ if (x_ := solution_[sym]) != 1 else '' }{r[i]}"
for i, sym in enumerate(rs)
]
)
p_f = " + ".join(
[
f"{ x_ if (x_ := solution_[sym]) != 1 else '' }{p[i]}"
for i, sym in enumerate(ps)
]
)
print(f"{r_f} == {p_f}")

完整代码

注意,此程序要求使用高于 Python3.93.9 的版本,否则会报错!(低于此版本的可以自行敲一下 lcm 函数,辗转相除法等可以手动实现~)

如果没有安装过 sympy,请先 pip 安装:

1
pip install sympy

以下为完整的程序:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
from math import lcm
from sympy import symbols, Eq, solve, fraction


def parse(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_ else 1

# 合并到 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_ else 1

# 合并到 e 字典中
e[el] = e.get(el, 0) + el_num

else:
i += 1

return e


def main():
# 输入并分割化学式
input_ = lambda x: [j for i in input(f"{x} >>> ").split("+") if (j := i.strip())]
r = input_("反应物")
p = input_("生成物")

# 生成未知数
lr = len(r)
lp = len(p)
rExp = ""
pExp = ""
for i in range(lr):
rExp += f"r{i} "
for i in range(lp):
pExp += f"p{i} "

# 设出未知数,确保返回的是元组形式
rs = symbols(rExp.strip())
ps = symbols(pExp.strip())
if not isinstance(rs, tuple):
rs = (rs,)
if not isinstance(ps, tuple):
ps = (ps,)

# 创建符号变量字典,供 eval 使用
sym_dict = {str(sym): sym for sym in (*rs, *ps)}

# 解析化学式,计算原子数
tot = {}
for i in range(lr):
z = parse(r[i])
for k, v in z.items():
if k not in tot:
tot[k] = [f"{v} * r{i}", ""]
else:
tot[k][0] += f" + {v} * r{i}"
for i in range(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]) != 1 else '' }{r[i]}"
for i, sym in enumerate(rs)
]
)
p_f = " + ".join(
[
f"{ x_ if (x_ := solution_[sym]) != 1 else '' }{p[i]}"
for i, sym in enumerate(ps)
]
)
print(f"{r_f} == {p_f}")


if __name__ == "__main__":
main()

附赠两个 AI 写的配平方程式的程序:

这个是使用正则表达式实现的:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
import sympy as sp
import re

def parse_formula(formula):
"""解析化学式,返回元素和对应原子数的字典"""
elements = {}

# 处理括号的情况,比如Ca(OH)2
while '(' in formula:
match = re.search(r'\(([A-Za-z0-9]+)\)(\d*)', formula)
if not match:
break

inner = match.group(1) # 括号里的内容
num = int(match.group(2)) if match.group(2) else 1 # 括号后的数字

# 递归处理括号里的内容
inner_elements = parse_formula(inner)

# 把括号内的元素乘以系数
for elem, count in inner_elements.items():
elements[elem] = elements.get(elem, 0) + count * num

formula = formula.replace(match.group(0), '')

# 处理普通元素
matches = re.findall(r'([A-Z][a-z]?)(\d*)', formula)
for elem, num_str in matches:
count = int(num_str) if num_str else 1
elements[elem] = elements.get(elem, 0) + count

return elements

def balance_equation(reactants, products):
"""配平化学方程式的主函数"""
all_compounds = reactants + products
all_elements = set()

# 解析所有化学式
parsed = []
for compound in all_compounds:
p = parse_formula(compound)
parsed.append(p)
all_elements.update(p.keys())

all_elements = list(all_elements)
num_compounds = len(all_compounds)

# 创建变量
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 in range(len(reactants)))
right = sum(vars[len(reactants)+i] * parsed[len(reactants)+i].get(elem, 0) for i in range(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 in vars]

# 把分数转为整数
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]

return coeffs

def main():
print("=== 化学方程式配平程序 ===")
print("输入说明:")
print("1. 多个反应物/生成物用逗号分隔")
print("2. 化学式区分大小写,如CO和Co不同")
print("3. 支持括号,如Ca(OH)2")
print("--------------------------")

reactants = input("请输入反应物: ").replace(" ", "").split(',')
products = input("请输入生成物: ").replace(" ", "").split(',')

try:
coeffs = balance_equation(reactants, products)

# 构建反应式
def format_side(compounds, start, end):
parts = []
for i in range(start, end):
coeff = coeffs[i]
parts.append(f"{coeff if coeff !=1 else ''}{compounds[i-start]}")
return " + ".join(parts)

left = format_side(reactants, 0, len(reactants))
right = format_side(products, len(reactants), len(coeffs))

print("\n配平结果:")
print(f"{left}{right}")

except Exception as e:
print(f"\n错误: {e}")

if __name__ == "__main__":
main()

另外一个是使用矩阵实现的,似乎是个数学概念。但是这个东西我还没系统地学,毕竟我也只是一个准高一学生嘛。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
import re
import numpy as np
from fractions import Fraction
from sympy import Matrix, lcm

def parse_compound(compound):
"""解析化合物,返回元素及其原子数的字典"""
# 处理括号,例如 (OH)2
while '(' in compound:
# 找到最内层的括号
start = compound.rfind('(')
end = compound.find(')', start)
if end == -1:
break # 不匹配的括号,直接返回

# 获取括号内的内容和系数
inner = compound[start+1:end]
# 查找括号后的数字
num_match = re.match(r'\d+', compound[end+1:])
if num_match:
num = int(num_match.group())
compound = compound[:start] + inner*num + compound[end+1+len(num_match.group()):]
else:
compound = compound[:start] + inner + compound[end+1:]

# 解析元素和原子数
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

def get_all_elements(reactants, products):
"""获取所有涉及的元素"""
elements = set()
for comp in reactants + products:
elements.update(comp.keys())
return sorted(elements)

def build_matrix(reactants, products, elements):
"""构建矩阵方程 Ax = 0"""
num_elements = len(elements)
num_compounds = len(reactants) + len(products)

# 初始化矩阵
matrix = [[0 for _ in range(num_compounds)] for _ in range(num_elements)]

# 填充反应物(正系数)
for i, elem in enumerate(elements):
for j, comp in enumerate(reactants):
matrix[i][j] = comp.get(elem, 0)

# 填充产物(负系数,因为移到等号左边)
for i, elem in enumerate(elements):
for j, comp in enumerate(products):
matrix[i][len(reactants) + j] = -comp.get(elem, 0)

return matrix

def balance_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)

# 构建矩阵
matrix = build_matrix(reactants, products, elements)
num_compounds = len(reactants) + len(products)

# 如果矩阵行数大于列数,取前num_compounds-1行
if len(matrix) > num_compounds - 1:
matrix = matrix[:num_compounds - 1]

# 使用sympy的Matrix求解
sym_matrix = Matrix(matrix)

# 计算零空间(齐次方程的解空间)
null_space = sym_matrix.nullspace()

if not null_space:
raise ValueError("无法配平方程式,可能输入有误")

# 取基础解系的第一个解
solution = null_space[0]

# 将解转换为分数
fractions = [Fraction(float(x)).limit_denominator() for x in solution]

# 找到公分母,将解转换为整数
denominators = [f.denominator for f in fractions]
least_common_multiple = lcm(denominators)

coefficients = [int(f * least_common_multiple) for f in fractions]

# 确保系数为正数
if any(c < 0 for c in coefficients):
coefficients = [-c for c in coefficients]

# 提取反应物和产物的系数
reactant_coeffs = coefficients[:len(reactants)]
product_coeffs = coefficients[len(reactants):]

return reactant_coeffs, product_coeffs

def main():
print("化学方程式配平器(矩阵法)")
print("请按照提示输入反应物和生成物,使用正确的化学符号(如H2O、CO2)")
print("注意:元素符号首字母大写,如'O'表示氧,'Fe'表示铁\n")

# 获取反应物
reactants_input = input("请输入反应物,用加号(+)分隔(例如:H2+O2): ").strip()
reactants_list = [r.strip() for r in reactants_input.split('+')]

# 获取生成物
products_input = input("请输入生成物,用加号(+)分隔(例如:H2O): ").strip()
products_list = [p.strip() for p in products_input.split('+')]

try:
# 配平方程式
reactant_coeffs, product_coeffs = balance_equation(reactants_list, products_list)

# 格式化输出
balanced_reactants = []
for coeff, reactant in zip(reactant_coeffs, reactants_list):
if coeff == 1:
balanced_reactants.append(reactant)
else:
balanced_reactants.append(f"{coeff}{reactant}")

balanced_products = []
for coeff, product in zip(product_coeffs, products_list):
if coeff == 1:
balanced_products.append(product)
else:
balanced_products.append(f"{coeff}{product}")

balanced_equation = " + ".join(balanced_reactants) + " = " + " + ".join(balanced_products)

print("\n配平结果:")
print(balanced_equation)

except Exception as e:
print(f"\n配平过程中出现错误:{str(e)}")
print("请检查输入的化学式是否正确(注意元素符号大小写和括号使用)")

if __name__ == "__main__":
main()

更新

8 月 22 日更新:
找到括号的那一步可以利用 rindex() 方法简化。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
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
+ j = s.rindex(")") + 1
# 递归处理括号内
...