「告白」（湊かなえ）を読んだ。

映画化されるということで最近よくCMが放映されるのを見てて、かつ文庫化されていたので読んでみた。

これは。。。確かに重いなあ。救いが全くないけれど、何かまた読み返したくなるようなお話。読むこと自体は1時間半くらいですっと読めるんですが。

しかし映画のCMの携帯電話でどっかーんというのは完全にネタばれではないのか、あるいは結末が違うのか。。。

「月光ゲーム Yの悲劇'88」を読んだ。

だいぶ前に買っていたのだが、オリジナルのYの悲劇を読み返してからのほうがいいかな、ということで置いていた。最近オリジナルを読んだのでようやく読んだ、が、オリジナルとは特に関係なかったorz。

大学生が今と違ってみんな純朴なのがほほえましいというかちょっと恥ずかしいというか。作者の長編第一作ということですが、なかなかおもしろかった。でもダイイング・メッセージはなくてもいいんじゃないか、というか、小学生のなぞなぞというか。

アレンストーフ軌道を14次ルンゲクッタ法で計算（PARI/GP)

PARI/GPで35段14次ルンゲクッタを計算するシリーズ。アレンストーフ軌道も高精度に計算できるかな？と思ってやってみた。

y1" = y1 + 2* y1' - μ'* (y1 + μ)/D1 - μ*(y1 - μ') / D2

y2" = y2 - 2* y2' - μ' *y2/D1 - μ*y2 / D2

D1 = ((y1+μ)^2+y2^2)^(3/2)

D2 = ((y1-μ')^2+y2^2)^(3/2)

μ=0.012277471

μ' = 1-μ

y1(0) = 0.994, y1'(0)=0,

y2(0)=0,y2'(0)=-2.00158510637908252240537862224

計算結果はこちら。

ソースはこちら。↓

「aren.gp」をダウンロード 

「トイレの神様」は反則です！

さっき探偵！ナイトスクープを見ていたら植村花菜さんのトイレの神様で号泣する4歳児が出ていた。

こういう曲があるのは聞いたことがあったけど、聴くのは初めてだった。で、私も号泣。。。私も小さいころはおばあちゃん子だった。両親が共働きでほとんど構ってもらえず、ご飯もおばあちゃんと二人で食べることが多かった。いつも話を聞いてもらって、にこにこしてもらって、本当に大好きだった。おばあちゃんが若いころ、洋裁の勉強に京都にいたことがあって、そのとき当時の京大生に言い寄られた（けれど振った）ことも聞いたことがあって、私が京大に入ったらだれよりも喜んでくれた。おばあちゃんは非常に高貴な家の出だったけれど、だまされてうちのじいさん（会ったことなし）のところに嫁いで、じいさんが稼がなくて超貧乏で苦労した。あのとき京大生とくっついていたら苦労せずにすんだ。。。（でもそれなら私生まれないんですが）。

おばあちゃんは私が大学生のころ90歳くらいで死んでしまったんですが、やっぱりおばあちゃん孝行まったくできなかった。超後悔。京大に入ったことだけ。。。

こんなの反則だよ。泣くにきまってるじゃないか。

誕生寺→音声館あたりをぶらぶら

ならまちあたりです。

音声館。

ここらへんはものすごくお寺が多いです。写真撮り忘れたけれど。。。

ルンゲクッタ14次でJapanese Attractorを描いてみる（PARI/GP)。

Duffing方程式のストロボプロットでジャパニーズ・アトラクタを描いてみよう。35段14次のルンゲクッタ法で、PARI/GPを使う（グラフはExcelで）。60桁精度で。

計算している式は

dx/dt = y

dy/dt = -k*y - x^3 + B*z

z=cos(t)

k=0.05, B=7

初期値はx(0)=2.1, y(0)=0で。

こういうカオスを生じるものはどんなに精度上げて計算してもいつかは破綻するわけですが、4段4次のルンゲクッタで倍精度よりは（このくらいの時間なら）圧倒的に正確に描けているはず。とはいえ当時、上田さんがこれを出したときは精度がどうのこうの、とか言って日本の当時の権威のような老人方からは黙殺されたそうですが（で海外で認められて逆輸入ということで、情けない話）。

はねるのトびらのしりたしゲームコーナーの中近東風のBGMはアバのヴーレ・ヴー

たまたまテレビ見てたら、ABBAのVoulez　Vousがかかっていた。

いつ聴いてもいいアレンジ。でも何故これを選んだんだろうか。

「鳴風荘事件 殺人方程式II」を読んだ。

これはですね、同じ症状の人が私の同期にいまして、仕事をお話と同じように変わったんです。それですぐわかったりして。これだけではネタばれではないだろう。

館シリーズと比べてさわやかなんですが、やっぱりおどろおどろしいものも読みたいですね。新しい館シリーズいつでるんだろうか。

ルンゲクッタ法4次、8次、14次の比較

先日、PARI/GPを導入して35段14次のルンゲクッタ法のプログラムを書いてみた。でも全然精度のチェックとかしていなかったのでここでやってみる。

比較するのはこのブログでは最多登場の13段8次のルンゲクッタであるDormand-Princeと、ルンゲクッタと言えば普通はこれ！という4段4次のもの。

計算したのは単振動の方程式。

dx/dt = v

dv/dt=-x

時間刻みをかなり大きくΔｔ＝0.1ととった。

で厳密解のcos(t)と比べたのがこちら。

想像通り、14次が一番精度高い（ただしExcelに移したときに桁が減っていて、本当はもっといい）。でも時間は想像以上にかかる。

35段14次のルンゲクッタ法をPARI/GPで計算(Lorenz方程式)

先日導入したPARI/GPに35段14次のルンゲクッタ法を実装してみた。

係数はこちら↓60桁分あります。

http://sce.uhcl.edu/rungekutta/

まずはいつものようにローレンツ方程式をやってみよう。式はこちら。

dx/dt = -σ*(x-y)

dy/dt=-y-x*z+r*x

dz/dt=x*y-b*z

でプログラムのソースコードはこちら。60桁で計算してますよ。

「Lorenz.gp」をダウンロード

使い方は、PARI/GPを立ち上げて、

\r Lorenz

でロード。

rk14("ファイル名")

で実行されて、ファイル名のところに保存される。散布図がどうも描きにくそうなので、ファイルに落としてExcelで描いてみた。それがこれ。

すばらしく安定感のある計算ができるなあ。これからはこれでしばらくいろいろ計算していってみよう（続く）。

任意精度実数計算用にPARI/GPを導入してみる。

最近、35段14次の高次ルンゲクッタ法をExcel VBAで計算するシリーズをやってみたけれど、VBAが倍精度実数しか使えないので精度ではいまいちだった。そこで多倍長実数（任意精度実数）か、そこまで言わなくても4倍精度実数が計算できるプログラミング言語がないかな、と少し探す。CのライブラリならすぐGMPとか

http://gmplib.org/

出るんですが、言語のほうがおもしろそう（うちのブログでやっているのは単に私の暇つぶしなので、覚えることが多いほうがいいということで。。。）

昔から知っているのはUBASIC

http://www.rkmath.rikkyo.ac.jp/~kida/ubasic.htm

なんですが、Windows版がないしちょっと古くなった感はどうしても否めない。同じくBASIC系で、（仮称）十進BASIC

http://hp.vector.co.jp/authors/VA008683/

も249桁できるということ。でも日本語の文献があり過ぎる！（いや、暇つぶしだから、文献はないほうが。。。）

で、見つけたのがPARI/GP。

http://pari.math.u-bordeaux.fr/

ダウンロードしてインストールはすぐにできた。

\p 100

とかすると、100桁の精度になります。

（表示が長すぎるんで、表示のみ短くするのは、default(format, "g0.28");とか。）

でPiと打つと、

3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068

と表示。簡単。

繰り返しがCと似てる分混乱するんだが、

for (i=1, 5, print(i))

見たいな感じ。最初for(i=1, i<=5,...)とか書いて全然動かなかった。ネストも（）で行うんで、

for (i=1,5,

   for (j=1,5,

      ...

  )

)

みたいな。関数の定義も特にreturnで返すんじゃなくて

f(x) = {

x;

}

なんて感じ。この関数内だけで使う変数はlocal(a,b);とか。

ベクトルと行列は

v=vector(5,i,0.);         m=matrix(5,5,i.j,0.);

で定義して、v[1]やm[1,2]とか。これも微妙にCと違うので混乱。また添え字は1から始まるのも微妙に違う。

でもこれだけあればもうルンゲクッタくらい書けそうだ。やってみよう（続く）。

「Yの悲劇」を読んだ。

この前帰省した時に大昔（小学生のころ）、買っていたこの本を読んでみた。とはいえ全く内容覚えてない。というかこれを小学生が読むのは無謀だったのでは。

しかし今読むとだれも狂ってるようには見えないなあ。それだけ世間の方がおかしくなってしまったのか。。。犯人も当時は衝撃的だったと思うけれど、今見るとあってもおかしくないというか、途中でわかった（もしかしたら昔の記憶が戻っただけかもしれないが）。

生田神社でお参り。

あの生田神社ですよ。あの。

芦屋の交通安全地蔵尊

突然現れます。お参りしてきました。

JAFの看板＠芦屋

おなじみのマークがでかくなるとちょっと違和感が。。。

バーガース方程式をExcel VBA&擬スペクトル法&高次ルンゲクッタ法で計算

KdV方程式、蔵本・シバシンスキー方程式に続いて、Burgers方程式だ（順番が前後したような気がするが。。。）

空間部分は擬スペクトル法、時間部分は高次ルンゲクッタ法（5次のCash-Karp法）で計算してみる。

結果はこちら。だんだん突っ立ってくる。

とにかく今まで計算してきたものでは一番精度が高い方法だと思うが、何せ時間がかかる。。。

蔵本・シバシンスキー方程式をExcel VBA&擬スペクトル法&高次ルンゲクッタ法で計算

昨日はKdV方程式を計算してみたけれど、今日は Kuramoto-Sivashinsky方程式をやってみよう。空間部分は擬スペクトル法、時間部分は高次ルンゲクッタ法（5次のCash-Karp法）で計算してみる。

これはさすがに時間がかかった。参考にしたのは「常微分方程式の数値解法II」に出ていた方法（コレ-エックマン-エプスタイン-シュトゥッベ）です。文献ではN=1024になっているんだが、そんなことしたらすさまじく時間がかかるのでN=200くらいで。

結果はこちら。なかなかきれいな絵がかけた。

KdV方程式をExcel VBA＆擬スペクトル法＆高次ルンゲクッタ法で計算

この前、新PCを購入したので今まで時間がかかるかと思って敬遠していたような計算をしてみよう。まずは偏微分方程式を空間部分は擬スペクトル法、時間部分はルンゲクッタ法（ここでは6段5次のCash-Karp法）という手法を使ってExcel　VBAでプログラムを書く。題材は一番ポピュラーなKdV方程式にしてみる。

で早速計算した結果がこちら。

オリジナルのザブスキーとクルスカルの数値計算は蛙飛び（の工夫したもの）で、それより今回の結果はやっぱりパルスがシャープに計算できてる。でも逆にあの時代にああいった計算をしたのはすごいな。前にうちのブログでも3次精度の風上差分で計算したときはパルス鈍りまくりだったし。

時間発展を追っていったものはこちら。これもきれいに計算できてる。

「どんどん橋、落ちた」を読んだ。

5つの短編集ですが、何を言ってもネタばれになりそうな話。新本格派は「人が書けてない」とか批判されることのパロディにもなってます。中でもやっぱりすごいのは「伊園家の崩壊」で、これは「サザエさん」の殺人事件。サザエさんにこういうパロディがあるのは聞いたことがあったけど、まさか綾辻さんが書いていたとは！後味はひじょーーに悪いですが、面白いのは面白い。しかし、これで推理小説を読むときは殺”人”事件とか、犯”人”とかが書いてあるかどうか細かく読むようになるかも。

中山寺に久々に行ってきた。

仁王像が改修されてきれいになっていました。

でいちばんてっぺんのここにもお参り。

焼肉・小倉優子を今度は山科・醍醐付近で見つけた。

京都に異常に多いなあ。5件くらい通りかかった。一件も入ったことはないけれど。。。

亀が多すぎ

小倉付近の川ですが、ここから100メートル以上に渡って亀が異常発生。何があったんだ？

DellのデスクトップPC(Studio XPS 8100)を買った。

前に使っていた古いEPSONのデスクトップPCが壊れて、ずっとノートPC（Let's note B5サイズ）を使っていたのだが、さすがにメインマシンとしては使いにくいので思い立って購入。

CPUは Core i7 2.8GHz, OSはWindows 7 Ultimate, メモリは8GB, ハードディスクは1.5TB。その他NVIDIAのGeForce GTS 240,ブルーレイディスクドライブ、とかとりあえずそれなりのものは全部つけてみた。で203,470円なり（ディスプレイは別）。

基本的にはネットワークに繋いで立ち上げたらもう設定とかなにもいらなかった。楽だけど、昔の大変さを知っている身にはちょっと物足りない。

ではこれでExcelで計算シリーズとかの規模が倍増できそうだなあ。すこし考えてみよう。

「時間封鎖」を読んだ。

まったく予備知識なしで、Amazonでお勧めされていたこの本読んでみた。最初は「宇宙消失」+「竜の卵」+「幼年期の終わり」みたいな話かな、と思ってたんだけど（確かにそういうところはあるけれど）これはまた別の面白さですよ。

特に時間の進み方が遅れている理由が、人間側の問題でなくて仮定体側の問題（ってネタバレ？）ってのが面白かった。いいやつじゃないか！仮定体は。それに比べて人間のあさましいことというか。。。

あと、薬を飲んで。。。というところは双子の弟が飲んだら姉になるんだとずっと思ってた。一回生まれ変わって。これじゃメルモちゃんですね。まったく違ってました。

もっとも意外なのが母親に手紙を書いた主で、こういうところをちゃんと書くから深くなるんですね。まだ第二部、三部と続くそうなんで楽しみ。

カレーに生卵（続き）

昨日、ドラマ「怪物くん」でカレーに生卵を入れる話が出ていてうちの母親も同じことをしていたのを書いたけれど、生卵と言えばいろいろあることを思い出した。

まずは、有名どころで「オロナミンC+生卵（の黄身）＝オロナミン・セーキ」です。大村崑さん（実はこの方、今はうちの地元にすんでいるそうです）のCMでおなじみ。公式サイトでも紹介されている。

http://www.otsuka.co.jp/orc/contents/recipe.html

次は「黒ビール+生卵（の黄身？）」です。これはうちの親父が大昔にこれがうまい、と言っていたもの。しかし飲んだことがないし、うちの親父も飲んでいるのをみたことがない（いわく、黒ビールは高い）。調べると

http://portal.nifty.com/2007/12/26/a/

とかあった。でも合わないって？イエローサンセットって名前だそうですが。

しかしTKG（たまごかけごはん）も贅沢に黄身しか使わない人がいる。大反対。白身も大事にしよう。小栗旬くんが味の素のCMで食べているちゃんと白身も使ったたまごかけごはんは最高にうまそうです。

http://www.ajinomoto.co.jp/cm/tvcm/ajinomoto.html

あと卵酒があるんですが。。。

私これ今まで何十年か生きてきて一回も飲んだことがない。作り方は

http://portal.nifty.com/koneta06/01/30/01/

とかありましたが、うちは風邪ひくともう生トマトとかしかなかった。。。

あと納豆に卵いれます？私は昔あまり納豆が得意じゃなかったころは入れていたけれど、好きになってからは納豆の風味が減るからあんまり入れなくなった。ネギも。しょうゆとからしだけでOK。あとは釜玉うどん最強とか。

三栖神社でお参り。

中書島付近にて。結構さびしいところにあります。とりあえずお参り。

カレーに生卵

さっき「怪物くん」を見ていたら、カレーに生卵を混ざるエピソードがあった。これ、実はうちの母親が昔、全く同じことをしていた。子供のころ、「なんで？」と聞いたところ、味がマイルドになる、私の親戚のおばさんもやっていたのを真似た、ということらしい。関西では結構やっている人が多いらしい。でも私はパス。

私の知り合いで納豆をかける人もいたが、これもパス。。。

でも有名なカレー店、バーグでは生卵をかけるらしい。一度行きたいがちょっと遠い。。。

http://r.tabelog.com/kanagawa/A1403/A140308/14000923/

相対論的水素原子のエネルギー準位（Dirac方程式の解）をkeisan.casio.jpにUP!

実は細かい数字まで追ったことがなかったりして。とりあえず自作式をUPしてみた↓

相対論的水素原子のエネルギー準位（Dirac方程式の解）

内容は以下に：

相対論的量子力学の基礎方程式であるDirac方程式による水素原子(Z=1)のエネルギー準位を計算します。非相対論的なシュレーディンガー方程式と比較しています。

Javascriptで複素数を扱うライブラリ

※この記事は古すぎてもう動かないのでこちらを参照。

TypeScriptで複素数クラスを作って複素数の四則演算、exp、べき、logの計算ができるようにした。大昔にJavaScriptで作ったり、数値計算ライブラリmathjsを使ったりしたがもうちょっとお手軽なものが欲しかったので。 

以下は記録としてだけ。

ーーーーーーーーーーーーーーーーーーーー

練習を兼ねて作ってみた。ソースリストは最後のほうにつけてます。使い方は

(作成） a = new Complex(1,1);

(足し算)　c = Complex.cadd(a, b);

（引き算）　c = Complex.csub(a, b);

（掛け算） c = Complex.cmul(a, b);

（割り算） c = Complex.cdiv(a, b);

（べき乗） c = Complex.cpow(a, b);

（exp）    c = Complex.cexp(a);

(log)     c = Complex.clog(a);

とかでOK。

---------------------

function Complex(x, y) {
    this.x = x;
    this.y = y;

    return this;
}

Complex.cadd = function (a, b) {
    return Complex(a.x+b.x, a.y+b.y);
}

Complex.csub = function (a, b) {
    return Complex(a.x-b.x, a.y-b.y);
}

Complex.cmul = function (a, b) {
    return Complex(a.x*b.x-a.y*b.y, a.x*b.y+a.y*b.x);
}

Complex.cdiv = function (a, b) {
    r2= b.x*b.x + b.y*b.y;
    return Complex((a.x*b.x+a.y*b.y)/r2, (-a.x*b.y+a.y*b.x)/r2);
}

Complex.cconj = function (a) {
    return Complex(a.x, -a.y);
}

Complex.cexp = function (a) {
    return Complex(Math.exp(a.x)*Math.cos(a.y),Math.exp(a.x)*Math.sin(a.y));

}

Complex.clog = function (a) {
    r2= a.x*a.x + a.y*a.y;
    return Complex(0.5*Math.log(r2), Math.atan2(a.y,a.x));
}

Complex.cpow = function (a, b) {
    return Complex.cexp(Complex.cmul(Complex.clog(a),b));
}

function cabs(a) {
    return Math.sqrt(a.x*a.x+a.y*a.y);
}

function carg(a) {
    return Math.atan2(a.y,a.x);
}

「平家伝説殺人事件」を読んだ。

うちの実家も私が生まれるはるか前、伊勢湾台風で流されたということで人ごとでないお話。。。

しかし浅見さんのシリーズは、最後は寅さんのように女性とうまくいかないはずなのに、このお話では120％うまくいってるじゃないですか（笑）。これはなかったことになってるのだろうか？

こいのぼり列

こどもの日ということで、こいのぼりの写真。うちの地元のそうとう田舎のほうなんですが、毎年すごい数のこいのぼりが一列に並んでいます。

「サマーウォーズ」をBluray版で見た。

PS3を買ったので、ようやく我が家もブルーレイディスクが使えるようになった。自慢の42インチのプラズマテレビ（ビエラ）とHDMIケーブル（高かった。。。）で接続。で最初に観るものとして、前から気になっていたサマーウォーズを選んだ。

これは面白いですよ。おばあちゃんがとにかく最高。泣かせどころも何箇所か用意されていて、劇場で観ていたらヤバかったかも。。。あとカズマは話の後半までずっと女の子だと思っていた。

おまけディスクではCG担当ディレクターのインタビューがありましたが、アバターは100人のデザイナーが250体くらい本当にデザインしたとのこと。またラブマシーン巨大版も本当に粒子シミュレーションでたくさん描いてるそうです。大変だ。。。

では他ではあまり語られていないようなお話を。

主人公が電車の中で読んでいたのが暗号の伏線でして、Shorの因数分解アルゴリズム。これは量子コンピュータ向けのもんなのですがそれを暗算でやるとは。。。スーパーコンピュータが出てきたのも暗号を解くためでなく、格闘に勝つためというのもすごいな。

だから暗号が解けるわけないという突っ込みはもう最初の段階で成り立ってないですよ。なんせ頭が量子コンピュータ（というかグレッグ・イーガンの小説で出てきた波束の収束マシーンか？）なんだから。

http://ja.wikipedia.org/wiki/%E9%87%8F%E5%AD%90%E3%82%B3%E3%83%B3%E3%83%94%E3%83%A5%E3%83%BC%E3%82%BF%E3%83%BC

彼が数学オリンピックに出られなかったのは、計算力だけではだめという示唆？（多分違う。。。）

通信回線もミリ波でということですが、あれは距離飛ばないし指向性があるから。。。超ハイパワーで出力しているという設定かな。だとするともはや兵器と呼ぶほうがふさわしい。

和歌山駅近くのオークワはでかい。

こんな大きいの初めて見た。

となりに研修センターまで。本社？（調べる気なし。。。）

JR紀三井寺駅でお帰り。

行きは和歌山駅からがんばって歩いたけれど、帰りは階段がきつくて体力の限界。。。

しかたなく紀三井寺駅から帰りました。でも、ここICOCA使えないよ。。。切符に日付を印刷する機械が改札に置いてありました。こんなの初めて使った。

紀三井寺に行ってきた（その２）。

大きな観音様像のある建物は展望台にもなってますよ。すごい眺め。マリーナシティーのところにあるあの塔って何だろう？左側ね。右側は普通の塔っぽいけど。

桜がきれい。