数学研发论坛

 找回密码
 欢迎注册
123
返回列表 发新帖
楼主: zhouguang

[提问] 立方体中的最大内接正方形

[复制链接]
发表于 2018-10-25 18:59:32 | 显示全部楼层
楼上矩阵最优解要求各行元素平方和相等,各列元素绝对值和为1,并且任意不同两行之间正交。可以看出,对于任意这样的矩阵,交换任意两行两列或改变任意行列所有数据的符号,还是满足条件的。利用这些性质,至少对于3*4阶矩阵,我们可以穷举所有情况,其中非零数不超过9个,反复利用上面性质变换,最后可以得出总共只会余下四种:
下面所有字母代表[0,1]区间中实数:
\(\begin{pmatrix}-a & d& 1-g& 1\\b& -e&g & 0\\1-a-b&1-d-e&0&0\end{pmatrix}\)
对应方程:
\(\begin{cases}a(1-a-b)=d(1-d-e)\\b(1-a-b)=e(1-d-e)\\ab+de=(1-g)g\\a^2+d^2+(1-g)^2+1=b^2+e^2+g^2=(1-a-b)^2+(1-d-e)^2\end{cases}\)
这个方程的关键是前两式相加,得a+b=c+d或a+b+c+d=1,其中第二中解有负数出现,而第一种情况得b是方程$16*b^4 + 8*b^3 - 23*b^2 + 14*b - 2$的根,选择$e=b=0.20490155350665129314330190383234265145$得$a=d=0.082734498297453867827044842101191253555,g=0.96486035057909231410839945902474434917$
对应边长的平方为1.0149247893784870917727027004682468871


\(\begin{pmatrix} -a&d&g&0\\ b&-e&1-g&0\\1-a-b&1-d-e&0&1\end{pmatrix}\)
对应方程
\(\begin{cases}a(1-a-b)=d(1-d-e)\\b(1-a-b)=e(1-d-e)\\ab+de=g(1-g)\\a^2+d^2+g^2=b^2+e^2+(1-g)^2=(1-a-b)^2+(1-d-e)^2+1\end{cases}\)
这个方程解法和上面方程类似,没有得出更好的解

\(\begin{pmatrix}-a&d&1-f&0\\b&1-d&0&h\\c&0&f&h-1\end{pmatrix}\)
对应方程
\(\begin{cases}ab=d(1-d)\\bc=h(1-h)\\ac=f(1-f)\\a^2+d^2+(1-f)^2=b^2+(1-d)^2+h^2=c^2+f^2+(1-h)^2\\a+b+c=1\end{cases}\)
这个解出来结果比较可怕,比如h需要满足方程:
$(h-1)*h*(h^2-h+1)*(2h^2-2h-1)(9h^2-9h+1)(96*h^4 - 176*h^3 + 92*h^2 + 6*h - 9)
   (4*h^6 - 12*h^5 + 16*h^4 - 12*h^3 + 6*h^2 - 2*h + 1)
  (1458*h^12 - 8748*h^11 + 18711*h^10 - 13365*h^9 - 13797*h^8 + 39150*h^7 - 29122*h^6 - 9564*h^5 + 25793*h^4 - 11355*h^3 + 555*h^2 + 284*h + 8)
  (4608000*h^14 - 32486400*h^13 + 106297280*h^12 - 208927040*h^11 + 273352416*h^10 - 252211768*h^9 + 175199545*h^8 - 104473284*h^7 + 54392414*h^6 - 14730172*h^5 - 1523377*h^4 + 436852*h^3 + 334918*h^2 + 46456*h + 3721)
  (165888*h^16 - 1327104*h^15 + 4783104*h^14 - 10167552*h^13 + 13055616*h^12 - 7273152*h^11 - 6556608*h^10 + 17191200*h^9 - 12065016*h^8 - 5695176*h^7 + 16083319*h^6 - 10125760*h^5 + 1059958*h^4 + 1398650*h^3 - 611586*h^2 + 150880*h + 51253)
(5308416*h^26 - 69009408*h^25 + 525975552*h^24 - 21776007168*h^23 + 279629595648*h^22 - 1755045642240*h^21 + 6774140980224*h^20 - 17688419204736*h^19 + 32024536790592*h^18 - 38324963108544*h^17 + 23050469214944*h^16 + 11097784110800*h^15 - 38178355003336*h^14 + 32961273735048*h^13 - 2566180746576*h^12 - 19901614525832*h^11 + 16553725604058*h^10 - 2217221289542*h^9 - 3810021021401*h^8 + 1602353883447*h^7 + 439842772479*h^6 - 325579847011*h^5 - 32978746973*h^4 + 41942923279*h^3 - 3484899190*h^2 - 547935318*h + 357216)
(7644119040000*h^39 - 199860269875200*h^38 + 2437351880785920*h^37 - 18759902729011200*h^36 + 103074607169667072*h^35 - 435377615689838592*h^34 + 1491209724313731072*h^33 - 4298851157566006272*h^32 + 10644024801160240128*h^31 - 22725803111402993472*h^30 + 41473159626822259712*h^29 - 63493882655040338784*h^28 + 78994789689497720864*h^27 - 74719500922192891184*h^26 + 42940885445432024752*h^25 + 9534775436809759632*h^24 - 61502974905383268024*h^23 + 89004913392991538292*h^22 - 80304104536362644692*h^21 + 43080843348906315238*h^20 + 836842045447274914*h^19 - 29487275391337400990*h^18 + 34209282358889258953*h^17 - 22554620937493320006*h^16 + 8738926791071167945*h^15 - 1443514589525468923*h^14 + 251119191467260523*h^13 - 920986793774588347*h^12 + 995856852080848484*h^11 - 531495875439132419*h^10 + 157583299487498045*h^9 - 23090925236014258*h^8 + 443993273326087*h^7 + 309966779646301*h^6 - 44102946274893*h^5 - 5733815275467*h^4 + 3295423854898*h^3 + 188330030823*h^2 + 5250435188*h + 110766728)$
最终只有合法解
$9h^2-9h+1=0=>d=f=h => b=1/3,a=1/3$也不是最优解


\(\begin{pmatrix}0&f&f&0\\1-f&1-f&f-1&1-f\\ -f&0&0&f\end{pmatrix}\)
对应解满足$f^2-4f+2=0,f=0.58578643762690495119831127579030192143$边长平方为$2f^2=0.68629150101523960958649020632241537145$不是最优解
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2018-11-3 11:51:18 | 显示全部楼层
根据20#的公式,我们知道对于m阶情况,我们需要构造$m\times n$阶矩阵,每行范数相等,任意两行正交。而每列绝对值和不超过1。
但是对于最优解,必然每列绝对值和等于1,而且矩阵中所有非零数字数目不超过${(m+2)(m-1)}/2+n$个。
首先由于每列绝对值和不超过1,所以每列数字平方和不超过1,所以$m\times n$个数字平方和不超过$n$,于是每行元素平方和不超过$n/m$,
我们得出最大内接m维超立方体棱长不超过$\sqrt(n/m)$,其中只有n正好是m的倍数的时候才能够取到等号,这时要求每列正好一个数绝对值是1。
而另外一方面,我们总是可以选择$m[n/m]$列,分成$[n/m]$组,每组$m$列取成一个$m\times m$单位阵,余下$n-m[n/m]$列全部为0,得到一个非最优解可以取到棱长为$\sqrt([n/m])$的m维超立方体。所以最优解棱长必然大于$\sqrt([n/m])$而不超过$\sqrt(n/m)$。
而对于$m=3$情况,我们应该可以计算出最优解
我们假设某个最优解有k列不止一个元素绝对值大于0,那么对于$m=3,k<=5$
那么对于这k列构成$3\times k$阶的子矩阵,它的任意两行还是正交,但是任意两行数字平方和的差值是整数。
而$k=0$对应n是3的倍数时的最优值,已经讨论过了,所以只要讨论$1<=k<=5$的情况。
第一种情况:
这$3\times k$阶的子矩阵有一行全部是0,那么余下两行每列都是两个绝对值和为1的数,相当于只有k个变量,额外约束是这两行数平方和是整数,而且这两行正交,
  k个变量只有3个约束。由于全局最优解必然取边界条件,对于$k>3$时必然还要取更多边界条件直到$k=3$
  而对于$k=1$和$k=2$,容易看出这种情况无解。所以这种情况只能$k=3$
于是对应矩阵\(\begin{pmatrix}a&b&c\\1-a&1-b&c-1\\0&0&0\end{pmatrix}\)
其中由于$3/4<=a^2+b^2+c^2+(1-a)^2+(1-b)^2+(c-1)^2<3$,而且$a^2+b^2+c^2>=1,(1-a)^2+(1-b)^2+(c-1)^2>=1$,而且它们都是整数,所以只能
\(\begin{cases}a^2+b^2+c^2=1\\(1-a)^2+(1-b)^2+(c-1)^2=1\\a(1-a)+b(1-b)=c(1-c)\end{cases}\)
前两式相减得出$a+b+c=3/2$,得出a,b满足方程$8x^2-8x+1=0$,即$a=0.14644660940672623779957781894757548036,b=0.85355339059327376220042218105242451964,c=1/2$
对应$n=3h+4$时,可以得到棱长平方为$h+1={n-1}/3$
第二种情况:
这$3\times k$阶的子矩阵三行都有非零数。其中如果某一行只有一个非零数,由于这一行和其它行正交,我们可以得出这个非零数所在列其它数都是0,与每列至少两个非零数矛盾。
所以每行至少两个非零数。另外每两行正交,所以三行线性无关,得出$k>=3$.
2.1 如果k=5, 由于非零数比列出最多多5,所以只能每列都正好两个数。
   其中如果一行5个非零数,那么另外两行就分别3个和2个非零数,对应矩阵
\(\begin{pmatrix}a&b&c&d&e\\1-a&1-b&c-1&0&0\\0&0&0&1-d&e-1\end{pmatrix}\)
由于5个变量其中只有4个约束条件(每两行平方和之差是整数两个条件,两行之间两两正交还是两个条件),所以还得有个变量取边界值,所以退化到$k<5$的情况
  另外两种模式4+4+2和4+3+3都容易利用正交性得出无解。
比如4+4+2,说明2个非零元那行都和另两行都只有一列同时有非零数,和它们正交矛盾。4+3+3得出3个非零元的两行只有一列同时有非零数,和它们正交矛盾。
所以k=5淘汰
2.2 如果k=4
  那么最多9个非零数(不超过k+5),最少8个非零数(不少于2k), 对应三行包含的非零数数目有
    4+3+2
     对应矩阵\(\begin{pmatrix}a&b&c&d\\1-a&1-b&e&0\\0&0&1-c-e&d-1\end{pmatrix}\),和第二第三行正交矛盾,淘汰
    4+2+2
       对应矩阵\(\begin{pmatrix}a&b&c&d\\1-a&b-1&0&0\\0&0&1-c&d-1\end{pmatrix}\)
    3+3+3
      对应矩阵\(\begin{pmatrix}a&b&c&0\\d&b-1&0&e\\1-a-d&0&c-1&e-1\end{pmatrix}\)
    3+3+2
      对应矩阵\(\begin{pmatrix}a&b&c&0\\1-a&b-1&0&e\\0&0&c-1&d-1\end{pmatrix}\),其中第三行和前两行都无法正交,淘汰。
  四种模式。
2.3 如果k=3
   那么最多8个非零数(不超过k+5),最少6个非零数(不小于2k),对应三行包含的非零数数目有
   3+3+2
     对应矩阵\(\begin{pmatrix}a&b&c\\d&e&c-1\\1-a-d&b+e-1&0\end{pmatrix}\)
   3+2+2
     对应矩阵\(\begin{pmatrix}a&b&c\\d&0&c-1\\1-a-d&b-1&0\end{pmatrix}\),其中第二第三行无法正交,淘汰
   2+2+2
     对应矩阵\(\begin{pmatrix}a&b&0\\1-a&0&d\\0&b-1&1-d\end{pmatrix}\),其中各行之间无法正交,淘汰
三种模式

毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2018-11-3 15:23:12 | 显示全部楼层
2.2 4+2+2模式对应矩阵
\(\begin{pmatrix}a&b&c&d\\1-a&b-1&0&0\\0&0&1-c&d-1\end{pmatrix}\)
有$a(1-a)=b(1-b),c(1-c)=d(1-d)$,所以$a=b$或$a=1-b$,$c=d$或$c=1-d$
于是有
\(\begin{pmatrix}a&a&c&c\\1-a&a-1&0&0\\0&0&1-c&c-1\end{pmatrix}\)\(\begin{pmatrix}a&1-a&c&c\\1-a&-a&0&0\\0&0&1-c&c-1\end{pmatrix}\)\(\begin{pmatrix}a&1-a&c&1-c\\1-a&-a&0&0\\0&0&1-c&-c\end{pmatrix}\)
其中后面两行平方和相等或差1,同样它们和第一行的平方和也相差整数
第一个矩阵如果后面两行平方和相等,有$a=c$,所以分别得出
$a^2+2*a-1/2=0, a=0.22474487139158904909864203735294569598,4a^2=0.20204102886728760721086370117643443213,2(1-a)^2=1.20204102886728760721086370117643443213$
    对应n=3h+5, 立方体棱长平方$h+1.20204102886728760721086370117643443213=(n-2)/3+0.20204102886728760721086370117643443213$
$a^2+2a-1=0, a=0.41421356237309504880168872420969807857,4a^2=2(1-a)^2=0.68629150101523960958649020632241537145$
   对应n=3h+4,立方体棱长平方$h+0.68629150101523960958649020632241537145=(n-1)/3-0.31370849898476039041350979367758462856$不是最优情况 (n=1(mod 3))
$a^2+2a-3/2=0,a=0.58113883008418966599944677221635926686,4a^2=1.3508893593264826720044258222691258651,2(1-a)^2=0.3508893593264826720044258222691258651$
   对应n=3h+2+4,不是最优情况
$a^2 + 2*a - 2=0,a=0.73205080756887729352744634150587236694,4a^2=2.1435935394489816517804292679530210644,2(1-a)^2=0.1435935394489816517804292679530210644$
   对应n=3h+4+4, 立方体棱长平方为$h+2.1435935394489816517804292679530210644=(n-2)/3+0.1435935394489816517804292679530210644$,不如前面方案 (n=2(mod 3))
$a^2 + 2*a - 5/2=0,a=0.87082869338697069279187436615827465088,4a^2=3.0333704529042344576650050707338027930,2(1-a)62=0.0333704529042344576650050707338027930$
   对应n=3h+6+4, 立方体棱长平方为$h+3.0333704529042344576650050707338027930=(n-1)/3+0.0333704529042344576650050707338027930$ (n=1(mod 3))
第二个矩阵比较前两行平方和相差$2c^2$是整数,只能$2c^2=1$,得出无解,同样第三个矩阵比较前两行平方和要求$c^2+(1-c)^2$是整数同样不行
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2018-11-3 16:44:41 | 显示全部楼层
2.3的   3+3+2模式
     对应矩阵\(\begin{pmatrix}a&b&c\\d&e&c-1\\1-a-d&b+e-1&0\end{pmatrix}\)
对应方程组部分为
a(1-a-d)=b(1-b-e)
d(1-a-d)=e(1-b-e)
其中两式相加得到$a+d=b+e$或$a+d+b+e=1$
如果$a+d=b+e$,还可以得到$a=b,d=e$, 最后得到合法解
   $8*c^4 - 32*c^3 + 45*c^2 - 30*c + 1=0$
   c=0.035139649420907685891600540975255650831
   a=0.082734498297453867827044842101191253550,
   d=0.20490155350665129314330190383234265145,
   对应棱长平方为$(n-1)/3 + 0.014924789378487091772702700468246887087$ (n=1(mod 3))
$64*c^8 - 256*c^7 + 304*c^6 - 32*c^5 - 327*c^4 + 296*c^3 - 152*c^2 - 32*c + 16=0$
   c=0.27558461705654999784812328936037690179
   a=0.17987838803224237835524423692465721565
   d=0.55492418539619261753733901460781193790
  对应棱长平方为$(n-2)/3+0.14065935012036122417723850465213714931$ (n=2(mod 3))
  $16*c^8 - 64*c^7 + 64*c^6 + 24*c^5 - 96*c^4 + 24*c^3 + c^2 - 22*c + 3=0$
   c=0.13860489686232080596276176669196028601
   a=0.074380948371657708078488878468073906554
   d=0.80258172315533576161036783428400273187
$n=3h+3+4,  r^2=(n-1)/3 + 0.030276368395549006911606670186456882565$ (n=1(mod 3) 而且$n>=7$)
不过$a+d+b+e=1$情况的计算已经很复杂了,计算下来好像无解
而余下一个2.2的3+3+3模式矩阵对应情况会更加复杂 (如果这个方程解决了,那么m=3时就完全解决了)
除了这种情况,其余情况的结果
$n -=0(mod 3)$, 最大棱长平方 $n/3$
$n -=2(mod 3)$,  已知最大棱长平方${n-2}/3+0.20204102886728760721086370117643443213$
$n -=1(mod 3)$,  
  当$n<=4$时,已知最大棱长平方${n-1}/3+0.014924789378487091772702700468246887087$,
  而在$n=7$时,已知最大棱长平方为${n-1}/3+0.030276368395549006911606670186456882565$
  而在$n>=10$时,已知最大棱长平方${n-1}/3+0.0333704529042344576650050707338027930$
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2018-11-4 14:32:10 | 显示全部楼层
esol(0,0)
得到
9*c^2 - 9*c + 1
a=0.12732200375003505059847105521145396076;b=0.12732200375003505059847105521145396076;c=0.12732200375003505059847105521145396076;x=0.33333333333333333333333333333333333333;y=0.33333333333333333333333333333333333334;z=0.33333333333333333333333333333333333333
9*c^2 - 9*c + 1
a=0.87267799624996494940152894478879859054;b=0.87267799624996494940152894478796257612;c=0.87267799624996494940152894478854603924;x=0.33333333333333333333333333333202867014;y=0.33333333333333333333333333333431385034;z=0.33333333333333333333333333333365747953, 这个解不够好
$1327104*c^22 - 14598144*c^21 + 78520320*c^20 - 277198848*c^19 + 718562304*c^18 - 1445202432*c^17 + 2327525376*c^16 - 3056990976*c^15 + 3306544768*c^14 - 2954054848*c^13 + 2169816440*c^12 - 1288672208*c^11 + 592924988*c^10 - 185582420*c^9 + 14218336*c^8 + 27499858*c^7 - 22520955*c^6 + 10772277*c^5 - 3524433*c^4 + 663759*c^3 - 19311*c^2 - 11715*c + 1300$
$a=0.61645575611595019564366164451969012446;b=0.99964861889131459720628119890601018649;c=0.99914013369073816425105409912190416390;x=0.87543900054754931571057569585193373024;y=0.0028868720701670442017409435891338646630;z=0.12167412738228364008768336055893240510$
\(\begin{pmatrix}x&a&0&c-1\\y&a-1&b&0\\z&0&b-1&c\end{pmatrix}\)
    这组很好,对应$n -=1(mod 3)$, 棱长的平方$r^2=(n-1)/3+0.1464118822982497217046378348507596815$ 但是好像是增根?有条件检查好像通过不了,后面一个解也类似
11943936*c^30 - 214990848*c^29 + 1827090432*c^28 - 9753716736*c^27 + 36696038400*c^26 - 103401232128*c^25 + 226155106176*c^24 - 392322766464*c^23 + 547110063040*c^22 - 620405629664*c^21 + 584627278416*c^20 - 485442765424*c^19 + 400839008496*c^18 - 365668024160*c^17 + 346843699576*c^16 - 294398468288*c^15 + 201109515718*c^14 - 103025529162*c^13 + 35508759545*c^12 - 4751415874*c^11 - 3359452963*c^10 + 3163768959*c^9 - 1546309579*c^8 + 462530161*c^7 - 66873410*c^6 + 3478243*c^5 - 908424*c^4 - 214127*c^3 - 2115*c^2 - 2700*c - 200
a=0.71883970112472172592436784012526963004;b=0.27384239774228909248921926723362136821;c=0.68740721959507614742307169636665297282;x=0.028075493203415081208096501315243155710;y=0.67911134810728449916208563638479935889;z=0.29281315868930041962981786229995748540
增根

esol(0,-1)得到
1327104*c^24 - 15925248*c^23 + 91348992*c^22 - 335812608*c^21 + 889150464*c^20 - 1797152256*c^19 + 2860263936*c^18 - 3639095040*c^17 + 3706622592*c^16 - 2969565504*c^15 + 1755586808*c^14 - 577548952*c^13 - 195332444*c^12 + 488517392*c^11 - 449148924*c^10 + 281428234*c^9 - 126051203*c^8 + 33066423*c^7 + 4604354*c^6 - 11159397*c^5 + 6784002*c^4 - 2240179*c^3 + 367313*c^2 - 35036*c + 500
a=0.83075415866767342164081349566679015338;b=0.95646500394930519144717044967819810047;c=0.97650580758675700622895686025787971476;x=0.51496528479610405439835147108507665117;y=0.11146405329919368235620643373084308116;z=0.37357066190470226324544209518408026768
增根

esol(-1,-2)得到
1679616*c^28 - 20155392*c^27 + 109921536*c^26 - 357011712*c^25 + 746013888*c^24 - 952093440*c^23 + 386020224*c^22 + 1241707968*c^21 - 3242195168*c^20 + 4094308096*c^19 - 2678362336*c^18 - 536662928*c^17 + 3707328032*c^16 - 5104501336*c^15 + 4379168340*c^14 - 2497810860*c^13 + 717860116*c^12 + 277914444*c^11 - 507176717*c^10 + 333946194*c^9 - 106430119*c^8 - 11392092*c^7 + 26972507*c^6 - 9778522*c^5 + 103708*c^4 + 742752*c^3 - 100736*c^2 - 20824*c + 4600
a=0.022395238513895352908459101701502169208;b=0.16167701978121523069178099777660409280;c=0.91540024756583824365232716591630919668;x=0.11184581172946412415880434687636092276;y=0.19574887487748029835324145420547932654;z=0.69240531339305557748795419891815975070
对应$n -=1 (mod 3), n>=7$,对应边长平方为${n-1}/3+0.020168150421438454409348065456978559301$
esol(-1,-1);esol(0,-2);esol(0,-3)没有找到合法解

  1. reduce_by_z(f1,f2)={
  2.    local(f3,f4,d1,d2,r,h);f3=f1;f4=f2;
  3.    d1=poldegree(f3,z);d2=poldegree(f4,z);h=polcoeff(f4,d2,z);
  4.    while(d1>d2,
  5.         f3=f3-polcoeff(f3,d1,z)/h *z^(d1-d2)*f4;
  6.         d1=poldegree(f3,z)
  7.    );
  8.    if(d1 == d2,
  9.         f3 = f3-polcoeff(f3,d1,z)/h*f4;
  10.    );
  11.    f3
  12. }

  13. reduce_by_b(f1,f2)={
  14.    local(f3,f4,d1,d2,r,h,df);f3=f1;f4=f2;
  15.    d1=poldegree(f3,b);d2=poldegree(f4,b);h=polcoeff(f4,d2,b);
  16.    while(d1>d2,
  17.         df=polcoeff(f3,d1,b)/h;
  18.         f3=denominator(df)*f3-numerator(df) *b^(d1-d2)*f4;
  19.         d1=poldegree(f3,b)
  20.    );
  21.    if(d1 == d2,
  22.         df=polcoeff(f3,d1,b)/h;
  23.         f3 = denominator(df)*f3-numerator(df)*f4;
  24.    );
  25.    f3
  26. }

  27. esol(N1,N2)={
  28. local(f1,f2,f3,f4,f5,f6,r1,r2,r3,r4,r5, tmp,equ,rts);
  29. local(x0,y0,z0,a0,b0,c0);
  30. f1 = x+y+z-1;
  31. f2 = x*y-a*(1-a);
  32. f3 = y*z-b*(1-b);
  33. f4 = z*x-c*(1-c);
  34. f5 = x^2+a^2+(1-c)^2-y^2-(1-a)^2-b^2-N1;
  35. f6 = x^2+a^2+(1-c)^2-z^2-(1-b)^2-c^2-N2;
  36. if(poldegree(f5,a) !=1, print ("Error in r1"),
  37.     r1 = - polcoeff(f5,0,a)/polcoeff(f5,1,a);
  38.     f1 = subst(f1, a, r1);
  39.     f2 = subst(f2, a, r1);
  40.     f3 = subst(f3, a, r1);
  41.     f4 = subst(f4, a, r1);
  42.     f6 = subst(f6, a, r1);
  43.     if(poldegree(f1, x) != 1, print("Error in r2");print(f1),
  44.         r2 = - polcoeff(f1, 0, x)/polcoeff(f1,1,x);
  45.         f2 = subst(f2, x, r2);
  46.         f3 = subst(f3, x, r2);
  47.         f4 = subst(f4, x, r2);
  48.         f6 = subst(f6, x, r2);
  49.         if(poldegree(f3, y)!=1, print("Error in r3");print(f3),
  50.               r3 = - polcoeff(f3, 0, y)/polcoeff(f3, 1, y);
  51.               f2 = subst(f2, y, r3)*4*z;
  52.               f4 = subst(f4, y, r3);
  53.               f6 = subst(f6, y, r3)*4*z^2;
  54.               f2 = reduce_by_z(f2,f4);
  55.               f6 = reduce_by_z(f6,f4);
  56.               if(poldegree(f2,z)!=1, print("Error in r4"),
  57.                  r4 = - polcoeff(f2, 0, z)/polcoeff(f2, 1, z);
  58.                  f4 = numerator(subst(f4, z, r4));
  59.                  f6 = numerator(subst(f6, z, r4))-f4;
  60.                  while(poldegree(f6, b)>0,
  61.                       tmp = reduce_by_b(f4,f6);
  62.                       f4=f6;f6=tmp
  63.                  );
  64.                  f6=component(f6,1);
  65.                  if(poldegree(f4, b)!=1, print("Error in r5"),
  66.                  r5 = -polcoeff(f4,0,b)/polcoeff(f4,1,b);
  67.                  equ = factor(f6);
  68.                  for(s=1, length(equ[,1]),
  69.                     rts=polroots(equ[s,1]);
  70.                     for(t=1,length(rts),
  71.                          if(0< real(rts[t]) && real(rts[t])<1 && abs(imag(rts[t]))<0.00001,
  72.                              c0=real(rts[t]);
  73.                              b0=subst(r5,c, c0);
  74.                              z0=subst(r4,c, c0); z0=subst(z0, b, b0);
  75.                              y0=subst(r3,c,c0);y0=subst(y0,b,b0);y0=subst(y0,z,z0);
  76.                              x0=subst(r2,c,c0);x0=subst(x0,b,b0);x0=subst(x0,z,z0);x0=subst(x0,y,y0);
  77.                              a0=subst(r1,c,c0);a0=subst(a0,b,b0);a0=subst(a0,z,z0);
  78.                             a0=subst(a0,y,y0);
  79.                             a0=subst(a0,x,x0);
  80.                              if(0<b0&&b0<1&&0<z0&&z0<1&&0<y0&&y0<1&&0<x0&&x0<1&&0<a0&&a0<1,
  81.           print(equ[s,1]);                       print("a="a0";b="b0";c="c0";x="x0";y="y0";z="z0";x^2+a^2+(1-c)^2="x0^2+a0^2+(1-c0)^2)
  82.                              )
  83.                          )
  84.                     )
  85.                  )
  86.                  )
  87.               )
  88.         )
  89.     )
  90. )
  91. }
复制代码
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2018-11-5 20:21:52 | 显示全部楼层
为了确认上面的结果,我又对上面部分方程采用密集采样点加牛顿迭代法进行数值解验算,结果一致,这说明24#得出的结果对于m=3已经是最优解了
比如下面代码数值验算了25#的方程:
  1. #include <stdio.h>
  2. #include <math.h>

  3. #define MAXSOL   50
  4. #define N 6
  5. double solbuf[MAXSOL][N];
  6. int solc;
  7. #define A(i,j) a[(i)*N+(j)]
  8. void solve(double *a, double *b)
  9. {
  10.     int i,j,k;
  11.     for(i=0;i<N;i++){//in i column
  12.        double maxv = -1.0;
  13.        int index=-1;
  14.        for(j=i;j<N;j++){
  15.           if(fabs(A(i,j))>maxv){
  16.               maxv=fabs(A(i,j));
  17.               index = j;
  18.           }
  19.        }
  20.        if(index != i){//swap row index and i
  21.          double tmp;
  22.          for(j=i;j<N;j++){
  23.              tmp = A(i,j);
  24.              A(i,j)=A(index,j);
  25.              A(index,j)=tmp;
  26.          }
  27.          tmp = b[i];
  28.          b[i]= b[index];
  29.          b[index]=tmp;
  30.        }
  31.        //divid by A(i,i) in row i
  32.        for(j=i+1;j<N;j++)A(i,j)/=A(i,i);
  33.        b[i]/=A(i,i);A(i,i)=1.0;
  34.        for(j=0;j<N;j++){
  35.           if(j==i)continue;
  36.           for(k=i+1;k<N;k++){
  37.              A(j,k)-=A(i,k)*A(j,i);
  38.           }
  39.           b[j]-=b[i]*A(j,i);
  40.           A(j,i)=0.0;
  41.        }
  42.     }
  43. }

  44. void init_b(double *b, double x0, double y0, double z0, double a0, double b0, double c0, int N1, int N2)
  45. {
  46.     b[0]=x0+y0+z0-1.0;
  47.     b[1]=x0*y0-a0*(1-a0);
  48.     b[2]=y0*z0-b0*(1-b0);
  49.     b[3]=z0*x0-c0*(1-c0);
  50.     b[4]=x0*x0+a0*a0+(1-c0)*(1-c0)-y0*y0-(1-a0)*(1-a0)-b0*b0-N1;
  51.     b[5]=x0*x0+a0*a0+(1-c0)*(1-c0)-z0*z0-(1-b0)*(1-b0)-c0*c0-N2;
  52. }

  53. void init_a(double *a, double x0, double y0, double z0, double a0, double b0, double c0,int N1, int N2)
  54. {
  55.     A(0,0)=1.0;A(0,1)=1.0;A(0,2)=1.0;A(0,3)=0.0;A(0,4)=0.0;A(0,5)=0.0;
  56.     A(1,0)=y0;A(1,1)=x0;A(1,2)=0.0;A(1,3)=-1.0+2*a0;A(1,4)=0.0;A(1,5)=0.0;
  57.     A(2,0)=0.0;A(2,1)=z0;A(2,2)=y0;A(2,3)=0.0;A(2,4)=-1+2*b0;A(2,5)=0.0;
  58.     A(3,0)=z0;A(3,1)=0.0;A(3,2)=x0;A(3,3)=0.0;A(3,4)=0.0;A(3,5)=-1+2*c0;
  59.     A(4,0)=2*x0;A(4,1)=-2*y0;A(4,2)=0.0;A(4,3)=2.0;A(4,4)=-2*b0;A(4,5)=-2+2*c0;
  60.     A(5,0)=2*x0;A(5,1)=0.0;A(5,2)=-2*z0;A(5,3)=2*a0;A(5,4)=2-2*b0;A(5,5)=-2.0;
  61. }

  62. #define MAXITER 20
  63. #define MINERR  0.00001
  64. #define DELT    (4*MINERR)
  65. void pushs(double x0, double y0,double z0, double a0,double b0,double c0)
  66. {
  67.     int i;
  68.     for(i=0;i<solc;i++){
  69.        if(fabs(x0-solbuf[i][0])<DELT&&fabs(y0-solbuf[i][1])<DELT&&
  70.           fabs(z0-solbuf[i][2])<DELT&&fabs(a0-solbuf[i][3])<DELT&&
  71.           fabs(b0-solbuf[i][4])<DELT&&fabs(c0-solbuf[i][5])<DELT){
  72.             break;//find a matched solution
  73.        }
  74.     }
  75.     if(i>=solc){
  76.        if(solc>=MAXSOL){
  77.           fprintf(stderr, "buffer overflow\n");
  78.           return;
  79.        }
  80.        solbuf[solc][0]=x0;solbuf[solc][1]=y0;solbuf[solc][2]=z0;
  81.        solbuf[solc][3]=a0;solbuf[solc][4]=b0;solbuf[solc][5]=c0;
  82.        printf("x=%f,y=%f,z=%f,a=%f,b=%f,c=%f\n",
  83.             x0,y0,z0,a0,b0,c0);
  84.        solc++;
  85.     }
  86. }

  87. int main(int argc, char *argv[])
  88. {
  89.     int N1=atoi(argv[1]);
  90.     int N2=atoi(argv[2]);
  91.     int count;
  92.     double a1,b1,c1;
  93.     double a[N*N],b[N];
  94.     for(a1=0.1;a1<=0.9;a1+=0.1)for(b1=0.1;b1<=0.9;b1+=0.1)for(c1=0.1;c1<=0.9;c1+=0.1){
  95.         double a0=a1,b0=b1,c0=c1,x0,y0,z0;
  96.         double r=sqrt(a0*(1-a0)*b0*(1-b0)*c0*(1-c0));
  97.         z0=r/(a0*(1-a0));x0=r/(b0*(1-b0));y0=r/(c0*(1-c0));
  98.         if(0<x0&&x0<1&&0<y0&&y0<1&&0<z0&&z0<1){
  99.           for(count=0;count<MAXITER;count++){
  100.             init_a(a,x0,y0,z0,a0,b0,c0,N1,N2);
  101.             init_b(b,x0,y0,z0,a0,b0,c0,N1,N2);
  102.             solve(a,b);
  103.             double maxerr=-1.0;
  104.             int u;
  105.             for(u=0;u<N;u++){
  106.                if(fabs(b[u])>maxerr)maxerr=fabs(b[u]);
  107.             }
  108.             x0-=b[0];y0-=b[1];z0-=b[2];a0-=b[3];b0-=b[4];c0-=b[5];
  109.             if(0<x0&&x0<1&&0<y0&&y0<1&&0<z0&&z0<1&&
  110.                 0<a0&&a0<1&&0<b0&&b0<1&&0<c0&&c0<1){
  111.                  if(maxerr<MINERR){
  112.                      
  113.                      pushs(x0,y0,z0,a0,b0,c0);
  114.                      break;
  115.                  }
  116.             }else{
  117.                 break;
  118.             }
  119.           }
  120.         }
  121.     }
  122. }
复制代码
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

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

GMT+8, 2019-1-23 05:57 , Processed in 0.052687 second(s), 16 queries .

Powered by Discuz! X3.4

© 2001-2017 Comsenz Inc.

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