Saturday, December 10, 2016

平均寿命の国際比較

平均寿命というものが、個々人にとってどういう意味があるのかははっきりしないですが、平均寿命の国際比較(厚生労働省)というのがあります。

ニュースでは「日本が一位」「長寿日本」という報道の仕方をするので、なんとなく日本の医療のお陰かなと思っているのです。上の厚生省のデータを見ると、確かに日本は順序の上で一位が続いているのですが、実際の意味はあるんでしょうか。

例えば、男の平均は80歳程度なわけですが、自分がどうかというと、70で死ぬかもしれないし、90まで生きるかもしれないと思う。少なくとも±5歳くらいの誤差というか幅があるように感じますね。

そこで他の国を見ると、先進国、少なくともヨーロッパとは大差ないじゃないですか。どうしたニッポン。だいたいが男の場合70後半。誤差の範囲ですよね。まあロシアは低そうですね。不摂生だからでしょうか。アメリカも少し低いようです。

Aの平均と、Bの平均を比べる時に、AとBの差に比べてAやBの分散が大きいと、A、Bの区別はつかないというのが統計学のココロのはず(分散分析)。平均寿命に対する分散が幾つなのかは知りませんが、直感的に誤差の範囲と言って良いんじゃないでしょうか。

ヨーロッパのある国では、病院が混雑していて、救急車で運ばれないかぎりお医者さんにj見てもらえないだの、待合室で待っている間に死んだだのと悪口を聞きますが、その国の寿命も特に低くないですね。不思議です。イギリスの話として聞きましたが、ブリティッシュジョークなのかな。

もう一つは、寿命ってのは人種と関係ないように思えますね。体の大きさや肌の色は違うのに、寿命ってのは人間の仕組みが持ってる根源的なものなのですかね。


Monday, October 31, 2016

Python3とPython2違いメモ

join
Python2では";".join(['a','b'])という奇怪な文法だったが、Python3ではは以下のようになっている。
 
>>> str.join(',', ['a','b'])
'a,b'
strはbuiltinのクラスなのでインポートする必要はない。

reduce
Python3ではfunctoolsの一部になった。functoolsをインポートする必要がある。
 
>>> from functools import reduce
>>> reduce(lambda x,y:x+y, [1,2,3])
6
map
mapが返すのはiteratorに変わった。listが欲しい場合はlistで囲む。
>>> map(str, [1,2,3])
<map 0x205b20b0="" at="" object="">
>>> list(map(str, [1,2,3]))
['1', '2', '3']
>>> 

多次元配列のsort
やりたいことは2つ目の要素でソートする、などのこと。Python2について昔ここに書いたことがあったのだけど、Python3ではkeyで渡すように変わっている。関数の書き方も変わったようだ。
>>> a = [[1,2],[3,0]]
>>> a.sort()
>>> a
[[1, 2], [3, 0]]
>>> a.sort(key=lambda x: x[1])
>>> a
[[3, 0], [1, 2]]
数の型
$ python2.7
Python 2.7.10 (default, Jul 30 2016, 18:31:42) 
>>> type(1)

>>> type(1.2)

>>> 
$ python3.5
Python 3.5.2 (default, Aug 16 2016, 05:35:40) 
>>> type(1)

>>> type(1.5)


Tuesday, July 05, 2016

matplotlib入門:アニメーション。バブルソート

バブルソートって、なんだっけ?と思ってwikipediaを見てて、右下の方のアニメーションがいい感じで、自分でも作ってみたくなった。

matplotlibでやってみようと思った。いろいろ調べてこんな感じになった。
100個の乱数を並べてゆく。バブルソートは一番右(か左)から席が確定してゆくので面白いことになる。



ちなみに肝心のバブルソートのアルゴリズムはgen_points()の中の4行くらいの箇所。

できたmp4のアニメをyoutubeにアップロードした。


Saturday, June 18, 2016

Seabornでaxを設定する

MatplotlibだとダサいのでSeabornを使うという話があるのだけど、細かな設定をどうするのかわからんなあと思っていた。

どうやるかというと、matplotlibのplt.subplotsのaxを使って、axをseabornに渡せtば良いのであった。

以下の例ではy軸をカンマ区切りに設定したaxをpointplotの際にax=axとして渡している。
(df.passengerをわざと大きな値にしている)

結果は以下。

Saturday, May 28, 2016

人工知能とデバッグ

少し前に、NHKスペシャルの「天使か悪魔か羽生善治、人工知能を探る」というのを見たのだけど、将棋の羽生がAlphaGoのDeepMindに行ってDemiss Hassabisにインタビューする内容だった。

この番組の1シーンで、Lee SedolがAlphaGoに一勝したときの記者会見で、記者がHassabisに「碁の場合はマシンが暴走した、と言って許されるかもしれないが、ヘルスケアでそれで許されるのか」という趣旨の質問をしていた。これはなかなか良い質問だと思った。DeepMindはヘルスケアをやって行くと言っているのだ。Hassabisの答えは、AlphaGoは完成品ではなく、ヘルスケアの場合は完全にテストしたものを製品とする、という優等生的なものだった。

このやり取りについて考えていたのだが、人工知能で非常に正解率が高い判断をするものの、当然失敗がある。結果を人間が直接見るのが前提なケースは良いだろう。しかしこれを取り込むシステム全体はどうなんだろう?番組の中では「良心回路」を持った人工知能の研究が出てきていたが、私が言いたいのはもっと単純なエラーについてだ。道が右に曲がっているという判断に基づき、自動運転車が右に曲がろうとする。この成果率は非常に高い。しかし100%ではない。間違えるとしたらシステム全体からするとバグだ。複数のロジックを用いて確率を下げる、映画2001年のHALみたいに。でもエラーはゼロとは言えない。どう判断するのだろう。人間の正解率と同等なら良しとするのだろうか。

エラーからの回復方法は、テストを完全にするとかどうという以前に、どうやってシステムをつくるかの思想の問題になる。実際にはケースごとにいろいろ逃げ道あるのだろうが。記者がこのつもりで質問したのかどうかは定かではないが、私も聞いてみたい。

Monday, May 23, 2016

Raspberry Pi FreeBSD。ScanSnapのサーバにする

Saneという名前は聞いたことがあったが、ScanSnapのデータが読めるとは思わなかった。

まず、sane-backendsをインストールする。

pkgが最後に使い方についていろいろ言ってくるが、まず以下の出力を見るように言われる。

root@rpi-b:/home/freebsd # usbconfig list
(略)ugen0.5: at usbus0, cfg=0 md=HOST spd=HIGH (480Mbps) pwr=ON (98mA)




この0.5ということを調べて、デバイスの情報を得る。
root@rpi-b:/home/freebsd # usbconfig -d 0.5 dump_device_desc
ugen0.5: at usbus0, cfg=0 md=HOST spd=HIGH (480Mbps) pwr=ON (98mA)

(略)
 idVendor = 0x04c5
 idProduct = 0x11a2


そしてこの情報をsaned.conf というところに入れて、

root@rpi-b:/home/freebsd # cat /usr/local/etc/devd/saned.conf
notify 100 {
       match "system" "USB";
       match "subsystem" "INTERFACE";
       match "type" "ATTACH";
       match "cdev" "ugen[0-9].[0-9]";
       match "vendor" "0x04c5";
       match "product" "0x11a2";

       action "chown -L cups:saned /dev/\$cdev && chmod -L 660 /dev/$cdev";
};


devdをリスタートする。

root@rpi-b:/home/freebsd # /etc/rc.d/devd restart

すると、scanimageでスキャナが認識できる。
root@rpi-b:/home/freebsd # scanimage -L
device `fujitsu:ScanSnap S1500:14489' is a FUJITSU ScanSnap S1500 scanner

GUIを使うつもりはないので、コマンドラインで操作する。
テストの紙をスキャナに一枚入れて、以下のコマンドを実行すると、あら不思議。
root@rpi-b:/home/freebsd # scanimage --device-name='fujitsu:ScanSnap S1500:14489' --format tiff > test.tiff

紙が吸い込まれて、tiffのイメージができた。

なかなか愉快だ。
まあLinuxでも同じでしょう。

Tuesday, May 03, 2016

RaspberryPi B+, FreeBSDでワイヤレス接続

さて、先日作ったFreeBSD10.3のシステムで、ワイヤレス接続をしてみよう。

実はdmesgの出力を眺めていてwifiがまだだったなあと気づいたのでした。
念の為書いておくと、私がアマゾンでついでにポチったWifiアダプタはGW-USNANO2Aという、おそらくありふれた機種。dmesgにはライセンスを読んでオッケーだと思うならloader.confにlegal.realtek.license_ack=1と書け、とひっそりと書かれていたのでした。その辺のことはurtwnのマニュアルに書いてあります。

マニュアルだけでもわかるのかもしれないけれど、こちらの有識者の設定を参考にしました。やったのは以下。

/boot/loader.confに追加
legal.realtek.license_ack=1

create wlan0を作成
# ifconfig wlan0 create wlandev urtwn0

/etc/rc.confに追加
wlans_urtwn0="wlan0"
ifconfig_wlan0="WPA DHCP"
networkをスキャンして様子を見る
# ifconfig wlan0 up scan
/etc/wpa_supplicant.confに自分のネットワークのSSIDとパスワードを記入
network={
    ssid="your_ssid"
    psk="your_password"
}

# reboot
以上でwifiにつながりました。

Sunday, May 01, 2016

RaspberryPi B+でFreeBSD (続き)

raspberry pi B+でFreeBSDで動作確認は済んだ訳だけれど、その後、どうするか。

まず、前回ダウンロードはしたものの、使わなかったraspbianを試した。これは素晴らしい。RPB+たった512MBのメモリでありながら、頑張ってGUIとRPのロゴを表示した。Python2, Python3のIDLを表示する。えらい。Mathematicaもあるらしい。えらい。

しかし、違う。やろうと思っているのは人が作ったオモチャを使うことではないのだ。

で、まずはBSDにはNanoBSDというツールがあって、組み込み的な組み合わせが作れるらしいということを知った。でRP用のイメージを作ったつもりなのだが、うまく動かなかった。どう動かなかったかというと、画面が真っ暗だったのでどうしょうもなかった。詳細はおぼろげで、よく覚えていない。

次に、crochetというものがあり、これでRP用のイメージが作れるという話らしい。しかし、これは起動時にブートイメージは読めるらしいのだが、ディスク(SDカードね)のマウントで落ちる。mounting from ufs, failed with error 19。FreeBSD本家の情報も見たがよくわからなかった。検索でいろいろ有識者らしい人の古めの情報も見たがダメだった。

と、走馬灯のようにここ2日間のことを思い出す。

現在、最新の10.3Rのリリースノートを見て、この近辺からarm配下RPの公式ブートイメージがあるらしいことを知った。ftp://ftp.freebsd.org/pub/FreeBSD/releases/ISO-IMAGES/10.3/

ここから、以下のものを頂いてddでsdcardを焼いて普通にブートしたのが以下の状態。
FreeBSD-10.3-RELEASE-arm-armv6-RPI-B.img


$ df -h
Filesystem                Size    Used   Avail Capacity  Mounted on
/dev/ufs/rootfs           3.6G    393M    2.9G    12%    /
devfs                     1.0K    1.0K      0B   100%    /dev
/dev/msdosfs/MSDOSBOOT     17M    3.9M     13M    23%    /boot/msdos
tmpfs                      30M    4.0K     30M     0%    /tmp
$ sysctl hw.physmem
hw.physmem: 495087616
$ uname -a

FreeBSD rpi-b 10.3-RELEASE FreeBSD 10.3-RELEASE #0 r297264: Fri Mar 25 08:01:14 UTC 2016     root@releng1.nyi.freebsd.org:/usr/obj/arm.armv6/usr/src/sys/RPI-B  arm


私の懸案であった、既存のUSBドライブの確認まではできた。

root@rpi-b:/home/tsuchiya # dmesg | grep da0
da0 at umass-sim0 bus 0 scbus0 target 0 lun 0
da0: Fixed Direct Access SPC-4 SCSI device
da0: Serial Number 1111114000D1
da0: 40.000MB/s transfers
da0: 1907729MB (3907029158 512 byte sectors)

da0: quirks=0xa


(まだ終わってない)

Sunday, April 24, 2016

RaspberryPi B+ 動作確認 with RaspBSD

ラズパイなるものを買ってみた。動作確認してみよう。

買ったのは、Pi B+とmicroSD カード、USB WifiのPlanexなんたらというもの。Pi B+と一緒に買う人が多いという組み合わせで。

電源はガラケーのmicro usbの充電器を使い、画像出力はテレビのHDMIに出すので、ケーブルはもうある。キーボードはMac Miniのやつを一時的に使う。



OSは、http://raspbsd.org/raspberrypi.htmlからFreeBSD-armv6-11.0-RPI-B-292989.img.gzというのを選択。

Raspberry本家のRaspbian https://www.raspberrypi.org/downloads/raspbian/もダウンロード開始したのだが、BSDの方が先に終わったので、BSDを使うことにした。

Mac MiniのSDカードスロットにカードを入れて、以下でカードを焼きました。




電源を入れると、テレビにブート画面が出た。どうやらwifiは有効になっていないし、このwifiカードを有効にするにはrsuというモジュールを入れてカーネルをビルドしないといけないらしいので、wifiは後日考えることにした。

とりあえず有線のLANケーブルをMac Miniからうばって、接続し、別のMacからSSHで繋ぐとちゃんと入れた。



所要時間1H。




Saturday, April 16, 2016

「原子力政策研究会 100時間の極秘音源」

知らなかったですね。こういう経緯で日本で原子力が始まったとは。無知とは愚かなり。「原子力ムラ」の方々の証言をまとめた本。

日本の原子力導入で、立役者は当時議員だった読売新聞の正力松太郎、それと中曽根康弘。1950年代です。この偉い人たちのリーダーシップが強烈すぎたようだ。しかもあまり専門家の言うことは聞かなかったらしい。まあ、周りの人としては、偉い人のせいにしておけば良いというバイアスがあるのかもしれない。とにかく60年前に始まることなので、当時の責任者はもうこの世にいない。

1950年代の「科学技術」というのはまさに死活問題で、それで中曽根さんは科学技術の推進の為に原子力を推したということだ。現在の人工知能などの科学技術とはちょっとレベルが違うかなと思う。また、こんなにできるかどうか判然としない技術を確立するのは大変だったろうと思う。逆に昔だったからできたことなのかもしれない。

原子力発電所の建設の説明の時にはもちろん安全を謳う必要があった。ある程度商業ベースに乗ってきてからはコスト重視で稼働率を上げることが命題となり、やはり安全性は二の次だった。いやな感じなのは、安全に対する「懸念」が公表しにくいという点。「まだそんなことができてないのか」という批判を恐れて言えないことがあるという話。「原子力ムラ」って、原子力推進に批判的なことを言うと、村八分になる場所らしい。村八分、つまり次の仕事は来なくなる。

活断層と地震に関係があるということが定説となったのは、比較的最近のことなので、多くの原発は断層とは関係なく建っている。(本当だろうかと思って調べたら、どうやら本当で、活断層が注目されるようになったのは阪神淡路震災からで、「活断層学会」ができたのは21世紀のことらしい)

資源の無い国、日本。資源を東南アジアに求めて太平洋戦争を起こした。その反省。核武装を続けるためには、核の平和利用が同時に必要だったアメリカ。

廃棄物であるプルトニウムを燃料として使う核燃料サイクルが究極の夢だったが、これは破綻していて、ゴミはたまる一方となっている。総合的に見ると今まではなんとかやってきたけれど、今後のエネルギー問題は、別の方向で考えた方が良いのではないかと思う。




Saturday, April 02, 2016

python3メモ:unicode

python3になって、エンコード関係でエラーを受けるようになった。いろいろ読むと、python2からpython3で「正しく」なっているようだが、わかりづらい点があるようだ。私もようやく話が見えてきたので、理解を確かめるためにメモしておきたい。

ややこしくなるので、今回はpython2の話は書かないことにする。

まず、文字コードは内部の表現と、入出力が別であって、内部表現はstr型である。これはunicodeでもある。str型はunicodeと同じでpython3では区別がない。入出力できる文字コードはsjisやutf-8などである。これらの型がbytes型である。

つまり、ファイルの内容や、ファイル名などはbytes型としてプログラムに入ってくるので、必要に応じてstr型に変換する必要がある。

変換は、
str -> bytesとするのをencode、
bytes -> strとするのがdecode
となっている。内部表現のstr型が中心だと思えば覚えられるだろう。

要素として、
関数が受け付ける型
ファイルの中身
ファイル名
ロケールとprint()がやってくれる処理
ということが同時に登場するので、混乱するのだと(少なくとも私は)思う。



情報としては、まず、何と言っても公式のページが丁寧だった。

あと、こちらの偉い人が書かれた記事も有用だった。

先ほどこちらの記事を見つけたが、まだちゃんと読んではいない。
http://python-notes.curiousefficiency.org/en/latest/python3/text_file_processing.html

ここで例題をやってみよう。

# 基本編
string = 'アイウエオ'
type(string) # str
print(string) # アイウエオ
ustring = u'アイウエオ' #python3では'u'をつける必要はない。つけても怒られない
type(ustring) # str
print(string==ustring) # True

#変換
bytestr = bytes(string) # TypeErrorで失敗
# エンコードを指定してエンコードする必要
bytestr = bytes(string, 'utf-8')
bytestr2 = string.encode('utf-8')

print(bytestr==bytestr2) # True この2つは同じことのようだ。

print(bytestr) #b'\xe3\x82\xa2\xe3\x82\xa4\xe3\x82\xa6\xe3\x82\xa8\xe3\x82\xaa'

# ファイルへの出力
# デフォルトのコーディングで出力
f=open('Desktop/test_text.txt','w')
# writeに渡すのはstr
f.write(string)
# f.write(bytestr) # TypeError: write() argument must be str, not bytes
f.close()

# 確認してみる。
# $ file Desktop/test_text.txt
# Desktop/test_text.txt: UTF-8 Unicode text, with no line terminators

# 以下のように出力のコーディングを指定できる。
f=open('Desktop/test_text_sjis.txt','w',encoding='sjis')
f.write(string)
f.close()


# str型をprint()に渡して、処理してくれる場合と、できない場合がある。
print(string) #この場合はOK. アイウエオ

# utf-8を手で入力してみる。
T = b'\xe3\x80\x92' #e3 80 92
print(T.decode()) # 郵便マーク 〒

#サマリア文字というのがあるらしいので、ついでにそれも入力してみる。
samaria = b'\xe0\xa0\x80'
print(samaria.decode()) # Samaritan サマリア文字  ࠀ

# デフォルトのコーディングが上記のようにutf-8ならば変換してくれるが、以下のようにロケールがASCIIになっていると、UnicodeEncodeErrorで失敗する。
$ LANG=C
$ python3
Python 3.5.1 (default, Dec 31 2015, 09:25:09)
[GCC 4.2.1 Compatible Apple LLVM 7.0.0 (clang-700.1.76)] on darwin
Type "help", "copyright", "credits" or "license" for more information.

>>> t = b'\xe3\x80\x92' #e3 80 92 郵便マーク
>>> print(t.decode())
Traceback (most recent call last):
  File "", line 1, in
UnicodeEncodeError: 'ascii' codec can't encode character '\u3012' in position 0: ordinal not in range(128)

>>> samaria = b'\xe0\xa0\x80'
>>> print(samaria.decode())
Traceback (most recent call last):
  File "", line 1, in
UnicodeEncodeError: 'ascii' codec can't encode character '\u0800' in position 0: ordinal not in range(128)

# decodeでのコード指定。間違ったコードを指定すると、エラーが返る。
print(T.decode('sjis'))

UnicodeDecodeError: 'shift_jis' codec can't decode byte 0x92 in position 2: incomplete multibyte sequence


まだ他にもあるような気がするけれど、とりあえず今日はここまで。

Saturday, February 06, 2016

texの数式のテスト

このブログで数式を書きたいことが今までそう沢山あったかどうか定かでないのですが、ある記事を見ていたら絵でない数式が書かれていたので、ああそういうことができるんだ、と思って探してみたらこのblogger.comでのやり方を
こちらの方の記事元記事
で知りました。ブログのテンプレートにコードを埋めるんですね。あとは\$でくくってtexの数式を書くだけ。

$\cos^2\theta-\sin^2\theta=\cos 2\theta$

$e^{i\pi} = -1$

$e{^x} = \sum_{n=0}^{\infty} \frac{x^n}{n} $


前に書いた「円周率」の記事の数式のところを直してみた。

Saturday, January 23, 2016

matplotlib入門

これまでグラフを描く時にはRを使ったりpythonでmatplotlibでやったりと色々かじっていたのですが、2ヶ月くらい前にRで2万枚くらいのグラフを描いていた結果、この言語で続けるのはもう無理だという気分になりました。だって変数のスコープがないんだもん。

それでmatplotlibに移ることにしたのだけれど、Rでやってた基本的なことするのにいちいちstackoverflowの記事を読むのが大変だったので、ごく簡単なところをここにメモすることにしました。

情報としては、matplotlib.org本家があるのと、Wes McKinneyの"Python for Data Analysis"を持っています。本の方はMatplotlibに限らずnumpy, pandasなど色々幅広いです。

使い始めてすぐ気になるのは、plt.plotでできることと、subplotのax.plotでできることが似ている一方、違っている点で、結局matplotlibのソースを見て納得しました。pltの方はpyplotモジュールの関数で、axの方はaxesの下の_subplots.AxesSubplotあたりのクラスの関数なのですね。要するに別のもの。

1. 数字のカンマ区切り(1000ごと)。みんなグラフの数字はカンマで区切りましょう!excelでも「数値」を選べばそうなります。カンマをつけてくれないと、いちいちゼロの数を数えるのが面倒です!私は以下のようにFuncFormatterを使ったcommaformを渡すようにしています。

2. 色の変え方。系列ごとにblue, redなどで変える方法はあちこちに出てます。数字によって連続的に変える場合は以下のようにnormalizationと、colormapを使います。

3. colorbarで色の意味を示すのも以下のようにします。例では2の累乗で表示していますがこの例だとデフォルトは10の累乗になります。コメントアウトするとわかります。colorbarの書き方は以下のようにする方法と、コメントアウトしてあるplt.colorbar(sc)
の方法で設定する方法があるようです。

4.  最近もっとも驚いたのは、グラフの上にマウスを持ってきて、'l'(エル)とか'L'と打つと画面上でy軸x軸がログスケールに変わることです。みんな知ってるのかなあ。それはそうとして以下の例ではコメントアウトしましたが、ax.get_xaxis().set_major_locatorでbase=2でx軸をlogに変えることができます。

5. matplotlibはスタイルシートをサポートするようになりました。冒頭コメントアウトしてありますが、plt.style.availableで使えるスタイルシートの一覧が見られます。自分のスタイルシートは
$HOME/.matplotlib/stylelib/foobar.mplstyle
のようなファイルに入れておくと読んでくれます。

 6. matplotlibをいじってみたい気になったので、ソースをgit cloneしてdevelopモードで動くようにしました。ダウンロードしたソースに対してpython setup.py developとすると
/usr/local/lib/python3.5/site-packages/matplotlib.egg-link
というファイルが出来、そこにcloneしたソースへのリンクが書かれています。このローカルコピーで動くことになります。





以下が例。背景が黒いのはdark_backgroundスタイルシートのせい。



Friday, January 22, 2016

壁にぶつかるボールが円周率を作る

2015年2月号の数学セミナーの時枝先生という方の「ちいさい数理みつけた」という記事の中に、壁にぶつかる2つの物体の話が出てくる。

最初、mは静止しており、Mがぶつかってくる。左端には壁があって、床はツルツルしている。mは弾き飛ばされて壁とMに何度もぶつかる。その回数が問題。

質量m, Mの間に $\frac{M}{m} = 100^d,  (d=0,1,2...) $の関係がある時に、衝突の回数が円周率のd桁になるという不思議な話だ。
この記事にはそのことしか書いていないので、自分で確かめよう。

久しぶりに高校の物理を思い出して、m, Mそれぞれの速度をu, vとし、エネルギー保存の法則と、運動量保存の法則を解けば良いことに気づく。衝突後のu,vをそれぞれu',v'とすると、
$mu + Mv = mu' + Mv'$
$\frac{1}{2} mu^2 + \frac{1}{2} Mv^2 = \frac{1}{2} mu'^2 + \frac{1}{2} Mv'^2$
壁に向かう方の符号を正としている。
これを解くと、以下のようになる。

$ \begin{pmatrix}
u' \\
v'
\end {pmatrix}
= \frac {1}{m+M}
\begin{pmatrix}
2 M & m - M \\
M-m & 2 m
\end {pmatrix}
\begin{pmatrix}
u \\
v
\end {pmatrix}
$

mの速さは壁にぶつかっても変わらず、符号だけ変わる。m,Mが衝突すると、mは壁にもぶつかってふたたびMに衝突する。

問題なのは停止の条件で、Mが押し返されて右向きに動き始めてからも、Mの速さよりmが上まわればまた後ろから衝突する。mの方がMより遅くなっても、一回だけ壁にぶつかる。

回数をカウントするプログラムを書いたのが以下。$r = \frac{M}{m} $ とし、$ v = 1$ とした。

結果、以下のように円周率に従っている。

0 3
1 31
2 314
3 3141
4 31415
5 314159
6 3141592
7 31415926
どこにも円周率らしきものがないのに、結果として円周率になっているのが不思議である。

停止の条件が自分ではよくわからず、時枝先生の記事にはロシア人らしき人の名前が引用されているだけだったので、元ネタの所在も分からなかったのだけれど、「衝突 円周率」で調べているうちにNY Timesの記事が見つかって、そのアニメを見て停止の条件がわかった。というかこのアニメ、前に見たような気がする。


-----
(2/6/2016)
2つ後の記事に書いたがブログにtexで数式が書けるようになったので、式のところを直した。

Saturday, January 16, 2016

山の高さの測り方

太陽系最大の火山は火星のOlympus Monsで、こちらを読むと以下のように書いてあります。
Olympus Mons is a shield volcano 624 km (374 mi) in diameter (approximately the same size as the state of Arizona), 25 km (16 mi) high, and is rimmed by a 6 km (4 mi) high scarp.
 高さが25kmあって、広がりは624km。だいたいアリゾナ州と同じ大きさというのですが、その例えだと大きいのか小さいのかわからないですね。しかし日本の本州の幅が200kmだということを考えれば、尋常でない大きさだとわかります。高さが25kmというのも想像を絶しています。

この記事のもうちょっと下を見ると、地球における最大の火山が紹介されていて、それはハワイのMauna Loaだというのです。
To compare, the largest volcano on Earth is Mauna Loa. Mauna Loa is a shield volcano 10 km (6.3 mi) high and 120 km (75 mi) across. 
 高さが10kmというわけですが、ちょっと考えればこれはエベレストより高いことがわかる。なぜそうなのかを考えると、通常我々が知っている山の高さは海抜であって、火星には海がないので、測り方が違っているのであろうということにすぐ気づくわけです。

Mauna Loaが10kmは海底からからだろうなと思って探すと、こういうのが見つかります。理屈としては、この火山は海の底から続いているので、底から測るのが正しいというわけです。で、その数字は以下のように17km。
Thus, the total relief of Mauna Loa, from its true base to its summit, is about 17,170 m (56,000 ft).
10kmという話はどうなってしまったのでしょう。不思議なものだと思って、さらにここの「世界で一番高い山」という記事ですが、

1. エベレスト 海抜8850m
2. Mauna Kea 海底から10km
3. Chimborazo 地球の中心から6384km

と書いてある。 Mauna Loaの立場はどうなってしまったのだろうと思いますが、地球の中心から、っていうのも面白かった。