Oct 26, 2008

Blosxom環境構築(2/2)

今回このblogのBlosxom版を構築するために行ったことをまとめる。

●Blosxom本体
SourceForge.netのBlosxom projectから最新版を入手した。
なお、同ページは"The Unofficial Blosxom User Group"では"SF.net Project"と書かれてリンクされている。

●プラグイン
Blosxomは本体が小規模に作られているので、欲しい機能のプラグインを入れる所から始まる。
今回はblosxom starter kit(bsk)とBlosxom Plugin Registryを参考にプラグインを選び、"The Unofficial Blosxom User Group"のV2 Plugin Registoryから最新版を取得していった。
以下、インストールしたプラグイン名と、インストールに伴って改造した点を記す。
・archives
月毎のアーカイブページと、それらへの階層状のリンクを生成する。
特に改造無し。
・back_and_forth
個別エントリーのページから、更新日付が1つ前のページ、1つ後のページへのリンクを生成するプラグイン。
bskの"«前 | Main | 後»"という並びが気に入ったので、bskのを使わせてもらった。(V2 Plugin RegistoryのはMainの部分が出せない)
・categories
カテゴリー毎のアーカイブページと、それらへの階層状のリンクを生成するプラグイン。
ディレクトリ名と表示文字列の対応($friendly_name)を埋めたのと、'all entries'を'全てのエントリー'に書き換えたのみ。
※report_dir_start呼び出しがコメントアウトされてるバージョンだと、'all entries'を表示するにはコメントアウトを解除する必要あり。
・entries_index
エントリーの作成日時を保存する。これが無いとエントリーの作成日時がエントリーファイルの更新日時になるので、FTPでwebサーバーにアップロードするとファイルの日時がアップロード時刻になり、エントリーの作成日時が失われてしまう。
特に改造無し。(インデックスファイルは別途作成して一緒にアップロードした)
・entry_title
個別エントリーページのタイトル(<TITLE>タグの中身)をそのエントリーのタイトルにする。
"ブログ名 セパレータ タイトル"という並びが気に入ったので、Kyo Nagashimaさんのバージョンを使用。特に改造無し。
・find
ブログ内をテキスト検索して、マッチするエントリーだけを出力するプラグイン。
特に改造無し。
・geek
(出力例)。いらないけど消すには惜しいので残しておく。
・seemore
複数エントリーが表示されるページで、各エントリーを丸ごと表示するのでなく、部分的に表示して、残りの部分へのリンクを付けるようにする。
"See more ...."を「続きを読む」に変えかけたが、"See more ...."が格好よくて気に入ったので、改造無し。
・titles_index
複数エントリーが表示されるページで、各エントリーをタイトルのみで表示するようにする。
seemoreで表示するトップページ以外では表示エントリー数の制限を緩和するため、blosxom本体のgenerate{}の"# Stories"の所の

my $ne = $num_entries;
my $ne = $currentdir ? 100 : $num_entries;
に変更。
・writeback
後述。

●Writebackスパム対策
Blosxom公式サイトのWritebackプラグインを使うと、スパムの餌食になる。対策として、以下の5つのプラグインを見つけた。
(1) WritebackPlus
(2) Feedback
(3) spam_blocker
(4) hail2u.netのWriteback
(5) Writeback with Akismet
まずデファクトスタンダードと思われる(1)のWritebackPlusを試そうとしたが、その中で使われているMT-blacklistが既に配布されておらず(MovableType 3.2で他の方法が使われるようになったため)、スパム対策にならなかった。
ネット上で流し読みする限り(2)はマイナーで、(3)と(5)は導入が面倒くさそうで、(4)の方法で十分効果がありそうなので、今回はWritebackPlusに(4)の方法を適用することにした。
WritebackPlusの改造結果
(WritebackPlusはVersion 0.2.1aを使用)
改造点(コメントで"Y.Nomura"とある部分)
・250行目付近の$pathを修正
・257行目付近、(4)のWritebackの100~131行目をコピー
・305行目付近、blacklist無しでも動くようにコメントアウト
・375行目付近、Net::SMTP無しでも動くようにコメントアウト(Hi-Hoのサーバーに無かったため)

●MovableTypeのエントリーのインポート
MovableTypeの管理ページの「エントリーの書き出し」を使って全エントリーを1つのテキストファイルに出力し、Perlスクリプトをガリガリ書いてBlosxom用のファイルに変換した。

●Hi-Ho依存の設定
参考までに、Hi-HoのレンタルWebサーバーでBlosxomを動かすために設定した、blosxom.cgi内の設定項目の一部を挙げておく。(一部伏字)
$dataurlは、テンプレートHTML(フレーバー)内からの他ファイル(スタイルシートや画像)参照で必要になる絶対パス。フレーバー内の相対パスは特殊なCGIディレクトリからの相対パスになり、データディレクトリに辿り着けない。

# Where are this blog's entries kept?
$basedir = "$ENV{'HOME'}/html/*************";
$datadir = "$basedir/entries";

# What's my preferred base URL for this blog (leave blank for automatic)?
$url = "http://www.sam.hi-ho.ne.jp/cgi-bin/user/ynomura/blosxom.cgi";
$dataurl = "http://www.sam.hi-ho.ne.jp/~ynomura/blosxom";

# Where are my plugins kept?
$plugin_dir = "$basedir/plugins";


See more ...

WriteBacks

Oct 25, 2008

Blosxom環境構築(1/2)

先月、このweblogが動く自宅サーバーへのアクセス手段でお世話になってる、無料のDDNSサービス、ieServer.Netが不調で、しばしば外からここにアクセスできなくなっていた。
当初から大変ありがたく使わせてもらっているが、無料のサービスなので十分な設備で運営されるとは限らないし、いつサービス終了となるかもわからない。大してアクセス数がある訳でも無いが、せっかくネチネチと綴ってきた現役のweblogなので、まだ消滅させるには惜しい。そこで、予備のweblog環境を検討することにした。

このサーバーのblogは、MovableTypeを使って記事を静的なHTMLに変換して作っており、一度HTMLに変換すれば、WebサーバーにMovableTypeが無くても閲覧可能である(但しコメント/トラックバックの受信やコンテンツ内検索等はできない)。従って、MovableTypeを使って作成した静的なHTMLをプロバイダーが提供するWebサーバーにアップロードするだけで、公開はできる。
…と考えて作業を始めると、1つのMovableType環境で複数のURL向けに静的なHTMLを出力するのは容易でないことがわかった。MovableTypeの設定ファイルを2つ用意して手で切り替えるか、データベース内に全データを丸ごと2つ持たないといけないのだ。あまりにもすっきりしないし、不便である。

そこで、プロバイダーのWebサーバー側にブログシステムを構築することを考えた。
世の中にはMovableType以外にNucleus、tDiary、華式、blosxom等様々なブログツールや日記ツールがあるが、私が使用しているプロバイダーであるHi-HoのレンタルWebサーバーは持ち込みのプログラムを動かす制約が厳しく、ほとんど素のPerlのCGIしか動かせない。従って、多数のPerlモジュールを必要とするMovableTypeも、PHPで動く華式も、SQLデータベースを必要とするNucleusも動かない。
そんな中、blosxomは動いた。

実はこのblogを始める以前に、Kyo Nagashimaさんのblosxom starter kit(以下bsk)を使ってHi-Ho上でblosxomを動かして、少しだけblogしていた。そのbskをインストールした4年前、Hi-Hoのサーバー上で動かせたフリーウェアのブログ/日記ツールはblosxomだけだったが、現在やはり動かせるのはblosxomだけのようである。

Blosxomの公式サイトのページは更新が止まっており、そのサイトにあるblosxom本体も5年前からバージョンアップされていない。ユーザーも減っているようで、ユーザー自身が「なぜ今更blosxomを使っているかというと…」と自嘲気味に言い訳しているのを見かける。ユーザーがいなければノウハウを公開する人もいない。もうblosxomを使うのは厳しいか…と思いかけた所で、"The Unofficial Blosxom User Group"というサイトを見つけた。どうやら今年に入っても有志によってバージョンアップされ、プラグインも追加されていってるようだ。
そして、日本でもInfoseekの無料Webレンタルサーバーでblosxomしている人がいた。InfoseekでもHi-Hoでも動くのなら、レンタルウェブサーバー間のポータビリティは十分期待できる。一旦システムを構築したら、無料のレンタルサービスででもその環境を入れれば全く同じブログを使い続けられるのだ。エコロジーだ。

ひと安心してblosxomをバージョンアップしようとすると、以前にインストールしたblosxomが動かなくなっていた。Perlスクリプトの途中でエラーになってたようで、サーバーのPerlの仕様が変わったのか?とげんなりしたが、仕方なく解析すると、spamを受けまくってwriteback(コメント、トラックバック)ファイルが巨大になってたのが原因だと判明した。インストールしていたのはbsk 1.02そのままの環境で、当時spamは全く無かったので、spam対策はしていなかった。しかも、パスワードを設定していたはずなのだが、有効にしていたwikiedish(ブラウザ上でエントリーを編集するためのプラグイン)を使って、全てのエントリーがspamによって書き換えられていた。
辛うじてbsk 1.02を初期状態に戻して動かすことはできたが、blosxom本体だけをバージョンアップする方法と、spam対策されたWritebackPlusを組み込む方法がわからなかったので、bskを参考にして1からblosxom環境を構築することにした。

See more ...

WriteBacks

Maximaでマンデルブロー図形

お盆なので、特にお盆とは関係ないが、ついでに定番のマンデルブロー集合も描いておく。

f(c,z):=c-z^2$
density(x,y):= block([c,z,d,i],
c:x+y*%i,
z:c,
d:0,
for i:0 while (abs(z)<2 and d<40) do (z:f(c,z),d:d+1),
d
)$
plot3d(density,[x,-0.5,1.5],[y,-1.5,1.5],[grid,100,100],
[gnuplot_pm3d,true],[gnuplot_preamble,"set view map;unset surface;set xrange [-0.5:1.5];set yrange [-1.5:1.5];"]);

See more ...

Posted at 22:50 in PC一般 | WriteBacks (0)
WriteBacks

Maximaで微分方程式

せっかくMaxima使い始めたので、慣れるためにさらにMaximaのマニュアルを読み進めて遊んでると、Maximaは微分方程式も解けることがわかった。難しそうな響きを擁して実際難しい、物理学の基本ツールである。一瞬気が進まなかったが、大したものだと思ったので、敬意を表して試してみた。
すると、微分方程式が見えた気分になるようなグラフが出せることがわかり、面白かったので、ここに書き記すことにした。

まず、コマンドラインで微分させてみる。diff()が微分演算だ。

例によってTeXmacsを使って出力を整形している。(%i??)に続く青い部分が入力した部分、(%o??)に続く黒い部分がMaximaの出力だ。diff()の2つ目の引数は何で微分するか、3つ目の引数(省略可)は何回微分するかを意味する。

次に、desolve()を使ってシンプルな微分方程式\frac{d^2}{dx^2}y+\frac{d}{dx}y+y=0を解かせてみる。

desolve()の2つ目の引数は、どの関数について解くかを示す。diffの前の'(アポストロフィー)は、そのLisp関数を評価(実行&展開)せず、そのまま渡すことを指示するものである(desolve()の場合はこれを省略して評価してしまっても良いようだが、後述のode2()では不可)。y(x)の(x)はyがxの関数であることを示し、これを省略するとdy/dxが0になってしまう。
atvalue()を使って、y,y'の初期値を与えてから再度解かせてみる。

入力行の行末に;の代わりに$を付けると、計算結果を出力しない。
念のため、元の方程式に当てはまるかどうか、検算してみる。

%は直前の計算結果を表す。%o7は見るからに0なのだが、なぜか計算してくれなかったので、ratsimp()で整理させた。

常微分方程式(未知関数が1つだけの微分方程式)はode2()を使って解く方法もある。ode2()を使う場合、初期値は一般解の計算結果に対してic1(),ic2(),bc2()などを使って与える。同じ方程式を解かせてみる。

ode2()の2つ目の引数は常微分方程式の未知関数、3つ目の引数はその関数が従属する変数だ。上の場合は、引数からyがxの関数であることがわかるので、yをy(x)と書かなくても良い。

desolve()とode2()の能力を比較するため、解くには少し難しそうな微分方程式xy+\frac{d}{dx}y=0を解かせてみる。

ode2()は難なく解いた。念のため検算しておく。

それに対して、同じものをdesolve()に解かせると、下のような奇妙な式になった。

おそらく、ラプラス変換してx,yについて解いてラプラス逆変換したもの、という意味なのだろうが、これでは使えないだろう。これでも使える解なのかどうかを一応確認してみる。

やはり解けていない、使えない解であることがわかる。常微分方程式を解くのはdesolve()よりode2()の方が優れているようだ。

しかし、desolve()は複数の未知関数を含む連立微分方程式も解ける。
\left\{ \begin{array}{lll} \frac{d}{dx}f(x) & = & \frac{d}{dx}g(x)+\sin{x} \\ \frac{d^2}{dx^2}g(x) & = & \frac{d}{dx}f(x)-\cos{x} \end{array} \right.

e1とe2が方程式、f(x)とg(x)が未知関数だ。desolve()には、上のように、方程式と未知関数をそれぞれリスト形式で与える。

次に、plotdf()を使って微分方程式のグラフを出してみる。plotdf()に与える微分方程式は、dy/dx=F(x,y)または[dx/dt, dy/dt]=[F(x,y), G(x,y)]に限定されている。

load("plotdf");
plotdf(x);
これを実行すると、下のような画面が開く。(FreeBSD 6.xではTcl/Tkをインストールしていても「wishが無い」というエラーが出たが、wishという名のwish?.?へのリンクを作ればOKだった)

たくさんある矢印は、その地点(x,y)におけるdy/dxの向きと大きさを表している。このグラフ上のどこかをクリックすると、その座標を初期値とするグラフが描かれる。下の図は(0,0)辺りをクリックしたもの。

さらに色々な所をクリックする。

次の2つの図は、

plotdf(sin(x));
を実行したものである。

xとyを含む式も試してみる。

plotdf(sin(x)+cos(y));

放物線を横にした、x=y^2を出そうとした。両辺をxで微分すると1=2y(dy/dx)なので、dy/dx=1/(2y)である。よりシンプルに、dy/dx=1/yで試そうとして、

plotdf(1/y);
とすると、原因がよくわからないが、Tclのエラーになった。しかし、次のように小さなオフセットを与えると、何か描画された。2つ目以降の引数に色々なオプションを指定しているが、これらが無くてもエラーにはならない。
plotdf(10^-10+1/y,[trajectory_at,0.5,1],[nsteps,200],[xradius,3],[yradius,6],[xcenter,3],[ycenter,0]);

trajectory_atで初期値を(0.5,1)として描画したものだが、左側が何かおかしい。さらに何ヶ所かクリックしてみると、

y=0の近辺でおかしくなるようだ。1/(y+⊿y)-1/yの計算でもしているのだろうか。

x=cos(t), y=sin(t)のグラフは円になるので、[dx/dt, dy/dt]=[-y,x] (=[-sin(t), cos(t)])のグラフも円になる(図は省略)。それよりもっと場を複雑にさせようと考えて、dx/dtとdy/dtを周期関数にしてみた。次のグラフは[dx/dt, dy/dt]=[sin(y),sin(x)]である。

plotdf([sin(y),sin(x)]);

色を変えながら初期値をクリックした。セル状のグラフが無限に続く。
dx/dtとdy/dtを入れ替えると、次のようなグラフになる。
plotdf([sin(x),sin(y)]);

こんなのも作ってみた。

plotdf([
10^-10 - sin(atan(y/x))*cos(%pi/4*sin(8*atan(y/x))) - cos(atan(y/x))*sin(%pi/4*sin(8*atan(y/x))),
10^-10 - sin(atan(y/x))*sin(%pi/4*sin(8*atan(y/x))) + cos(atan(y/x))*cos(%pi/4*sin(8*atan(y/x)))
]);

[dx/dt, dy/dt]=[F(x,y), G(x,y)]で表せない(dx/dtにtが入るとか)方程式も、desolve()で解けるものなら
、いつものplot2d()でグラフにできる。
\left\{\begin{array}{l}\frac{dx}{dt}=3x-7y\\ \frac{dy}{dt}=5x-4y \end{array}\right.
(この方程式は[dx/dt, dy/dt]=[F(x,y), G(x,y)]で表せるが、plotdf()ではおかしなグラフになる)

atvalue(x(t),t=0,1);
atvalue(y(t),t=0,1);
Z:desolve([diff(x(t),t)=3*x(t)-7*y(t), diff(y(t),t)=5*x(t)-4*y(t)],[x(t),y(t)]);
plot2d([parametric, x(t), y(t), [t,0,6]],[nticks,200]),Z;

アステロイド曲線\left\{\begin{array}{l} x''(t)=2y'(t)-3x(t)\\ y''(t)=-2x'(t)-3y(t) \end{array}\right.(実はx=cos(t)^3, y=sin(t)^3から逆算したもの)を描いてみる。

atvalue(x(t),t=0,1);
atvalue(y(t),t=0,0);
atvalue(diff(x(t),t),t=0,0);
atvalue(diff(y(t),t),t=0,0);
Z: desolve([
diff(x(t),t,2)=2*diff(y(t),t)-3*x(t),
diff(y(t),t,2)=-2*diff(x(t),t)-3*y(t)
],[x(t),y(t)]);
plot2d([parametric,x(t),y(t), [t,0,2*%pi]],[nticks,200]),Z;

[parametric, ...]の部分を、それらのリストにすることにより、複数のグラフを1枚に描ける。

eq: [
diff(x(t),t,2)=2*diff(y(t),t)-3*x(t),
diff(y(t),t,2)=-2*diff(x(t),t)-3*y(t)
];
atvalue(x(t), t=0, 1);
atvalue(y(t), t=0, 0);
atvalue(diff(x(t), t), t=0, 0);
atvalue(diff(y(t), t), t=0, 0);
Z1: desolve(eq, [x(t), y(t)]);
atvalue(x(t), t=0,2);
Z2: desolve(eq, [x(t), y(t)]);
atvalue(x(t), t=0,3);
Z3: desolve(eq, [x(t), y(t)]);
L1: [parametric, x(t), y(t), [t, 0, 2*%pi]], Z1;
L2: [parametric, x(t), y(t), [t, 0, 2*%pi]], Z2;
L3: [parametric, x(t), y(t), [t, 0, 2*%pi]], Z3;
plot2d([L1,L2,L3], [nticks,200]);

See more ...

Posted at 22:50 in PC一般 | WriteBacks (0)
WriteBacks

Maximaでドラゴンカーブ

MathematicaとかMapleとかの数式処理ソフトは、何であんなに高いのだろう。2008/7/26現在、個人向けライセンス価格はMathematicaは40万円超え、Mapleも30万円超えだ。教育機関にも政府機関にも属していない、日曜数学ファンの私にはとても手が出せない。
会社でも研究職とかでなく、ソフトウェア開発の付帯業務をやってるので、数学は一切使っておらず、買ってもらえる見込みは無い。バグ収束曲線が云々で必要とか言っても、そんな統計処理ならExcelで十分と言われてしまうであろう。
従って、MathematicaやMapleに触れる機会が無い。

そんな折、Maximaというフリーの強力な数式処理ソフトがあることを知ったので、今回ヤボ用があって、インストールして使ってみた。
とりあえず、テーラー展開させてみる。

 KURO-BOX> maxima

Maxima 5.10.0 http://maxima.sourceforge.net
Using Lisp GNU Common Lisp (GCL) GCL 2.6.7 (aka GCL)
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
This is a development version of Maxima. The function bug_report()
provides bug reporting information.
(%i1) taylor(sin(x),x,0,10);
                          3    5      7       9
                         x    x      x       x
(%o1)/T/             x - -- + --- - ---- + ------ + . . .
                         6    120   5040   362880
(%i2) taylor(sqrt(r^2-y^2),y,0,10);
                    2      4      6         8       10
                   y      y      y       5 y     7 y
(%o2)/T/       r - --- - ---- - ----- - ------ - ------ + . . .
                   2 r      3       5        7        9
                         8 r    16 r    128 r    256 r
(%i3)
まあ言いたいことは十分わかるが、表示がなんか残念だ。
そこで、満足のいく表示を得るためにTeXmacsを使ってみる。以下の画面はMaxima 5.11.0とTeXmacs 1.0.6.11を使って実際に表示されたものである。%i??という青い部分が入力した部分だ。(表示環境はFreeBSD 6.3+Xorg 1.4.0)

積分もしてみる。

e^(-x^2)を-∞から+∞まで積分すると√πになるのは、手計算でやると面積分にして極座標変換してと技を尽くして偶然解ける感じの有名な例題だが、Maximaは区分求積による近似値計算ではなく、あっさり解析的に計算したようだ。

せっかくTeXmacsを使ってるので、見た目が派手な連分数を表示してみる。

cf(continued fraction)の出力は連分数のリスト表記になっていて、それをcfdisrepすると連分数表記になる。それをratで実数にすると(無限の連分数でないので元々実数だが)約分され、floatで小数にすると元の√2に限りなく近いことがわかる。
蛇足だが、cfdisrepしないと実数や小数に変換できない。

√3についてもやってみる。

方程式を解いてみる。solveで方程式を解析的に解ける。

解けない場合は、xの解にxが混ざったりする。(例は省略)

allrootsを使うと、多項式の方程式を数値的に計算できる。

指数の方程式を渡すと文句を言われる。

find_rootを使うと、(おそらくニュートン法による)近似計算で、解が1つだけ含まれる範囲の解を計算できる。

x^3=3^xのe~100の間の解は当然3なのだが、誤差が出た。


さて、せっかくフリーソフトなので、参考資料も無料のものを、と思って、図書館でMaximaの本を探したら、無かった。その代わり、あるMathematicaの本で、わずか十数行のコードでドラゴンカーブ(ドラゴン曲線)を描くMathematicaのコードを見つけた。マンデルブロー集合とかペアノ曲線とかのフラクタル図形の簡単なやつだ。Maximaの練習とMathematicaとの比較を兼ねて、その本を参考にMaximaで似たようなコードを書いてみた。

dragon(segments_):= block([p, q, r, s],
 p: map(
  lambda([x], [
   [x[1][1], x[1][1] + (x[1][2] - x[1][1]) * (1-%i)/2, x[1][2]],
   [x[2][1], x[2][1] + (x[2][2] - x[2][1]) * (1+%i)/2, x[2][2]]
  ]), segments_),
 q: map(
  lambda([x], [
   {[x[1][1], x[1][2]], [x[1][2], x[1][3]]},
   {[x[2][1], x[2][2]], [x[2][2], x[2][3]]}
  ]), p),
 r: map(lambda([x], listify(x)), flatten(q)),
 s: map(lambda([x], [realpart(x), imagpart(x)]), flatten(r)),
 plot2d([discrete, float(s)],
  [gnuplot_preamble, "set xrange [-1.2:1.2];set yrange [-1.2:1.2]"]),
 r
)$

p: [[[0, 1/2-%i/2], [1/2-%i/2, 1]]]$
for i:1 thru 10 do p:dragon(p)$
gnuplotが使える環境でこれを実行すると、以下のようなドラゴンカーブが描かれる。TeXmacsは必要ない。順に、L字型に配置した2本の線から始めて、全体を90°回転させて先端に繋げる、というのを1回、2回、…10回繰り返したものである。
アルゴリズムとしては、全ての線をL字に折り曲げる、というのを繰り返せばドラゴンカーブと同じ形になるというのを使っている。必要最小限にコードを簡略化して難解になったし、同様のアルゴリズムの説明は多くあると思うので、ここではコードの説明は省く。

See more ...

Posted at 22:50 in PC一般 | WriteBacks (0)
WriteBacks

円盤の衝突計算

弾性衝突円盤アプレットの衝突計算について、悲しいことに式を求めるのにかなり時間を費やしたので、忘れないように記録する。これが無いと、自分でもソースコードを理解できなくなってしまう。

まず直線上の衝突について考える。物体1,2の質量をm1,m2、衝突前の速度をv1,v2、衝突後の速度をv1',v2'とすると、完全弾性衝突の場合、運動量保存の法則より
m_1v_1+m_2v_2=m_1v'_1+m_2v'_2
-\frac{v'_1-v'_2}{v_1-v_2}=1
が成り立つ。速度変化に注目してこれを整理すると、
v'_1-v_1=\frac{-2m_2(v_1-v_2)}{m_1+m_2}
v'_2-v_2=\frac{-2m_1(v_2-v_1)}{m_1+m_2}
が得られる。これは、
v_{tmp}=\frac{m_1v_1+m_2v_2}{m_1+m_2}---(1)
と置くと、
v'_1-v_1=2(v_{tmp}-v_1)---(2)
v'_2-v_2=2(v_{tmp}-v_2)---(3)
と整理できる。

次に平面上の衝突について考える。
ball_hit.png
図のように円盤と円盤が衝突する場合の速度変化については、円盤の中心同士を結ぶ線上の方向の
速度が衝突によって変化し、衝突方向でない速度成分は変化しないので、衝突方向の速度V1r、V2rについて、上記の式(1)~(3)により速度の変化を求めれば良い。つまり、 Pを円盤の中心間を結ぶベクトル、V1,V2を円盤1,2の速度ベクトル、v1r,v2rを円盤1,2のP方向の速度、V1,v1rの衝突後の速度をそれぞれV1',v1r'とすると、円盤の速度変化は
\vec{V'_1}-\vec{V_1} = \frac{\vec{P}}{|\vec{P}|}(v'_{1r}-v_{1r})
である。

v1rの計算については、θをPとV1の間の角度とすると、
v_{1r} = |\vec{V_1}|\cos{\theta}
であり、PとV1との内積は
\vec{V_1}\cdot\vec{P} = |\vec{V_1}||\vec{P}|\cos{\theta}
だから、
v_{1r} = \frac{\vec{V_1}\cdot\vec{P}}{|\vec{P}|}
で求められる。同様に
v_{2r} = \frac{\vec{V_2}\cdot\vec{P}}{|\vec{P}|}
となる。

直線の場合の式(1)と同様に
v_{tmp}=\frac{m_1v_{1r}+m_2v_{2r}}{m_1+m_2} \left (=\frac{1}{|\vec{P}|}\frac{m_1<br>(\vec{V_1}\cdot\vec{P})+m_2(\vec{V_2}\cdot\vec{P})}{m_1+m_2}\right )
と置くと、
v'_{1r}-v_{1r}=2\frac{\vec{P}}{|\vec{P}|}(v_{tmp}-v_{1r})
となる。vtmpとv1rの両方に1/|P|が入っているので、計算上の便宜のため、
v_{1r} = \vec{V_1}\cdot\vec{P}
v_{2r} = \vec{V_2}\cdot\vec{P}
と置き直すと、
v'_{1r}-v_{1r}=2\frac{\vec{P}}{|\vec{P}|^2}(v_{tmp}-v_{1r})
v'_{2r}-v_{2r}=2\frac{\vec{P}}{|\vec{P}|^2}(v_{tmp}-v_{2r})
が得られる。2次元のベクトルをX軸成分とY軸成分に分けて
\vec{P}=(P_x,P_y), \vec{V_n}=(V_{nx},V_{ny})
とスカラー表記にすると、
|\vec{P}|^2 = P_x^2+P_y^2, \vec{V_n}\cdot\vec{P}=V_{nx}P_x+V_{ny}P_y
なので、円盤同士の衝突による速度変化は、
\vec{V'_1}-\vec{V_1} = \left (\frac{2P_x}{P_x^2+P_y^2}(v_{tmp}-v_{1r}),\frac{2P_y}{P_x^2+P_y^2}(v_{tmp}-v_{1r}) \right )
\vec{V'_2}-\vec{V_2} = \left (\frac{2P_x}{P_x^2+P_y^2}(v_{tmp}-v_{2r}),\frac{2P_y}{P_x^2+P_y^2}(v_{tmp}-v_{2r}) \right )
v_{tmp}=\frac{m_1v_{1r}+m_2v_{2r}}{m_1+m_2}
v_{1r} = V_{1x}P_x+V_{1y}P_y
v_{2r} = V_{2x}P_x+V_{2y}P_y
で求められることがわかる。

See more ...

Posted at 22:49 in 雑記 | WriteBacks (0)
WriteBacks