找回密码
 欢迎注册
查看: 20953|回复: 1

[分享] mathematica 机器证明平面几何定理【例 1】

[复制链接]
发表于 2017-5-23 11:07:27 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?欢迎注册

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

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

泰博定理图.png

  1. u =.; v =.; a =.;
  2. XI = 0; YI = 0; XA = u; YA = -1; XB = v; YB = -1;
  3. ZA = XA + YA I; ZB = XB + YB I; ZI = XI + YI I;

  4. XC =.; YC =.; aa = {XC, YC} /.
  5.   Solve[{(
  6.      XA Sqrt[(XB - XC)^2 + (YB - YC)^2] +
  7.       XB Sqrt[(XA - XC)^2 + (YA - YC)^2] +
  8.       XC Sqrt[(XA - XB)^2 + (YA - YB)^2])/(
  9.      Sqrt[(XA - XB)^2 + (YA - YB)^2] +
  10.       Sqrt[(XB - XC)^2 + (YB - YC)^2] +
  11.       Sqrt[(XC - XA)^2 + (YC - YA)^2]) == 0, (
  12.      YA Sqrt[(XB - XC)^2 + (YB - YC)^2] +
  13.       YB Sqrt[(XA - XC)^2 + (YA - YC)^2] +
  14.       YC Sqrt[(XA - XB)^2 + (YA - YB)^2])/(
  15.      Sqrt[(XA - XB)^2 + (YA - YB)^2] +
  16.       Sqrt[(XB - XC)^2 + (YB - YC)^2] +
  17.       Sqrt[(XC - XA)^2 + (YC - YA)^2]) == 0,
  18.     Abs[XA (YB - YC) + XB (YC - YA) + XC (YA - YB)]/(
  19.      Sqrt[(XA - XB)^2 + (YA - YB)^2] +
  20.       Sqrt[(XB - XC)^2 + (YB - YC)^2] +
  21.       Sqrt[(XC - XA)^2 + (YC - YA)^2]) == 1}, {XC, YC }];
  22. bb = aa[[1]];
  23. XC = bb[[1]];
  24. YC = bb[[2]];
  25. ZC = Simplify[XC + YC I];
  26. Print["C = ", ZC]

  27. (* 下面求三角形 ABC 外接圆圆心 O 的坐标及其半径 R  *)
  28. XO = ((XA^2 (YB - YC) + XB^2 ( YC - YA) +
  29.      XC^2 ( YA - YB)) - (YA - YB ) ( YB - YC) (YC - YA))/(
  30.   2 XA (YB - YC) + 2 XB ( YC - YA) + 2 XC ( YA - YB));
  31. YO = -(((YA^2 (XB - XC) + YB^2 ( XC - XA) +
  32.       YC^2 ( XA - XB)) - (XA - XB ) ( XB - XC) (XC - XA))/(
  33.    2 XA (YB - YC) + 2 XB ( YC - YA) + 2 XC ( YA - YB)));
  34. ZO = Factor[XO + YO I];
  35. Print["O = ", ZO]
  36. (* 下面求三角形 ABC 外接圆的半径 R  *)
  37. R = Assuming[u < 0 && v > 0 && 1 + u v < 0,
  38.    Factor[Simplify[
  39.      Sqrt[((XA - XB)^2 + (YA - YB)^2)  ((XB - XC)^2 + (YB -
  40.           YC)^2) ((XC - XA)^2 + (YC - YA)^2)]/(
  41.      2 (XA (YB - YC) + XB (YC - YA) + XC (YA - YB)))]]];
  42. Print["R = ", R];
  43. (* 下面求 OF 线与 \[CircleDot]O 的交点 E 的坐标  *)
  44. XE = Factor[ComplexExpand[Re[ZO]]];
  45. YE = Factor[ComplexExpand[Im[ZO]] - R];
  46. ZE = Simplify[XE + YE I];
  47. Print["E = ", ZE]
  48. ZF = (ZA + ZB)/2; (* F 是 AB 的中点  *)
  49. ZK = Simplify[(1 - a) ZA +
  50.     a ZF] ;  (* 在 AB 线上求一个 K 点,其位置由 a 决定。a 是介于 0 与 1 的实数。  *)
  51. Print["K = ", ZK]
  52. XK = Factor[ComplexExpand[Re[ZK]]];
  53. YK = -1;
  54. (* 下面求 EK 的延长线与 \[CircleDot]O 的交点 L 的坐标  *)
  55. x =.; y =.; aa = {{x, y}} /.
  56.   Simplify[Assuming[u < 0 && v > 0 && 1 + u v < 0 && x < XK,
  57.     Solve[{(x - XK)/(y - YK) == (XK - XE)/(
  58.        YK - YE), (x - XO)^2 + (y - YO)^2 == R^2}, {Factor[
  59.        Simplify[x]], Factor[Simplify[y]] }]]];
  60. bb = aa[[2]];(* 交点有两个,第一组解不符合要求,舍去。取第二组解  *)
  61. cc = bb[[1]];
  62. XL = cc[[1]];
  63. YL = cc[[2]];
  64. ZL = Simplify[XL + YL I]   ;
  65. Print["L = ", ZL]

  66. (* 下面求 \[CircleDot]I1 圆心 I1 的复坐标,I1 点是 OL 线与通过 K 点的 AB 的垂线的交点  *)
  67. (* x=XK 是通过 K 点的 AB 的垂线 KK1 的直线方程,KK1 与 OL 的交点即是 I1 点  *)
  68. XL = Factor[ComplexExpand[Re[ZL]]];
  69. YL = Factor[ComplexExpand[Im[ZL]]];
  70. x =.; y =.; aa = {{x, y}} /.
  71.   Simplify[Solve[{(x - XO)/(y - YO) == (XO - XL)/(YO - YL),
  72.      x == XK}, {Factor[Simplify[x]], Factor[Simplify[y]] }]];
  73. bb = aa[[1]];(* 交点只有一个,取第一组解  *)
  74. cc = bb[[1]];
  75. XI1 = cc[[1]];
  76. YI1 = cc[[2]];
  77. ZI1 = Simplify[Factor[XI1 + YI1 I] ];
  78. Print["I1 = ", ZI1]

  79. (* 从 C 点向 \[CircleDot]I1 作切线,切点为 M, 下面求 M 点的坐标  *)

  80. x1 = XI1; y1 = YI1;(* 圆心坐标 *)
  81. x2 = XC; y2 = YC;     (* 向圆引切线的那个定点的坐标 *)
  82. r = YI1 - YK; (* 左边 \[CircleDot]I1 圆的半径 *)

  83. (*  所求的切点是 M1(XM1,YM1) 和 M2(XM2,YM2),已知有以下公式:  *)
  84. XM1 = x1 - (
  85.    r^2 (x1 - x2) - (y1 - y2) r Sqrt[(x1 - x2)^2 + (y1 - y2)^2 -
  86.       r^2])/((x1 - x2)^2 + (y1 - y2)^2) ;
  87. YM1 = y1 - (
  88.    r^2 (y1 - y2) + (x1 - x2) r Sqrt[(x1 - x2)^2 + (y1 - y2)^2 -
  89.       r^2])/((x1 - x2)^2 + (y1 - y2)^2) ;
  90.           (* 上面这个 ZM1 切点不符合要求,舍去 *)
  91. XM2 = 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. YM2 = 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. ZM = Assuming[
  98.    1 + u v < 0 &&
  99.     1 > a > 0 && ((2 ((1 + u^2)) v + a ((u - v - u^2 v + u v^2)))) >
  100.      0, Simplify[Factor[XM2 + YM2 I]]];
  101. (*  上面这条语句中,((2((1+u^2))v+a((u-v-u^2v+u v^2))))>0 \
  102. 这个条件是程序调试中发现的,不可能事先分析出来。若不加上这个条件,表达式就不能充分化简  *)
  103. Print["M = ", ZM]

  104. (* 将 CM 延长交 AB 于 N, 下面求 N 点的坐标  *)
  105. XM = Factor[ComplexExpand[Re[ZM]]];
  106. YM = Factor[ComplexExpand[Im[ZM]]];

  107. x =.; y =.; aa = {{x, y}} /.
  108.   Simplify[Solve[{(x - XC)/(y - YC) == (XC - XM)/(YC - YM),
  109.      y == -1}, {Factor[Simplify[x]], Factor[Simplify[y]] }]];
  110. bb = aa[[1]];(* 交点只有一个,取第一组解  *)
  111. cc = bb[[1]];
  112. XN = cc[[1]];
  113. YN = cc[[2]];
  114. ZN = Simplify[Factor[XN + YN I] ];
  115. Print["N = ", ZN]

  116. (* 延长 NB 到 N1, 使 NN1=NC  *)
  117. NC = Assuming[
  118.    u < 0 && v > 0 && 1 + u v < 0 && 1 > a > 0 && (-2 + a) u - a v < 0,
  119.     Simplify[Factor[Sqrt[(XN - XC)^2 + (YN - YC)^2]]]];
  120. XN1 = Simplify[XN + NC];
  121. ZN1 = Simplify[Factor[XN1 - 1 I] ];
  122. Print["N1 = ", ZN1]
  123. (* 求出 N1C 的中点 N2  *)
  124. ZN2 = Simplify[Factor[(ZC + ZN1)/2]];
  125. Print["N2 = ", ZN2]

  126. (* 求出 NN2 与 I1I 的交点 I2, 此即 \[CircleDot]I2 的圆心 *)
  127. XN2 = Factor[ComplexExpand[Re[ZN2]]];
  128. YN2 = Factor[ComplexExpand[Im[ZN2]]];
  129. x =.; y =.; aa = {{x, y}} /.
  130.   Simplify[Solve[{(x - XN2)/(y - YN2) == (XN2 - XN)/(YN2 - YN), (
  131.       x - XI)/(y - YI) == (XI - XI1)/(YI - YI1)}, {Factor[
  132.       Simplify[x]], Factor[Simplify[y]] }]];
  133. bb = aa[[1]];(* 交点只有一个,取第一组解  *)
  134. cc = bb[[1]];
  135. XI2 = cc[[1]];
  136. YI2 = cc[[2]];
  137. ZI2 = Simplify[Factor[XI2 + YI2 I]];
  138. Print["I2 = ", ZI2];
  139. (* 求出 \[CircleDot]I2 圆的半径, 为此从 I2 点向 AB 作垂线, 垂足为 P, I2P 即是半径 *)
  140. I2P = YI2 + 1;
  141. R2 = I2P;
  142. XP = XI2;
  143. YP = -1;
  144. ZP = Simplify[Factor[XP + YP I]];
  145. Print["P = ", ZP];

  146. (* 证明 \[CircleDot]I2 圆与 \[CircleDot]O 圆相切,就完成了全部证明。 *)
  147. S = Assuming[
  148.    u < 0 && v > 0 && 1 + u v < 0 &&  
  149.     1 - 2 (-2 + a) a u v + (-1 + a)^2 v^2 +
  150.       u^2 (1 - 2 a + a^2 + v^2) > 0,
  151.    Simplify[Sqrt[(XI2 - XO)^2 + (YI2 - YO)^2]]];
  152. Print["S+R2 = ", Factor[Simplify[S + R2]]];
  153. Print["R = ", R];
  154. If[Factor[Simplify[S + R2]] == R, Print["由于 \!\(\*
  155. StyleBox["S",\nFontColor->GrayLevel[0]]\)\!\(\*
  156. StyleBox["+",\nFontColor->GrayLevel[0]]\)\!\(\*
  157. StyleBox["R2",\nFontColor->GrayLevel[0]]\)\!\(\*
  158. StyleBox["=",\nFontColor->GrayLevel[0]]\)R,因此 \[CircleDot]I2 圆与 \
  159. \[CircleDot]O 圆相切。"]]

复制代码


程序运行结果:

泰博定理图.png

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



毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2017-5-23 11:24:47 | 显示全部楼层
我觉得机器证明中的一个最困难的环节,就是如何进行构图,使得每个关键点在复平面上的位置都能够用不太复杂的复数代数式来表达。如果公式过于复杂,mathematica 就不能胜任繁重的计算任务,即使经过长时间运行,也不能给出证明结果;或者即使能给出结果,但是需要时间过长,那就不是一个好的程序。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|数学研发网 ( 苏ICP备07505100号 )

GMT+8, 2024-12-23 00:41 , Processed in 0.026151 second(s), 19 queries .

Powered by Discuz! X3.5

© 2001-2024 Discuz! Team.

快速回复 返回顶部 返回列表