スラックラインの数理 2 張力と沈み込み量

かなり間が空きましたが、『スラックラインの数理』2回目です。今回は、実際にスラックラインを張るときに役立つかもしれない話題、ラインの張力と沈み込み量の関係を見てみたいと思います。沈み込み量とは、ラインに人が乗ったときにラインと足の接点が元の高さから下がる長さです。足とラインの接点は、この沈み込み量を半径とする円周上で揺れるので、この量が変わると乗り心地が全然異なります。一般的に、沈み込み量が大きいほど慣れるのに時間がかかりますが、小さい場合もそれはそれで特有の難しさがあるのでどちらが難しいというわけではなく、とにかく感覚が全然違うとしか言いようがありません。

スラックライン業界ではすでに知られた式があって、例えばスラックラインの歩き方さんのホームページではその計算フォームも公開されています。今回は、これよりももう少し踏み込んだ分析をやってみますが、とりあえずはこの式を導きましょう。

数式が多いのでMathJaXを使用します。古いブラウザでは読めないかもしれません。
キャプチャ 続きを読む »

素数の下一桁の分布について

「素数の出方はランダムではなかった。1億個調べて浮かんだ奇妙な数」という記事があったので、Mathematicaで確認してみました。まず100万番目までの素数を求め、素数とその次の素数の間で末尾がどう変化するかの遷移行列を作ります。一桁の素数は末尾が特異なので、外しました。

n = 1000000;
primes = Table[Prime[i], {i, 5, n}];
pairs = Table[{Prime[i], Prime[i + 1]}, {i, 5, n}];
mat = Map[Normalize[#, Total] &,
Table[Cases[Mod[pairs, 10], {i, j}] //
Length, {i, {1, 3, 7, 9}}, {j, {1, 3, 7, 9}}] // N];
mat // MatrixForm

MatrixPlot[mat, FrameTicks -> {{{1, 1}, {2, 3}, {3, 7}, {4, 9}}, {{1, 1}, {2, 3}, {3,7}, {4, 9}}}]

sos1

sos4-2

続きを読む »

Mathematicaで重力波データを解析

重力波検出のニュースが話題ですが、その波形生データが公開されていたのでMathematica で分析してみました。取り扱うのは4096Hzでサンプリングされた4096秒分×測定値2箇所分(Hanford, Livingston)のHDF5形式データです。(他に16384Hzのデータ、32秒分のデータがあるので、測定地点×サンプリング周波数×測定時間で計8通りのデータが載っています。)

分析するには、理論上どんな波形が検出されるかを知っておいたほうが有利です。今回は都合良く、2つのブラックホールが衝突した時の理論曲線が公開されています

sr = 4096;
template =   Partition[ToExpression[StringReplace[StringSplit[Import[“https://losc.ligo.org/s/events/GW150914/GW150914_4_NR_waveform.txt”]], “e” -> “*^”]], 2];
template2 = #[[2]] & /@ template;
ListLinePlot[template2]

ligo09

これが観測データに含まれているはずなので、パターンマッチングで見つければ良いわけです。波形でパターンマッチングと言えば相互相関関数ですかね。相互相関関数とは、観測データからパターンと同じ秒数分を取り出して、かけ算して積分(つまり内積)した値を並べた関数です。パターンの長さをTp秒、観測データの長さをT秒とすると、まず観測データの0~Tp秒分を取り出して積分を計算し、これを少しずつずらしながら最後はT-Tp~T秒分の積分値まで計算します。するとTからTp秒分だけ少なくなった一連の積分値が得られるのですが、これを関数と見なすわけです。波形の大きさを正規化しておけば、パターンと完全に一致している箇所では1に、似ていないほど0に近くなり、パターンと完全に正負が逆転した箇所では-1になります。観測データとパターンを相互相関関数でつきあわせてみましょう。

dataH = Import[“https://losc.ligo.org/s/events/GW150914/H-H1_LOSC_4_V1-1126257414-4096.hdf5”, {“Datasets”, “/strain/Strain”}];
dataL = Import[“https://losc.ligo.org/s/events/GW150914/L-L1_LOSC_4_V1-1126257414-4096.hdf5”, {“Datasets”, “/strain/Strain”}];

corH = ListCorrelate[template2, dataH]/Sqrt[MovingAverage[dataH^2, Length[template2]] Length[template2] Total[template2^2]];
posH = Position[corH, Max[corH]][[1, 1]]
corL = ListCorrelate[template2, dataL]/Sqrt[MovingAverage[dataL^2, Length[template2]] Length[template2] Total[template2^2]];
posL = Position[corL, Min[corL]][[1, 1]]

ListLinePlot[Downsample[corH, 10], PlotRange->All
ListLinePlot[Downsample[corL, 10], PlotRange->All]

ligo01ligo02

さすがに、生データからは何も読み取れませんでした。(最初にピークがありますが、端点なので処理上異常値が出ているようです。)そこで観測データに何らかの処理をしなければならないわけですが、今回必要なことはある一瞬の変化を際立たせることで、逆に言えばデータ全体にまんべんなく入っているような動きは抑え込まなくてはなりません。これにはホワイトニングという手法を使います。一番簡単なのは、観測データをフーリエ変換して各周波数の成分(複素数)を求め、その絶対値で割り算し、フーリエ逆変換で元に戻すという方法です。これにより、観測データ全体に入っているような周波数ほど抑えられ、局地的にしか入っていないデータが強調されることになります。

dataHwhitened =   InverseFourier[Fourier[dataH]/(Fourier[dataH] // Abs)] // Chop;
dataLwhitened =   InverseFourier[Fourier[dataL]/(Fourier[dataL] // Abs)] // Chop;

これで、観測データの見かけはホワイトノイズになりました。これをもしフーリエ変換したら、全周波数で同じ強度のパワースペクトルが得られるはずです。このデータで相互相関関数を作ってみます。先ほどとほぼ同じコードです。

corH = ListCorrelate[template2, dataHwhitened]/Sqrt[MovingAverage[dataHwhitened^2, Length[template2]] Length[template2] Total[template2^2]];
posH = Position[corH, Max[corH]][[1, 1]]
corL = ListCorrelate[template2, dataLwhitened]/Sqrt[MovingAverage[dataLwhitened^2, Length[template2]] Length[template2] Total[template2^2]];
posL = Position[corL, Min[corL]][[1, 1]]

ListLinePlot[Downsample[corH, 10], PlotRange->All]
ListLinePlot[Downsample[corL, 10], PlotRange->All]

ligo03ligo04

これで、明確に検出されました。実は2つのデータでピーク位置がデータ 29個分(7ms)だけずれているのですが、これは観測地が離れているために重力波の到達時刻が異なったのが原因です。観測所同士は約10ミリ光秒だけ離れているので、(重力波の進行速度が光と同じだと仮定すると)重力波が斜めに入ってきたことがわかります。これを元に、ある程度重力波源の位置を絞れるのでしょう。

この相関関数のヒストグラムを描くと、正規分布にかなり近いように見えます。(±1の範囲しか動かないので厳密に正規分布ではありません。)

ligoligo7

正規分布だと仮定すると、信号位置のピークは10.93σと8.16σとなり、かなり強く検出されたように思われます。2箇所の観測データを突き合わせて、例えば0.05秒以内に2つのデータが6σ以上を閾値とすると、誤検知はもっと抑えられるのではないでしょうか。

(上の時系列グラフはダウンサンプリングしているので、本当のピーク値はもう少し大きい)

検出された付近の波形を見てみます。ピーク部分の観測データを、正規化してパターンと重ね合わせてみました。また、ピークが負だった方のデータはパターンと逆転しているということですから、正負を逆転させています。

ListLinePlot[
Normalize /@ {template2[[Floor[Length[template2] 8/10] ;; Length[template2]]], dataHwhitened[[posH + Floor[Length[template2] 8/10] ;; posH + Length[template2]]],

-dataLwhitened[[posL + Floor[Length[template2] 8/10] ;; posL + Length[template2]]]}, PlotRange -> All]

 

ligo07高周波のノイズがまだ残っていて、わかりにくいですね。データから10~280Hzの部分だけを取り出してみます。

dataHfiltered =BandpassFilter[dataHwhitened, {10, 280} 2 Pi, 300, HammingWindow,SampleRate -> sr];
dataLfiltered =BandpassFilter[dataLwhitened, {10, 280} 2 Pi, 300, HammingWindow,SampleRate -> sr];
ListLinePlot[Normalize /@ {template2[[Floor[Length[template2] 8/10] ;; Length[template2]]],
dataHfiltered[[posH + Floor[Length[template2] 8/10] ;; posH + Length[template2]]],
-dataLfiltered[[posL + Floor[Length[template2] 8/10] ;;posL + Length[template2]]]}, PlotRange -> All]

ligo08

かなり良い感じではないでしょうか。というわけで、家庭のコンピュータでも重力波検出は一応可能でした。ただ実際には、もっと多くのパターンとの照合や、突発的なノイズ(地震、周囲の交通、風など)への対応が必要なのでしょう。

スラックラインの数理 1

最近、スラックラインを練習しています。スラックラインとは、綱渡りのようなもの、というかほとんど綱渡りそのものをスポーツにしたもので、綱渡りとの違いは

  • 綱ではなく、少し伸び縮みする数センチ幅のテープを使う。
  • ヤジロベエのような棒は持たない。
  • 上級者は、ただ歩くだけではなく、跳ねたり飛び乗ったり、トランポリンのように使う。

といったところです。で、これのコツは足を少し曲げて重心を低くすることなんですね。実際、足を曲げると安定します。しかしこれっておかしくないですかね?

いや、だって倒立振子とかでは重心位置が高い方が安定するじゃないですか。登山でも、重いものはザックの上に来るように配置しろと言われるじゃないですか。不思議だ、というわけでモデルを作って分析してみました。
sl 続きを読む »

起震車の揺れをiPhoneで計測して震度を計算

iPhoneは3軸の加速度計を内蔵しています。iPhoneの本体を傾けると画面が縦モードになったり横モードになったりしますが、あれは加速度計で重力方向を検出して、向きを検知しているわけです。しかしせっかく加速度計があるので、もっといろいろな応用、例えば震度計としても使えるのではないでしょうか。そこで、夏休みの自由研究として、iPhoneの加速度計で地震の震度を計ってみました。そして、気象庁がどうやって震度を決めているのかも調べました。

とは言うものの、地震が来るまで待ってはいられません。運が良いことに、防災行事で起震車に乗る機会があったので、iPhoneのSensorLogというアプリで地震発生から終了までの加速度を計測してみました。震度7の揺れです。

起震車加速度波形
続きを読む »

Twitterデフォルトアイコンの卵を表す方程式が判明(?)

2010年9月以降、Twitterはデフォルトアイコンとして卵の画像を使っています。それ以前は、鳥をかたどったものでした。
larrybird-11引用元
この鳥は、円を組み合わせて作られたようで、数学的に見ても面白い形です。では、現在の卵アイコンはどうなっているのか。これもまた、数学を利用した形ではないだろうかと思い、調べてみました。
twitterlogo
結論から述べます。Twitterのデフォルトロゴは、伊藤忠夫氏の卵形曲線方程式を利用して作られていると推定しました。山本信雄氏のホームページによれば、これは次のような形の式です。係数$a, b$は自由に決めて良いのですが、$b/a=0.78$とすることが推奨だそうです。
\[x=a cos \theta, \hspace{5mm} y=b cos \frac{\theta}{4} sin \theta \hspace{8mm}(-\pi\le \theta \le \pi) \]
続きを読む »

スパゲッティの乾麺を折るとなぜ3つ以上に分裂するのか(1)

 最近、スパゲッティの乾麺を折ると3つ以上に分裂することを解明した論文を読んだのですが、いろいろと疑問が湧いてきたので、麺シミュレーションを作ってみました。
pas001
(クリックすると動きます。)

 この動画をtwitterに公開したところ、やんぼりさんと議論になりました。それは、1度目の破断が発生する前の形はどうなるのか、ということに関してです。それをずっと考えていたのですが、何とか解らしきものが形になったので、公開したいと思います。
続きを読む »

モンティー・ホール問題をMathematicaで解く

Mathematicaには、Ver8.0頃から強力な確率・統計計算機能が搭載されています。しかし、あまり利用されているという話を聞きません。そこで、とある確率上の難問を取り上げつつ、この機能を使ってみたいと思います。

モンティー・ホール問題とは

数学の分野で、定期的に話題になる確率の問題に、「モンティー・ホール問題」があります。

テレビ番組で、司会者のM氏(モンティー・ホール氏)が3つの扉を提示する。出場者は、そのうち1つだけを開けて良い。いずれかの扉の裏側に景品があり、当たった場合はそれを貰える。

出場者が左の扉を選んだところ、M氏は真ん中の扉を開け、そこに景品がないことを示した。出場者は、選んだ扉を変更すべきか否か。

続きを読む »

ビットコインの原論文を読む 12節 結論

12 結論

信頼関係を必要としない電子取引システムを提案した。まず、電子署名だけに頼った通貨を考えたが、所有権を強力に制御できるものの、多重使用を防げないことが判明した。この問題を解決するため、各利用者が取引を公開し、演算量証明を利用してその履歴を記録するP2Pネットワークを考案した。もし、善良なノードがCPU能力の大半を占めていれば、攻撃者による偽造は計算論的にほとんど不可能である。本ネットワークの形は決まっておらず、素朴な仕組みなので、堅牢性に優れている。さらに、各ノードは協調せずばらばらに動いて良い。そしてすべての情報について、必ず届かなければならないという場所を設ける必要は無く、またベストエフォートで広がれば良い。つまり、ノードはすべて平等に扱える。各ノードはいつネットワークから離れ、再接続しても構わない。その間の演算量証明は他のノードから得られる。ネットワークはCPU能力をもとに投票を実施し、現行のチェーンに対して、正しいブロックを追加し不正なブロックを拒否する。この同意メカニズムによりルールが守られ、ネットワークの維持に貢献した各利用者が報酬を得る。

続きを読む »

ビットコインの原論文を読む 11節 数学的根拠(後編)

引き続きMathJaXを使用します。古いブラウザでは読めないかもしれません。

$p>q$という本システムの前提下では、攻撃者が追いつく確率は、追いつくべきブロック数が増えるほど指数関数的に小さくなる。従って、もし最初に運良く相手との差を詰められなければ、差が急激に開き、追いつく可能性はほとんど無くなるだろう。

支払人による取引改ざんを事実上不可能にするために、受取人はどの程度待てば良いだろうか。攻撃者は、受取人に対してしばらくの間正しく支払われたと思い込ませ、その後に払った額を奪い戻す。受取人は、いずれ不正な取引であるとの警告を受けるが、支払人はその時刻をできるだけ遅らせたい。

続きを読む »

古い記事へ «