みなさん、あけましておめでとうございます。
　新年開始とともに、ニコニコに動画を投稿しました。

　おおむね好評だったようで嬉しい限りです。コメントをくださった方、ありがとうございます。
　もっと詳しく説明しろとのコメントがいくつか寄せられていたので、解説します。2と9の文字は、両方とも4つの点とその周りを公転する4点、更にそれを回る4点を元に作られます。以下、それぞれを恒星、惑星、衛星と呼ぶことにします。4衛星を2組にわけ2本の直線を作り、その交点に文字を描かせます。
　それぞれの点の座標は、時刻をtとして次のように表しました。

p1=(a1*cos(2(2t-e1)),a1*sin(2(2t-e1)))+p5
p2=(a2*cos(2(2t-e2)),a2*sin(2(2t-e2)))+p6
p3=(a3*cos(2(2t-e3)),a3*sin(2(2t-e3)))+p7
p4=(a4*cos(2(2t-e4)),a4*sin(2(2t-e4)))+p8
p5=(a5*cos(3(t-e5)),a5*sin(3(t-e5)))+(c5,d5)
p6=(a6*cos(3(t-e6)),a6*sin(3(t-e6)))+(c6,d6)
p7=(a7*cos(t-e7),a7*sin(t-e7))+(c7,d7)
p8=(a8*cos(t-e8),a8*sin(t-e8))+(c8,d8)

(c5,d5)～(c8,d8)が4恒星の座標、p5～p8が4惑星の座標、p1～p4が4
衛星の座標となります。a1～a8は公転半径の大きさ、e1～e8は公転での位相を決定します。またp1とp2を通る直線、p3とp4を通る直線を結びその交点を求め、これを点pとします。pの座標は、複雑な形になりますが方程式を解けば式として表すことが可能です。a*,c*,d*,e*の合計24変数を調整し、pが文字を描くようにすれば良いわけです。
　なおそれぞれの式のtには1,2,3の係数が掛かっていますが、これは各点が回転する速度を表します。t=2Piで全ての点が元に戻るようにしたかったので、係数は整数である必要があります。この値は、実は手動で調整しました。手で描くとわかるのですが、「2」も「9」も一度普通になぞって逆方向もなぞり元の位置に戻ると、X軸方向は3往復、Y軸方向は1往復します。このことから縦に動く線とその3倍の速さで横に動く線ができることを期待してp5～p8でのtの係数、すなわちは惑星の公転速度は1および3と決定しました。これは上手く働いたようです。p1～p4でのtの係数、これは衛星の公転速度になりますが、いろいろと試行した末に上の形になっています。今見ると変な係数の付け方で計算していますが、どうせ最適化すれば同じ式になる（ことが期待できる）ので気にしなくて良いかと思います。
　さて、pにどのような軌跡を描かせたいのかを決めます。これはフリーハンドで描いてその座標を読みました。まず「2」はこんな感じです。

t=0Pi (-0.437739, 0.64179)
t=1Pi/28 (-0.286784, 0.88061)
t=2Pi/28 (-0.0754772, 0.940305)
t=3Pi/28 (0.166043, 1.)
t=4Pi/28 (0.437739, 0.940305)
t=5Pi/28 (0.588657, 0.910438)
t=6Pi/28 (0.679222, 0.761181)
t=7Pi/28 (0.739611, 0.611923)
t=8Pi/28 (0.709435, 0.432838)
t=9Pi/28 (0.588657, 0.223886)
t=10Pi/28 (0.528304, 0.164191)
t=11Pi/28 (0.407526, -0.0149336)
t=12Pi/28 (0.226395, -0.164191)
t=13Pi/28 (0.0150882, -0.283581)
t=14Pi/28 (-0.166043, -0.432838)
t=15Pi/28 (-0.316997, -0.582095)
t=16Pi/28 (-0.467915, -0.731353)
t=17Pi/28 (-0.649046, -0.850743)
t=18Pi/28 (-0.8, -0.970133)
t=19Pi/28 (-0.679258, -1.)
t=20Pi/28 (-0.467915, -1.)
t=21Pi/28 (-0.256608, -0.970133)
t=22Pi/28 (-0.0754772, -0.970133)
t=23Pi/28 (0.13583, -0.940305)
t=24Pi/28 (0.286784, -0.940305)
t=25Pi/28 (0.437739, -0.970133)
t=26Pi/28 (0.618869, -0.970133)
t=27Pi/28 (0.739611, -0.940305)
t=28Pi/28 (0.8, -0.910438)
t=29Pi/28 (0.739611, -0.940305)
t=30Pi/28 (0.618869, -0.970133)
t=31Pi/28 (0.437739, -0.970133)
t=32Pi/28 (0.286784, -0.940305)
t=33Pi/28 (0.13583, -0.940305)
t=34Pi/28 (-0.0754772, -0.970133)
t=35Pi/28 (-0.256608, -0.970133)
t=36Pi/28 (-0.467915, -1.)
t=37Pi/28 (-0.679258, -1.)
t=38Pi/28 (-0.8, -0.970133)
t=39Pi/28 (-0.649046, -0.850743)
t=40Pi/28 (-0.467915, -0.731353)
t=41Pi/28 (-0.316997, -0.582095)
t=42Pi/28 (-0.166043, -0.432838)
t=43Pi/28 (0.0150882, -0.283581)
t=44Pi/28 (0.226395, -0.164191)
t=45Pi/28 (0.407526, -0.0149336)
t=46Pi/28 (0.528304, 0.164191)
t=47Pi/28 (0.588657, 0.223886)
t=48Pi/28 (0.709435, 0.432838)
t=49Pi/28 (0.739611, 0.611923)
t=50Pi/28 (0.679222, 0.761181)
t=51Pi/28 (0.588657, 0.910438)
t=52Pi/28 (0.437739, 0.940305)
t=53Pi/28 (0.166043, 1.)
t=54Pi/28 (-0.0754772, 0.940305)
t=55Pi/28 (-0.286784, 0.88061)

また「9」はこうなりました。

t=0Pi (0.65, 0.830183)
t=Pi/21 (0.273674, 0.924517)
t=2Pi/21 (0.034217, 1.)
t=3Pi/21 (-0.273674, 1.)
t=4Pi/21 (-0.444739, 0.924517)
t=5Pi/21 (-0.65, 0.716971)
t=6Pi/21 (-0.615804, 0.547154)
t=7Pi/21 (-0.513173, 0.396214)
t=8Pi/21 (-0.205282, 0.377363)
t=9Pi/21 (0.034217, 0.433969)
t=10Pi/21 (0.239478, 0.547154)
t=11Pi/21 (0.444739, 0.667552)
t=12Pi/21 (0.513132, 0.735849)
t=13Pi/21 (0.615762, 0.84906)
t=14Pi/21 (0.513132, 0.547154)
t=15Pi/21 (0.478935, 0.45282)
t=16Pi/21 (0.478935, 0.30188)
t=17Pi/21 (0.410501, 0.0188773)
t=18Pi/21 (0.273674, -0.226423)
t=19Pi/21 (0.136848, -0.509426)
t=20Pi/21 (-0.0000205261, -0.811306)
t=21Pi/21 (-0.136848, -1.)
t=22Pi/21 (-0.136848, -1.)
t=23Pi/21 (-0.0000205261, -0.811306)
t=24Pi/21 (0.136848, -0.509426)
t=25Pi/21 (0.273674, -0.226423)
t=26Pi/21 (0.410501, 0.0188773)
t=27Pi/21 (0.478935, 0.30188)
t=28Pi/21 (0.478935, 0.45282)
t=29Pi/21 (0.513132, 0.547154)
t=30Pi/21 (0.615762, 0.84906)
t=31Pi/21 (0.513132, 0.735849)
t=32Pi/21 (0.444739, 0.667552)
t=33Pi/21 (0.239478, 0.547154)
t=34Pi/21 (0.034217, 0.433969)
t=35Pi/21 (-0.205282, 0.377363)
t=36Pi/21 (-0.513173, 0.396214)
t=37Pi/21 (-0.615804, 0.547154)
t=38Pi/21 (-0.65, 0.716971)
t=39Pi/21 (-0.444739, 0.924517)
t=40Pi/21 (-0.273674, 1.)
t=41Pi/21 (0.034217, 1.)

　改めてみると字が下手すぎですね。フリーハンドで描いたのは、別に味のある文字を使いたかったとかいうのではなく、単に時間が無かったからに他なりません。
　時刻tにおけるpの座標をp(t)とします。「2」の場合はt=0～55Pi/28、「9」の場合はt=0～41Pi/21まで変化させ、上で求めた各点とp(t)との距離の2乗をそれぞれの時刻で計算し、合計します。この式は24変数を含むスカラーになることがわかります。この巨大な式を評価式として最小化すれば良いわけです。ただ評価式が複雑すぎて解析的な解は求まらないかと思います。そこで最小化には遺伝的アルゴリズムと滑降シンプレックス法を組み合わせた自作のプログラムを使用しました。なお滑降シンプレックス法は、名前こそ似ていますが線形計画法で登場するシンプレックス法とは直接の関係は無いので注意が必要です。
　結果はこのようになりました。
「2」の場合

a1=-0.281248,a2=-0.217597,a3=-0.133595,a4=0.135594,
a5=0.233587,a6=-0.533821,a7=-0.865466,a8=-0.774918,c5=0.296162,
c6=0.128995,c7=-0.0170449,c8=-0.310357,d5=1.06265,d6=-0.197444,
d7=-0.199307,d8=-0.308207,e1=0.311494,e2=0.218138,e3=0.0167031,
e4=1.52846,e5=0.952566,e6=0.0411405,e7=1.25651,e8=1.3886

「9」の場合

a1=-0.197302,a2=0.145299,a3=-0.197725,a4=-0.200862,
a5=-0.524991,a6=0.779246,a7=-0.472475,a8=1.29185,c5=0.117048,
c6=0.560647,c7=0.611663,c8=-0.875045,d5=0.415859,d6=-0.461184,
d7=1.32743,d8=-0.340324,e1=1.35564,e2=-0.552802,e3=-1.25281,
e4=0.787604,e5=0.866171,e6=-0.150847,e7=-0.496076,e8=4.18745

　円の半径a*の値が負になったりしていますが、p1～p8の式のまま素直に描画すれば良いので特に問題はありません。最適化の限界で、またそもそも変数の数の限界から、元の点とぴったり同じ軌跡を描くことはできません。できあがった動画を見ても、往復で軌跡がずれているのがわかります。なお、円の数を増やせばもっと綺麗な文字の形になるのかとのコメントがありましたが、確かに綺麗になる可能性はありますが変数の数が増えるので探索空間が爆発的に増大し最適化が困難になるかもしれません。この辺りのさじ加減にはなかなか難しいものがあります。
　以上が「2」と「9」の描き方です。「0」については楕円で代用しています。あのように棒を動かすとその棒上の点が楕円を描くことは簡単に証明可能です。
　好評だったので来年の元旦にもまた何か投稿しようと思います。今度は遺伝的プログラミングなどを使って、どの円とどの円を組にするかなども含めて最適化すると面白いかな、と思っています。あるいは全く違う生物的なネタも考えています。ただ、「2010」という文字が簡単すぎて今ひとつ創作意欲をかき立てられないですね。まだ1年間あるので、いろいろと考えて見ることにします。

　ご無沙汰にもほどがありますが、とりあえず生存報告ということで。
いろんなことを放ったらかしにしたままですみません。
　最近機械式計算機に凝っていて、リヒテンシュタイン製のクルタ計算機をいくつか手に入れました。分解してみるとその巧みな機構に感嘆するばかりです。この計算機が生まれた過程もとても興味深い。
　そのうち、ニコニコにクルタ関係の動画を投稿するかもしれません。もしかすると冬コミで同人誌作ろうか、という話もあります。
　期待せずに待っていてください。

これの更に続きです。WD5の位置計算は、太陽とWD5の2体問題でやっていたのですが、今回は火星がWD5に及ぼす引力の影響も考慮した微分方程式で計算してみました。

20000モンテカルロですが、うーん、結果は変わりませんね。平均、標準偏差もほぼ同等です。火星に接近するケースが増えるかと思ったのですが。ちなみに正規分布を仮定して平均・標準偏差から衝突確率を求めると、約0.0031%となります（前回も）。

　前回の日記の続きです。
　すみません、計算の途中で足し算と引き算を間違えていました。
改めて、正しい方法で再計算しました。モンテカルロ部分のみの修正なので、中心ケースの最接近時刻・距離は変わりません。
　CGはこんな感じです。誤差を考慮した100シミュレーションの軌跡を細い線で表示しています。

　モンテカルロは、50000回に増やしてみました。結果はこうなります。

　1000kmごとのヒストグラムも作ってみました。

　正規分布っぽい結果で良い感じです。このうち火星に衝突するのは2ケースでした。2/50000=0.004%ということでNASAの発表0.01%と比べると少ないのですが、ケースが少ないので何とも言えません。マルコフ連鎖モンテカルロ法とかありますが、知りたいケース付近を集中的に試行するのでこういう場合に使われるのでしょうね。いずれにせよ、NASAの発表と矛盾は無いと思われます。ただ、2体問題の計算しかやっていないので、若干距離が遠くなる方向に計算が偏っているはずです。次に同じような天文ショーがあるときは、N体問題で解きたいと思います。
　というわけで、NASA発表の追試は終了です。NASAの発表文を読むと、なんとなくモンテカルロじゃなくて解析的に確率を求めているような気がします。どんな式を使っているのでしょうか。
結論

　WD5 2007 が火星にぶつかる可能性はほとんどありません。

　初っぱなからVOCALOIDと何の関係もない話題で申し訳ないです。小惑星衝突のお話です。
　先月ぐらいから小惑星が火星にぶつかる、ぶつからないと天文ファンの間で話題になっていて、一時は1/25の確率で衝突するという計算が発表されたりしていました。
　私も独自に計算してニコニコ動画に投稿するためにCGも作りかけていたのですが、計算をやり直したりCGを作り直したりとぐずぐずしているうちに他の人に先を越され、更にNASAにも「衝突の可能性はほとんど無くなった」と発表されてしまいました。
　というわけで完全に時機を逸した感があるのですが、計算結果をメモ代わりに書いておこうと思います。WD5の軌道要素はNASA/JPLのページから持ってきました。最新版は1/9のデータで、NASAの発表もこれに基づいています。
　まず最接近時のCGはこんな感じになります。WD5の大きさは誇張して描いています。 実際は直径50m程度なので点にしか見えないはずですね。

　やはりぶつかりそうにありません。横軸に日付、縦軸に距離(天文単位)のグラフを作るとこうなります。

　最接近は1/30 20:59:55（日本時間）で、0.000175006AU=26180.5kmまで近づくということになりました。これはNASAの発表とほぼ一致しています。
※以下の計算は間違っています。正しい計算はこちら。
　次に、衝突する確率を求めてみました。NASAのページに各軌道要素の誤差も載っているので、とりあえず正規分布で振ってモンテカルロを回してみました。独立した乱数ではなく、ちゃんと要素間の共分散行列も考慮した乱数を使っています。その結果がこれ。

　最接近時の距離(km)を10000回分プロットしたのですが、あれ？なんだか随分とばらついますね。火星の半径(3390km)以下の点が35もあります。つまり私の計算では、35/10000=0.35%の確率で衝突する可能性が残っていることになるのですが。
　中心値の計算はたぶん正しいので、誤差の解釈がおかしいのでしょう。ちょっと見直してみます。