Top

個別メモ

Latest update on 2012年3月5日 (月) at 10:54:46.

【第1473回】 一木会とか（2009年12月3日）

• 6:30起床。往路あさま510号にはダッシュで何とか間に合った。昨日でコーヒー豆が尽きたので，新前橋から大学に行く途中で県庁前の珈琲問屋に寄って豆を買っていくつもり。
• 今日の予定は，昼過ぎに留学生の研究相談に乗るのと，夜の一木会以外は採点と来週の口演準備か。
• とりあえずEulerの最新版をダウンロードしてインストールしてみた。起動してみると，プロンプトはRと同じケット（半角の大なり記号，つまり>）である。Eulerの計算式はRと同じように使える。
• とりあえずζ(2)が円周率の2乗を6で割った値に等しいことをみるには，sum(1/(1:10000000)^2)とpi^2/6が，ともに1.644933966847と表示されることでわかる。一千万項の計算だが，Let's Note Y5でもほぼ一瞬で終わる。Rではデフォルト表示桁数が小数以下7桁なのだが，options(digits=13)とするとEulerと同じ表示になる。この計算のための入力は，RとEulerでまったく同じである。
• このとき，計算する項数xを増やすにつれてsum(1/(1:x)^2)はpi^2/6に漸近していくのだが，その様子を1000項の場合までグラフ表示するには，RでもEulerでも，plot(cumsum(1/(1:1000)^2))とすればいい。ただしRではこれだと線でつながれずに点々がプロットされるので，線でつなぐためには，plot(cumsum(1/(1:1000)^2),type="l")とすればいい。Eulerでやる方法はまだ知らないのだが，Rの場合は，縦軸や横軸を対数軸にするのも簡単なので，項数を増やしたときにpi^2/6に近づく様子を見やすくするために横軸を対数軸にしてみるには，plot(cumsum(1/(1:1000)^2),type="l",log="x")とする。その状態でlines(c(1,1000),c(pi^2/6,pi^2/6),col="red")を実行すれば，右の図ができる。何日か前に書いたgslライブラリをlibrary(gsl)として呼び出してから，zeta(2)とすればpi^2/6と打たなくてもいい。
• 次にゼータ関数の引数を2以上の整数で徐々に大きくしていくとどうなるかを試してみる。もちろん，賢いgslライブラリを使えばplot(2:1000,zeta(2:1000),type="l",log="x")などということも簡単にできるのだが（どうやっているのかわからないが非常に速い），定義式通りにx乗分の1を計算させてzeta()と同じ精度を出すために一千万項まで計算させてしまうと，ζ(100)まででもハングアップしたかと思うくらいの途方もない時間がかかるのでやめておく方が無難である。計算させる項数を1万項くらいにとどめてζ(2)からζ(100)までの近似をプロットさせるには次のようにする。まずべき乗関数powを定義する。pow<-function(x,y) x^yとすればいい。次にplot(2:100,rowSums(sapply(1:10000,pow,-(2:100))),type="l")とすると，横軸に2から100まで，縦軸にそれに対応するゼータ関数の近似値をとったグラフが描けるわけである。この状態でlines(2:100,zeta(2:100),col="red",lty=2)とすると，gslライブラリのゼータ関数で計算した値が重ね描きされる。パソコン画面くらいの解像度では，2つのグラフのズレはほとんどわからないくらいである。以上，とりあえずRやEulerでも実数域ならそこそこζ関数をいじくれるというメモ。
• 思い出したのでついでに書いておくと，『算数宇宙の冒険　アリスメトリック！』の算数宇宙杯の1回戦とかで，ユーキが不満をもらす普通の旅人算などが使われているところで，ぼくが以前思いついたダンゴムシクモ算とかフラミンゴウサギ算を取り上げてもらえたら嬉しかったなあ。あの小説の読者ならば，きっとこういうくだらない遊びは好きなんじゃないかと思うのでリンクしておく。
• 午前中に事務補佐のＹさんがマークシート解答の読み取り作業をしてくれて，助教のＫ先生が記述の採点をしてくれたので，これから自分が記述の採点をすれば，採点作業のゴールはかなり近づくことになる。頑張ろう……と思ったところで昼になった。メールのやりとりが時間をくうなあ。
• 漸く採点が終わったのが18:15だったので，結局口演準備には入れなかった。まあ1つでも仕事が片付いたのでいいことにしよう。一木会はまあ普通に終わった。簡単に漢方が効くタイプと効かないタイプのRAを見分けられるマーカがもし見つかるなら，数百万くらい安いものだと思うが，問題はその戦略で違いが見つからなかった場合とか，あるいは見つかったマーカでは感度や特異度が低いという場合だろう。
• 雨が降っていたので自転車が使えないため，バス停に直行して前橋へ。復路あさま549号。7月に発行された『たぁくらたぁ』の給食特集を読んだが，実に読み応えがあった。真田町や塩尻市や小諸市や宮田村は素晴らしい取り組みをしていると思う。
• うーん，青木先生からいただいたコメントを見て再度確認してみたら，上で書いたことはたぶん間違っていて，一千万項くらいではpi^2/6には収束せず微妙に差が出てしまうようだ。しかし不思議な差の出方だなあ。計算機特有の誤差なのか？
• 昨夜見つけたJavaのページに書かれていた式をそのままRに置き換えてZeta(s)を関数定義したら（但しもちろん無限和はできないので柏野さんに倣って500項の和にして），複素数を受け付けるゼータ関数ができた。もっとも，Zeta(-1)が-0.0833313となってしまうしZeta(2)は1.644930074848となってしまうので近似精度は低いし，実はRのガンマ関数gamma()が複素数を引数にとれないので実部が負の場合はまだ計算できないのだが，実部が正の場合，例えばZeta(1/2+1i)は計算できて， 0.144659-0.6974934iとなる。もし，これが正しければ，柏野さんのサポートサイトに載っていたグラフが作れるはずだが，どうなんだろう。とりあえず複素数対応版のガンマ関数を作らねば。
• 完成度が低いものをダラダラと載せて恐縮だが，ここはメモ書きなのでご容赦いただきたい。

▼前【1472】（卒業試験の日（2009年12月2日） ）　▲次【1474】（ひどい雨の朝（2009年12月4日） ）　●Top

△Read/Write COMMENTS

Notice to cite or link here | [TOP PAGE]