古今算法記の問題4をkeisan.casio.jpで解いてみる。

数学セミナー2008年12月号の「現代数学から見た関孝和」に出てた問題。

をx,y,zについて解くというもの。未知数を消去するとz^1/4の108次式になるそうです。でこれExcelのソルバーで解けないかな？とおもってやったら全然値が違う。あー、桁数が全然足りてないんだ。ということでここはkeisan.casio.jpでやってみよう（これはUPしてもしょうがないかな）。

解くと、30桁で

x=36.8641465377853041287298140705
y=44.3516761876674573034940403248
z=32.5563821163811031295834914427

x=51.4570363183687277311149213375
y=10.293599457590667227853861756
z=49.4144096308403396077899331929

の2つが得られた。解き方はニュートン・ラフソン法。でも微分して逆行列を求めるのは死ぬほどめんどくさい。なのでやっぱりMaximaでやってみた。keisan.casio.jpに対応したプログラムはこちら。ね？めんどくさいでしょ。

for(k=0;k<=1;k=k+1){
  if (k==0) {
    x=36;
    y=44;
    z=32;
  } else {
    x=51;
    y=10;
    z=49;
  }
  for(j=0;j<=100;j=j+1) {
    x2=x-(-(3*sqrt(x)*y^(8/3)*(z^3+y^3-121750))/((12*y^(8/3)-8*x^(5/2))*z^(11/4)+6*x^(5/2)*y^(8/3))-((y^3+x^3-137340)*(4*sqrt(x)*z^(11/4)-3*sqrt(x)*y^(8/3)))/((12*y^(8/3)-8*x^(5/2))*z^(11/4)+6*x^(5/2)*y^(8/3))+(12*sqrt(x)*y^(8/3)*(z^(1/4)+y^(1/3)+sqrt(x)-12)*z^(11/4))/((6*y^(8/3)-4*x^(5/2))*z^(11/4)+3*x^(5/2)*y^(8/3)));
    y2=y-((3*x^(5/2)*y^(2/3)*(z^3+y^3-121750))/((12*y^(8/3)-8*x^(5/2))*z^(11/4)+6*x^(5/2)*y^(8/3))-(12*x^(5/2)*y^(2/3)*(z^(1/4)+y^(1/3)+sqrt(x)-12)*z^(11/4))/((6*y^(8/3)-4*x^(5/2))*z^(11/4)+3*x^(5/2)*y^(8/3))+(3*y^(2/3)*(y^3+x^3-137340)*z^(11/4))/((6*y^(8/3)-4*x^(5/2))*z^(11/4)+3*x^(5/2)*y^(8/3)));
    z2=z-(((3*y^(8/3)-2*x^(5/2))*z^(3/4)*(z^3+y^3-121750))/((6*y^(8/3)-4*x^(5/2))*z^(11/4)+3*x^(5/2)*y^(8/3))+(12*x^(5/2)*y^(8/3)*(z^(1/4)+y^(1/3)+sqrt(x)-12)*z^(3/4))/((6*y^(8/3)-4*x^(5/2))*z^(11/4)+3*x^(5/2)*y^(8/3))-(3*y^(8/3)*(y^3+x^3-137340)*z^(3/4))/((6*y^(8/3)-4*x^(5/2))*z^(11/4)+3*x^(5/2)*y^(8/3)));
    x=x2;
    y=y2;
    z=z2;
  }
println(x2,y2,z2);
}