« 少し量子化学に近い話:さらにHaskell | トップページ | さらにまたoctave »

2009年7月23日 (木)

再びmaxima

今日は、昨日書いたorbの出力をもとに、plot3dのプロット範囲をどう決めるか、書いてみます。

6Sに対応するorbの出力は、こうでした。

$ ./orb 6 0
[1,-90,2700,-32400,145800,-174960]

これを基に、こういう式を作って、maximaに代入してみます。

diff((exp(-r/6)*(r^5-90*r^4+2700*r^3-32400*r^2+145800*r-174960))^2*r^2,r);

こういう結果になります。

2 5 4 3 2 2 - r/3
r (r - 90 r + 2700 r - 32400 r + 145800 r - 174960) %e
(%o1) - -----------------------------------------------------------------
3
5 4 3 2 2 - r/3
+ 2 r (r - 90 r + 2700 r - 32400 r + 145800 r - 174960) %e
2 4 3 2
+ 2 r (5 r - 360 r + 8100 r - 64800 r + 145800)
5 4 3 2 - r/3
(r - 90 r + 2700 r - 32400 r + 145800 r - 174960) %e

ちょっと見にくいので、ratsimp(%);と入力して整理してみます。

12 11 10 9 8 7
(%o2) - (r - 216 r + 19440 r - 955800 r + 28285200 r - 523480320 r
6 5 4 3
+ 6101205120 r - 44026934400 r + 188484408000 r - 442158912000 r
2 - r/3
+ 489776025600 r - 183666009600 r) %e /3

ペーストしたらべき乗の数値がずれているようなので、ずらして読みなおしてください。

なにをやっているのかというと、関数が極大、極小になるを求めるために、微分してみたわけです。

それで、このrの多項式の部分が0になるところを求めたいのですが、exp(-r/3)が余計なので、exp(r/3)をかけることにして、

realroots(%o2*exp(r/3));

と入力してみます。これは多項式のの実数解を求める関数です。こうなります。

24395313 62112359 130477793 212698117
(%o4) [r = --------, r = --------, r = ---------, r = ---------,
33554432 33554432 33554432 33554432
330139973 464141663 643351957 845477765 1118977307
r = ---------, r = ---------, r = ---------, r = ---------, r = ----------,
33554432 33554432 33554432 33554432 33554432
1435468975 1980516091
r = ----------, r = ----------, r = 0]
33554432 33554432

最後の、r:=0の前のものが求めたい極大の場所(6Sだと一番外側のリングのところ)になっています。やっぱりペーストすると、すっかりずれますね。

maximaは数値計算は苦手なようなので、r=1980516091/33554432を、電卓なり表計算なりで求めてみると、

59.023979

くらいになります。ただ、これをプロット範囲にすると、リングのてっぺんで描画が切れることになるので、ルート2くらいをかけてやって、

83.472512

が適当でしょう。

plot3d(
[p, q,
(exp(-sqrt(p*p+q*q)/6)*(((sqrt(p*p+q*q) -90) * (p*p+q*q) + 2700 * sqrt(p*p+q*q) - 32400) * (p*p+q*q) + 145800 * sqrt(p*p+q*q) - 174960))^2*(p*p+q*q)*0.000000000744904*0.5],
[p,-83.4725,83.4725],[q,-83.4725,83.4725],['grid,160,160]);

のp,qのプロット範囲は、こうやって求めていたわけです。

|

« 少し量子化学に近い話:さらにHaskell | トップページ | さらにまたoctave »

コメント

この記事へのコメントは終了しました。

トラックバック


この記事へのトラックバック一覧です: 再びmaxima:

« 少し量子化学に近い話:さらにHaskell | トップページ | さらにまたoctave »