mathematica 机器证明平面几何定理【例 1】

TSC999

李涛，张景中的高徒，专门研究机器证明。本人曾看过李涛多年前写的一篇论文。按此论文的思路，使用 mathematica 为平台，用机器证明的方法，试证了几个平面几何定理。
至于专用的机器证明软件，国内外应该都有，只是不容易下载得到。
本文的证明没有采用专用软件，通过阅读程序，也许可以有个机器证明的初步概念。

下面是编程证明平面几何中的泰博定理。

1. 
   
2. u =.; v =.; a =.;
   
3. XI = 0; YI = 0; XA = u; YA = -1; XB = v; YB = -1;
   
4. ZA = XA + YA I; ZB = XB + YB I; ZI = XI + YI I;
   
5. 
   
6. XC =.; YC =.; aa = {XC, YC} /.
   
7.   Solve[{(
   
8.      XA Sqrt[(XB - XC)^2 + (YB - YC)^2] +
   
9.       XB Sqrt[(XA - XC)^2 + (YA - YC)^2] +
   
10.       XC Sqrt[(XA - XB)^2 + (YA - YB)^2])/(
    
11.      Sqrt[(XA - XB)^2 + (YA - YB)^2] +
    
12.       Sqrt[(XB - XC)^2 + (YB - YC)^2] +
    
13.       Sqrt[(XC - XA)^2 + (YC - YA)^2]) == 0, (
    
14.      YA Sqrt[(XB - XC)^2 + (YB - YC)^2] +
    
15.       YB Sqrt[(XA - XC)^2 + (YA - YC)^2] +
    
16.       YC Sqrt[(XA - XB)^2 + (YA - YB)^2])/(
    
17.      Sqrt[(XA - XB)^2 + (YA - YB)^2] +
    
18.       Sqrt[(XB - XC)^2 + (YB - YC)^2] +
    
19.       Sqrt[(XC - XA)^2 + (YC - YA)^2]) == 0,
    
20.     Abs[XA (YB - YC) + XB (YC - YA) + XC (YA - YB)]/(
    
21.      Sqrt[(XA - XB)^2 + (YA - YB)^2] +
    
22.       Sqrt[(XB - XC)^2 + (YB - YC)^2] +
    
23.       Sqrt[(XC - XA)^2 + (YC - YA)^2]) == 1}, {XC, YC }];
    
24. bb = aa[[1]];
    
25. XC = bb[[1]];
    
26. YC = bb[[2]];
    
27. ZC = Simplify[XC + YC I];
    
28. Print["C = ", ZC]
    
29. 
    
30. (* 下面求三角形 ABC 外接圆圆心 O 的坐标及其半径 R  *)
    
31. XO = ((XA^2 (YB - YC) + XB^2 ( YC - YA) +
    
32.      XC^2 ( YA - YB)) - (YA - YB ) ( YB - YC) (YC - YA))/(
    
33.   2 XA (YB - YC) + 2 XB ( YC - YA) + 2 XC ( YA - YB));
    
34. YO = -(((YA^2 (XB - XC) + YB^2 ( XC - XA) +
    
35.       YC^2 ( XA - XB)) - (XA - XB ) ( XB - XC) (XC - XA))/(
    
36.    2 XA (YB - YC) + 2 XB ( YC - YA) + 2 XC ( YA - YB)));
    
37. ZO = Factor[XO + YO I];
    
38. Print["O = ", ZO]
    
39. (* 下面求三角形 ABC 外接圆的半径 R  *)
    
40. R = Assuming[u < 0 && v > 0 && 1 + u v < 0,
    
41.    Factor[Simplify[
    
42.      Sqrt[((XA - XB)^2 + (YA - YB)^2)  ((XB - XC)^2 + (YB -
    
43.           YC)^2) ((XC - XA)^2 + (YC - YA)^2)]/(
    
44.      2 (XA (YB - YC) + XB (YC - YA) + XC (YA - YB)))]]];
    
45. Print["R = ", R];
    
46. (* 下面求 OF 线与 \[CircleDot]O 的交点 E 的坐标  *)
    
47. XE = Factor[ComplexExpand[Re[ZO]]];
    
48. YE = Factor[ComplexExpand[Im[ZO]] - R];
    
49. ZE = Simplify[XE + YE I];
    
50. Print["E = ", ZE]
    
51. ZF = (ZA + ZB)/2; (* F 是 AB 的中点  *)
    
52. ZK = Simplify[(1 - a) ZA +
    
53.     a ZF] ;  (* 在 AB 线上求一个 K 点，其位置由 a 决定。a 是介于 0 与 1 的实数。  *)
    
54. Print["K = ", ZK]
    
55. XK = Factor[ComplexExpand[Re[ZK]]];
    
56. YK = -1;
    
57. (* 下面求 EK 的延长线与 \[CircleDot]O 的交点 L 的坐标  *)
    
58. x =.; y =.; aa = {{x, y}} /.
    
59.   Simplify[Assuming[u < 0 && v > 0 && 1 + u v < 0 && x < XK,
    
60.     Solve[{(x - XK)/(y - YK) == (XK - XE)/(
    
61.        YK - YE), (x - XO)^2 + (y - YO)^2 == R^2}, {Factor[
    
62.        Simplify[x]], Factor[Simplify[y]] }]]];
    
63. bb = aa[[2]];(* 交点有两个，第一组解不符合要求，舍去。取第二组解  *)
    
64. cc = bb[[1]];
    
65. XL = cc[[1]];
    
66. YL = cc[[2]];
    
67. ZL = Simplify[XL + YL I]   ;
    
68. Print["L = ", ZL]
    
69. 
    
70. (* 下面求 \[CircleDot]I1 圆心 I1 的复坐标，I1 点是 OL 线与通过 K 点的 AB 的垂线的交点  *)
    
71. (* x=XK 是通过 K 点的 AB 的垂线 KK1 的直线方程，KK1 与 OL 的交点即是 I1 点  *)
    
72. XL = Factor[ComplexExpand[Re[ZL]]];
    
73. YL = Factor[ComplexExpand[Im[ZL]]];
    
74. x =.; y =.; aa = {{x, y}} /.
    
75.   Simplify[Solve[{(x - XO)/(y - YO) == (XO - XL)/(YO - YL),
    
76.      x == XK}, {Factor[Simplify[x]], Factor[Simplify[y]] }]];
    
77. bb = aa[[1]];(* 交点只有一个，取第一组解  *)
    
78. cc = bb[[1]];
    
79. XI1 = cc[[1]];
    
80. YI1 = cc[[2]];
    
81. ZI1 = Simplify[Factor[XI1 + YI1 I] ];
    
82. Print["I1 = ", ZI1]
    
83. 
    
84. (* 从 C 点向 \[CircleDot]I1 作切线，切点为 M, 下面求 M 点的坐标  *)
    
85. 
    
86. x1 = XI1; y1 = YI1;(* 圆心坐标 *)
    
87. x2 = XC; y2 = YC;     (* 向圆引切线的那个定点的坐标 *)
    
88. r = YI1 - YK; (* 左边 \[CircleDot]I1 圆的半径 *)
    
89. 
    
90. (*  所求的切点是 M1(XM1,YM1) 和 M2(XM2,YM2)，已知有以下公式：  *)
    
91. XM1 = x1 - (
    
92.    r^2 (x1 - x2) - (y1 - y2) r Sqrt[(x1 - x2)^2 + (y1 - y2)^2 -
    
93.       r^2])/((x1 - x2)^2 + (y1 - y2)^2) ;
    
94. YM1 = y1 - (
    
95.    r^2 (y1 - y2) + (x1 - x2) r Sqrt[(x1 - x2)^2 + (y1 - y2)^2 -
    
96.       r^2])/((x1 - x2)^2 + (y1 - y2)^2) ;
    
97.           (* 上面这个 ZM1 切点不符合要求，舍去 *)
    
98. XM2 = x1 - (
    
99.    r^2 (x1 - x2) + (y1 - y2) r Sqrt[(x1 - x2)^2 + (y1 - y2)^2 -
    
100.       r^2])/((x1 - x2)^2 + (y1 - y2)^2);
     
101. YM2 = y1 - (
     
102.    r^2 (y1 - y2) - (x1 - x2) r Sqrt[(x1 - x2)^2 + (y1 - y2)^2 -
     
103.       r^2])/((x1 - x2)^2 + (y1 - y2)^2);
     
104. 
     
105. ZM = Assuming[
     
106.    1 + u v < 0 &&
     
107.     1 > a > 0 && ((2 ((1 + u^2)) v + a ((u - v - u^2 v + u v^2)))) >
     
108.      0, Simplify[Factor[XM2 + YM2 I]]];
     
109. (*  上面这条语句中，((2((1+u^2))v+a((u-v-u^2v+u v^2))))>0 \
     
110. 这个条件是程序调试中发现的，不可能事先分析出来。若不加上这个条件，表达式就不能充分化简  *)
     
111. Print["M = ", ZM]
     
112. 
     
113. (* 将 CM 延长交 AB 于 N, 下面求 N 点的坐标  *)
     
114. XM = Factor[ComplexExpand[Re[ZM]]];
     
115. YM = Factor[ComplexExpand[Im[ZM]]];
     
116. 
     
117. x =.; y =.; aa = {{x, y}} /.
     
118.   Simplify[Solve[{(x - XC)/(y - YC) == (XC - XM)/(YC - YM),
     
119.      y == -1}, {Factor[Simplify[x]], Factor[Simplify[y]] }]];
     
120. bb = aa[[1]];(* 交点只有一个，取第一组解  *)
     
121. cc = bb[[1]];
     
122. XN = cc[[1]];
     
123. YN = cc[[2]];
     
124. ZN = Simplify[Factor[XN + YN I] ];
     
125. Print["N = ", ZN]
     
126. 
     
127. (* 延长 NB 到 N1, 使 NN1=NC  *)
     
128. NC = Assuming[
     
129.    u < 0 && v > 0 && 1 + u v < 0 && 1 > a > 0 && (-2 + a) u - a v < 0,
     
130.     Simplify[Factor[Sqrt[(XN - XC)^2 + (YN - YC)^2]]]];
     
131. XN1 = Simplify[XN + NC];
     
132. ZN1 = Simplify[Factor[XN1 - 1 I] ];
     
133. Print["N1 = ", ZN1]
     
134. (* 求出 N1C 的中点 N2  *)
     
135. ZN2 = Simplify[Factor[(ZC + ZN1)/2]];
     
136. Print["N2 = ", ZN2]
     
137. 
     
138. (* 求出 NN2 与 I1I 的交点 I2, 此即 \[CircleDot]I2 的圆心 *)
     
139. XN2 = Factor[ComplexExpand[Re[ZN2]]];
     
140. YN2 = Factor[ComplexExpand[Im[ZN2]]];
     
141. x =.; y =.; aa = {{x, y}} /.
     
142.   Simplify[Solve[{(x - XN2)/(y - YN2) == (XN2 - XN)/(YN2 - YN), (
     
143.       x - XI)/(y - YI) == (XI - XI1)/(YI - YI1)}, {Factor[
     
144.       Simplify[x]], Factor[Simplify[y]] }]];
     
145. bb = aa[[1]];(* 交点只有一个，取第一组解  *)
     
146. cc = bb[[1]];
     
147. XI2 = cc[[1]];
     
148. YI2 = cc[[2]];
     
149. ZI2 = Simplify[Factor[XI2 + YI2 I]];
     
150. Print["I2 = ", ZI2];
     
151. (* 求出 \[CircleDot]I2 圆的半径, 为此从 I2 点向 AB 作垂线, 垂足为 P, I2P 即是半径 *)
     
152. I2P = YI2 + 1;
     
153. R2 = I2P;
     
154. XP = XI2;
     
155. YP = -1;
     
156. ZP = Simplify[Factor[XP + YP I]];
     
157. Print["P = ", ZP];
     
158. 
     
159. (* 证明 \[CircleDot]I2 圆与 \[CircleDot]O 圆相切，就完成了全部证明。 *)
     
160. S = Assuming[
     
161.    u < 0 && v > 0 && 1 + u v < 0 &&  
     
162.     1 - 2 (-2 + a) a u v + (-1 + a)^2 v^2 +
     
163.       u^2 (1 - 2 a + a^2 + v^2) > 0,
     
164.    Simplify[Sqrt[(XI2 - XO)^2 + (YI2 - YO)^2]]];
     
165. Print["S+R2 = ", Factor[Simplify[S + R2]]];
     
166. Print["R = ", R];
     
167. If[Factor[Simplify[S + R2]] == R, Print["由于 \!\(\*
     
168. StyleBox["S",\nFontColor->GrayLevel[0]]\)\!\(\*
     
169. StyleBox["+",\nFontColor->GrayLevel[0]]\)\!\(\*
     
170. StyleBox["R2",\nFontColor->GrayLevel[0]]\)\!\(\*
     
171. StyleBox["=",\nFontColor->GrayLevel[0]]\)R，因此 \[CircleDot]I2 圆与 \
     
172. \[CircleDot]O 圆相切。"]]
     
173. 
     

复制代码

程序运行结果：



上述运行结果，与李涛论文所给结果完全相同。李涛的结果是用他专门开发的机器证明软件做的，如果写出程序，必然是有许多子程序调用，因而程序看起来会简明得多。

TSC999

我觉得机器证明中的一个最困难的环节，就是如何进行构图，使得每个关键点在复平面上的位置都能够用不太复杂的复数代数式来表达。如果公式过于复杂，mathematica 就不能胜任繁重的计算任务，即使经过长时间运行，也不能给出证明结果；或者即使能给出结果，但是需要时间过长，那就不是一个好的程序。