2009年12月29日火曜日

ファインマンに学ぶ科学的良心

年末は,ファインマンに学ぶ科学的良心.

リチャード P. ファインマンの自伝「ご冗談でしょう、ファインマンさん」の「カーゴ・カルト・サイエンス」.これは彼のカリフォルニア工科大学1974年卒業式式辞を文書化したものですが,非常に素晴らしいことが書かれています.ということで一部抜粋.


 今あげたこういう教育上の、または心理学上の研究は、実は私が「カーゴ・カルト・サイエンス(積み荷信仰式科学)」と呼びたいと思っているえせ科学の例なのです。南洋の島の住民の中には積み荷信仰ともいえるものがある。戦争中軍用機が、たくさんのすばらしい物資を運んできては次々に着陸するのを見てきたこの連中は、今でもまだこれが続いてほしいものだと考えて、妙なことをやっているのです。つまり滑走路らしきものを造り、ヘッドホンみたいな恰好のものを頭につけた男(フライトコントローラーのつもり)をその中に座らせたりして、一心に飛行機が来るのを待っている。形の上では何もかもがちゃんと整い、いかにも昔通りの姿が再現されたかのように見えます。
 ところが全然その効果はなく、期待する飛行機はいつまで待ってもやってきません。このようなことを私は「カーゴ・カルト・サイエンス」と呼ぶのです。つまりこのえせ科学は研究の一応の法則と形式に完全に従ってはいるが、南洋の孤島に肝心の飛行機がやってこないように、何か一番大切な本質がぽかっとぬけているのです。



 諸君に第一に気をつけてほしいのは、決して自分で自分を欺かぬということです。己れというものは一番だましやすいものですから、くれぐれも気をつけていただきたい。自分さえだまされなければ、他の科学者たちをだまさずにいることは割にやさしいことす。その後はただ普通に正直にしていればいいのです。
 もう一つ、これは科学にとってはさほど重大なことではないが、私が信じていることを申し上げたい。それは諸君が科学者として話をしているとき、たとえ相手が素人であっても決してでたらめを言ってはならないということです。(中略)私が言わんとしていることは、嘘を言う言わないではなく、科学者として行動しているときは、あくまでも誠実に、何ものもいとわず誠意を尽して、諸君の説に誤りがあるかもしれないことを示すべきだということです。これこそ科学者同士の間ではもちろんのこと、普通の人たちに対するわれわれ科学者の責任であると私は考えます。



 あるとき私はラジオに出演する友人と話をしていて、少なからずびっくりしたことがあります。彼は宇宙論と天文学をやっている男でしたが、この研究がどのようにして実際に応用できるかを、いかに説明すればいいかというので、頭を悩ませていました。私がつい「そうさなあ。応用できることなんか何もないんじゃないか」と言いますと、彼は「うん、実はその通りなんだ。しかしそんなことではこの種の研究に研究費が出してもらえなくなるからな」と言いました。私はこれを一種の不正直さだと考えるのです。もし科学者として立つなら、素人に自分の研究していることをよく説明し、それでも支持してくれないのなら、それは向こうの決めることなのだからしかたがないのではないか。


最近(というかちょっと前)話題になっている「事業仕分け」.それで上の話を思い出しました.
内なる世界にいれば,どうしても盲目的になりがちな気がしますが,そういう人も他の人に対し自分が何をやっているか懇切丁寧に説明しなければならないのではないかと感じました.

2009年12月26日土曜日

カオス#30 リアプノフ指数

リアプノフ指数(Lyapunov exponent)

リアプノフ指数とは,力学系において接近した軌道が離れていく度合いを表す量を表します.
1次元写像に関して,この値は非常に簡単に求めることが出来ます.
具体的な式は以下の通り.



例:

テント写像

テント写像T(x)は以下で定義されます.



パラメータaを0から2まで変化させたときのリアプノフ指数を計算してみます.



Fig.1

0≦a<1 ⇒ λ>0
1≦a≦2 ⇒ λ<0
となっている所がポイントです.

ここで,テント写像の分岐図を見てみましょう.


Fig.2

分岐図からは,a>1のとき,カオスになっていることが分かります.
これより,(強引だけど)λ>0の時,写像はカオス性を持つと予想出来ます.

ロジスティック写像

ロジスティック写像L(x)=ax(1-x)において,aを変化させた時のリアプノフ指数を求めてみます.


Fig.3 0≦a≦4
3.5≦a≦4が特にゴチャゴチャしているので,拡大してみます.



Fig.4 3.5≦a≦4

所々λ<0となっています.

さて,ロジスティック写像の分岐図はFig.5です.


Fig.5

これらを比較してみると,λ<0の時,分岐図はカオスではありません,ある周期の周期軌道に落ち込んでいます.

そんなわけで,リアプノフ指数は写像がカオス性を持つかどうかを判断するのに非常に役立つといえるでしょう.

2009年12月21日月曜日

カオス #29 クリスマスツリーの製作

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

クリスマスが今年もやって来る
悲しかった出来事を 消し去るように
さあ パジャマを脱いだら 出かけよう
少しずつ白くなる 街路樹を駆け抜けて

クリスマスがもうじきやって来る
うれしさをかくせない 犬や猫まで
もみの木に灯る明かりをうけて
いつもよりやさしそうな パパの目が笑ってる

クリスマスは誰にもやって来る
もしひとりぼっちでも 淋しがらずに
心に住むサンタに呼びかけて
幼い頃の夢を 思い出してごらんよ

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

すてきなホリデイ/竹内まりや

やってきましたクリスマス.
今年は「カオスで作るクリスマスツリー」です.

■木
ツリーは以前作ったシダのパラメータを変えることで表現してみます.

初期座標を(x0=0, y0=0)とし,以下に示す4つの座標変換の内一つを選択し,反復適用していきます.Pは選択される確率です.

xn+1=-0.08xn
yn+1=0.72yn-1
P1=1%
xn+1=-0.16xn-0.32yn
yn+1=0.4xn-0.18yn+1.6
P2=7%
xn+1=0.16xn+0.32yn
yn+1=0.4xn-0.18yn+1.6
P3=7%
xn+1=0.86xn
yn+1=0.91yn+0.65
P4=85%


Fig.1

■星,雪
これらはStar juliaというフラクタルを用いてみます(詳細は後日).

・星
m=2, n=2

Fig.2

・雪
m=1, n=1

Fig.3
このパラメータだと三角形が出来あがるので,60°ずらしたものを重ねています.

■クッキー
Gingerbread man(カオス #7 GingerbreadMan参照)が正に其のものです.


Fig.4

■飾り
いいネタが無いので,良く分からない飾りを幾つか.
グモウスキー・ミラ写像というものを使いました(詳細は後日).


Fig.5
a=0
m=0.302


Fig.6
a=0
m=-0.36


Fig.7
a=0
m=-0.237

■仕上げ
これらをいい感じに組み合わせます.


Fig.8

完成.カオスとともにクリスマス・イブ…それは予測不可能の中にある秩序を見出す喜び…

2009年12月20日日曜日

数学 #13 数値計算 #6 Lotka-Volterra方程式

Lotka-Volterra(ロトカ=ヴォルテラ)方程式

Lotka-Volterra方程式とは,捕食者と被食者の増減関係をモデル化した非線形微分方程式です.

dx/dt=x(α-βy)
dy/dt=-y(γ-δx)

x:被食者
y:捕食者
t:時間の個体数
α, β, γ, δ:正の実数パラメータ

被食者の増殖速度
dx/dt=αx-βxy
αxは,被食者が多ければ多いほど,その増加速度は大きいことを意味します,つまりαは増加に関するパラメータ.
-βxyは,被食者,捕食者が多い程,被食者の個体数が減少することを意味します(被食者が食べられる量,捕食者が食べる量はそれぞれの個体数に比例),つまりβは減少に関するパラメータ.

捕食者の増殖速度
dy/dt=δxy-γy
δxyは,被食者,捕食者が多い程,捕食者の個体数が増加することを意味します(↑の-βxyと対になっている),つまりδは増加に関するパラメータ.
-γyは,捕食者の自然現象の速度を表します,つまりγは減少に関するパラメータ.

今回はルンゲ=クッタ法で数値的に解いてみましょう.
被食者-赤
捕食者-青


Fig.1
初期値:
x0=1
y0=1
パラメータ:
α=0.9
β=0.5
γ=0.75
δ=0.25


Fig.2
初期値:
x0=1
y0=1
パラメータ:
α=1.5
β=1.5
γ=0.1
δ=0.5


Fig.3
初期値:3
x0=γ/δ=3
y0=α/β=1.8
パラメータ:
α=0.9
β=0.5
γ=0.75
δ=0.25

Fig.1を見ると,周期解になっていることが分かります.また,被食者,捕食者の位相が若干ずれていることが分かります.これは,

魚が増える→サメが食べることの出来る量が増えるので,サメは増え,魚は急激に減る→サメが食べることの出来る魚の量が減ったので,サメは自然減少する他無い→サメが減ったので,魚は増える→…

というサイクルを表しているのです.
※ここで,
魚=被食者
サメ=捕食者

Fig.2では,周期がFig.1のおよそ2倍になっています.また,個体数のピーク時は,被食者よりも捕食者の方が多いことが分かります.

Fig.3では,未来永劫 個体数が変わりません.これは,初期値がx0=γ/δ,y0=α/βでは,

dx/dt=0, dy/dt=0

となるためです.

現実におこる現象は人智を超えた複雑怪奇なものですが,こんな簡単な式によって決まっているのかもしれません,しかしそこには無秩序と秩序の狭間にあるカオスが渦巻いているのでしょう.

2009年12月19日土曜日

TopCoder Member SRM 455

Member SRM455(12/18 1:00~3:00)の報告.

・MinimalDifference(DIV2 Easy)
ある大きさの長方形の各セルに蜘蛛を設置,各蜘蛛には1ターンで東西南北の4方向の内,どの方向に進めるかが与えられる.1ターン過ぎた後,蜘蛛がいないセルはいくつあるか.

Result:Passed System Test(221.87p)

・WordsGame(DIV2 Medium)
寝落ち….
Result:Opened(0p)

・NumbersAndMatches(DIV2 Hard)
寝落ち….
Result:Opened(0p)

・Challenges
寝落ち….

Rating
1153->1155

2pしか上がりませんでした.

次回は12/23…

2009年12月15日火曜日

カオス #28 Clifford Attractor

再びClifford Attractor.少しサイケな世界をお楽しみください.


Fig.1 a=1.3, b=1.7, c=0.5, d=1.4


Fig.2 a=-1.5, b=-3.0, c=1.5, d=1.5


Fig.3 a=1.4, b=-1.4, c=0.5, d=0.8

2009年12月13日日曜日

ファインマン伝説

・学生時代,ガリ勉野郎の部屋のドアを外して地下室に隠す
・マンハッタン計画中に機密の入った金庫を開けまくる
・余りに腕が良かったので,番号を忘れた仲間から依頼されることも
・検問にむかついたので抜け穴から出る→入るを繰り返してお説教を食らう
・と思ったら今度は検閲にむかついて手紙をジグソーパズル化する.当然怒られる
・当時から物理学の大家であったボーアに「何とち狂ったこと言ってるんだ」と発言
・打楽器ボンゴの名人
・ブラジルにいた頃「明日は名物のパレードが見られますよ」と言われたがそのパレードに参加していた
・精神科医が嫌いなのでいたずらしたら兵役を免れた
・絵も得意.個展を開いたこともある
・たまたま喧嘩を目撃し,皿が飛んでるのをみて皿の運動を計算する
・その理論を元にしてノーベル物理学賞を受賞
・受賞の知らせに「何て無礼な奴なんだ,今深夜じゃないか」
・「辞退する方法はないのか?」「辞退したらもっとややこしいことになりますよ」→諦めて受け取る
・初来日の時,日本に来たんだし日本の旅館に泊まりたい!というので無理を言って宿泊場所を変更
・作法が分からないのでとりあえずヒャホゥ!と風呂に入ったら先客の湯川博士に怒られた
・別の来日時は「不敗魔」と自己紹介した

ノーベル物理学賞受賞者ファインマン伝説より,一部修正.

2009年12月7日月曜日

TopCoder SRM 454

SRM454(12/6 2:00~4:00)の報告.

・MinimalDifference(DIV2 Easy)

正整数A,B,Cが与えられる.
A~Bの数の中で,Cとの「差」の絶対値が一番小さく,その中でも最も小さい値を返しなされ.
xとyとの「差」の絶対値は,例えばこんな感じ.
x=1234, y=3443
1+2+3+4=10
3+4+4+3=14
|14-10|=4

簡単な四則演算で実装可能.計算オーダーも少ないのでちゃっちゃっと提出.
Result:Passed System Test(229.46p)

・WordsGame(DIV2 Medium)
あるNxN文字行列(String[]),N文字からなる単語(String)がそれぞれ与えられる.
文字行列には,「ある2つの行or列を入れ替える」という操作が出来る.
縦or横方向に単語が並ぶように文字行列を操作した場合,最小の操作回数を求めよ.

「文字行列への操作」を読み間違えていたので,無駄な時間を掛けてしまったです.
提出したソースは本質的には問題無かったはずですが,System Testで落とされました.

Result:Failed System Test(0p)

・NumbersAndMatches(DIV2 Hard)

ある整数(long N)をマッチ棒で表す.この際最大K本のマッチ棒を動かせたとすると,何種類の整数を表すことが出来るか.但し,動かした後の整数の桁数は動かす前と同じであるとする.

WordsGameで時間(と気合)を取られてしまったので手付かず.

Result:Opened(0p)

・Challenges

WordsGameにてサンプルすら通らないと思われるソースがあったので,サンプルを落とし込んで撃墜.
Result:Challenge Successed(50p)

Rating
1163->1153

3回連続でレーティングが下がっています.500点問題も解かないと,レーティング上昇は厳しいです.

2009年12月6日日曜日

専門書を読む #16

プログラミングのための線形代数
~p147

TopCoder 練習

12/6
ToolsBox(SRM453.5 DIV2 Easy)
MazeMaker(SRM453.5 DIV2 Medium)

MazeMakerはProblem Archiveを見ながら.

カオス #27 畑政義写像

以前,大量のアトラクタを紹介した畑政義写像でしたが,今回は,それらをアニメーション化した物を.
数々のアトラクタのパラメータを連続的に変化させることで,非常に面白い動きを見ることが出来ます.

2009年12月3日木曜日

数学 #12 数値解析 #5 van der Pol方程式

van der Pol方程式

van der Pol方程式とは,次の微分方程式です.

y''+ε(y2-1)y'+y=0

Wikipediaより,
「ファン・デル・ポール振動子とは、非線形の減衰を受けた非保存系の振動子である。」
だそうです.

今回はこれを解いてみましょう.

数値計算で求める場合,色々な手法がありますが,ここでは,
v=y'
とおき,

y'=v
v'=ε(1-y2)v-y

とすることで,2変数にしてみましょう.
さらに数値計算手法として,3次のテイラー級数を用いてみます.
即ち,

yn+1=yn+hy'n+(h2/2)y''n+(h3/6)y'''n
vn+1=vn+hv'n+(h2/2)v''n+(h3/6)v'''n

y'=vやy''=v'は上にありますし,v''=y'''やv'''は気合で求められます(少しのdと公式 そして私はあなたと微分を覚えた).

初期条件の設定
基本的に,数値計算で微分方程式を解くときは初期条件が必要です.
具体的には最初に
y0と,v0=y'0
を決めれば良いのです.

実際の計算例

以下のグラフでは,横軸にy,縦軸にv=y'をとっています.
Fig.1にベクトル場((y, v)での傾き),Fig.2に数値計算で求めた解を示します.


Fig.1


Fig.2

パラメータ:
ε=0.3
h(刻み幅)=0.01
初期値
赤(y0, v0)=(2.00092238555422, 0)
青(y0, v0)=(0.5, 0)
緑(y0, v0)=(2.5, -2.5)

さて,計算結果から以下の事が分かります.
赤の軌道…周期解になっている(以後ずっと回り続ける).
青の軌道…閉曲線(赤軌道)の“内側”を初期値としているが,赤軌道に落ち込んでいる.
緑の軌道…閉曲線(赤軌道)の“外側”を初期値としているが,赤軌道に落ち込んでいる.

このようにvan der Pol方程式の解は,十分な時間が経つとこの閉曲線上(赤軌道)を回っている点の運動とすることが出来ます.このようなattractorを“limit cycle”というようです.

今回で使った3次のテーラー級数では,刻み幅hを大きくし過ぎると,発散してしまうようなので注意.

2009年12月2日水曜日

専門書を読む #15

最近溜めていたので更新.

プログラミングのための線形代数
~p57