古今算法記の問題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);
}
« 妖怪ラーメンは手作り感満載 | トップページ | 伏見ででかい釜を見た。 »
「学問・資格」カテゴリの記事
- 高周波(RF・マイクロ波・ミリ波・5G)関連ニュース2021年2月16日 IEEE Microwave Magazineの特集はオールデジタルのRFID、Microwave JournalはEバンド ミリ波通信に衛星や気球を使う話、アメリカの半導体企業がバイデンに投資を迫る、(2021.02.17)
- カオスを生じる電気回路、Chua’s circuitをLTspiceで回路シミュレーションしてみる。(2021.02.19)
- Labyrinth Chaos(迷宮カオス)を生むThomas-Rössler方程式のパラメータbを色々変えて、Python+Scipyでルンゲクッタ8次のDOP853(Dormand&Prince)を使って計算してGIFアニメ(2021.02.16)
- フィッツヒュー・南雲 (FitzHugh-Nagumo) 方程式をPython+Scipyでルンゲクッタ8次のDOP853(Dormand Prince)で計算。(2021.02.23)
- 「水晶振動子の等価回路計算」をカシオの高精度計算サイトkeisan.casio.jpの自作式としてUP! インピーダンスの大きさと位相がグラフ化できる。(2021.02.12)
コメント