数学セミナー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);
}
最近のコメント