３次元ベクトルを扱うための C++ class ファイルです。粒子計算等でよく用いる３次元ベクトル量を、int や double といった基本型と同じように扱います。
vector3.h ファイルを適切な場所に設置し、 C++ ソースコードの先頭でヘッダファイルを include します。 以下のどちらかの書き方で大丈夫でしょう。
#include "vector3.h" #include <vector3.h>
そして、下のように vector3 型の変数（オブジェクト）を宣言して使います。 （x,y,z）の３つの double 型のメンバ変数を備えています。 これらのメンバ変数は public 変数ですので、個別に値を参照することができます。 初期化方法はいろいろあります。 int 型の引数は double 型に変換されます。
vector3 a;
a.set( 2.0,0.0,0.0 ); // set メソッド
vector3 b( 1.0,0.0,0.0 ); // コンストラクタの引数
vector3 c( 1, 0, 0 ); // int 型を渡すことも可能
c.x = 0.0; c.y = 0.0; c.z = 1;
整数・実数を交えた演算やベクトル型どうしの演算は下のように簡単に記述できます。 配列を使って計算するより、プログラムの見通しは各段に良くなります。 他にもいくつかの演算やメソッドを備えています。
vector3 v; double d; v = 3 * a + b / 1.5 - c; // 四則演算 v += b; // 代入演算子（ベクトル足し算） v /= 1.5; // 代入演算子（実数との割り算） v = a * b; // a,b の外積 d = ( a % b ); // a,b の内積 if( a == b ) if( a != b ) // 論理演算 d = v.abs(); // 絶対値 d = v.abs2(); // 絶対値の２乗
内積の・(ドット)だけ不自然な表記になってしまいますが、 その他は自然な感じで記述できています。
注) 論理演算子は正しく動作する保証はありません。ベクトルの各成分の値は倍精度の浮動小数点ですから、数学的には同じ結果でも、計算機では値が一致しないことがあります。
以下のサンプルを vector3.h ファイルと同じディレクトリに適当な名前で保存してください。 この説明では、ファイル名を sample.cc とします。
#include <stdio.h>
#include <vector3.h>
int main( int argc, char* argv[] )
{
vector3 a,b,c,d;
a.set( 1.0, 0.0, 0.0 );
b.set( 0.0, 1.0, 0.0 );
c = a + b;
printf( "%lf %lf %lf\n", c.x, c.y, c.z );
d = a * b;
printf( "%lf %lf %lf\n", d.x, d.y, d.z );
return 0;
}
そして、次のようにしてコンパイルします。 -I. は、カレントディレクトリ内でヘッダファイル（vector3.h）を探すためのオプションです。 また、数学ライブラリを使っているので、オプション -lm も必要です。
$ g++ -I. sample.cc -lm
gcc の代わりに clang が入っている環境では次のようにすれば良いでしょう。
$ clang++ -I. sample.cc -lm
プログラムの実行結果は下のようになり、 ベクトリの和と外積を計算できていることがわかります。
$ ./a.out 1.000000 1.000000 0.000000 0.000000 0.000000 1.000000
ラグランジュ粒子の運動を扱うための C++ class ファイルです。 上記のベクトルクラスを使って構築されています。 質点にかかる外力を定義すると、 4段4次、あるいは7段6次の Runge=Kutta 法を使って 質点の運動を解き進めていくことができます。
vector3.h ファイルと同様に、ヘッダファイルを適切な場所に設置します。 particle クラスは、 ２つの外部ファイル（vector3 クラスと Runge=Kutta 法の定数ファイル）に依存しています。 vector3.h と RK.h の２つのファイルを パスの通ったディレクトリに置いてください。
実際の使い方は、サンプルファイル（test.cpp）を見た方がわかりやすいでしょう。 C++ ソースコードの先頭でヘッダファイルを include しています。
#include <particle.h>
そして、粒子に加わる外力を関数 F で記述します。 関数 F は、電荷、時間、位置座標（vector3 形式）を引数を取り、 それに応じた外力を vector3 形式で返します。 サンプルでは、ローレンツ力を３次元形式で記述しています。
vector3 F( const vector3& _r, const vector3& _v,
const double& _t, const double& _q ){
return _q * ( E(_r) + _v * B(_r) );
}
粒子は時間 t、質量 m、電気量 q といった内部変数を持っています。 これらの変数の値は、set...() と定義されたメソッド群で変更することができます。 同様に get...() メソッドを使って値を参照することもできます。 オブジェクト志向プログラミングでは クラスの継承時に内部変数を変えるバグができやすいので、 このようなアクセサメソッドを明示的に使って、 プログラムの安全性を高めることが好まれています。
particle p; p.sett(0); p.setm(1); p.setq(1); // p.setq(-1); p.setr(0.0,0.0,0.0); p.setv(0.3,0.0,0.1);
ただ、この粒子クラスはそれほど複雑なプログラムではありませんので、 ユーザーが十分注意すれなら、内部変数を直接参照しても良いと思います。 内部変数に直接参照する必要があれば、particle.h ファイルの
protected:
の部分で定義されている変数を
public:
の部分に移してください。
粒子のベクトル成分（座標・速度）は、 外部から直接参照できるようにしています。
コンパイル方法はベクトルクラスと同じです。 プログラムの振る舞いは、実行結果（粒子軌道と速度成分）を見れば一目瞭然でしょう。
$ g++ -I. test.cc -lm $ clang++ -I. test.cc -lm
プログラムの内部では、4段4次および 7段6次の Runge=Kutta 法を用いて 粒子軌道を計算しています。7段6次 Runge=Kutta 法の係数は、 係数が簡単な Butcher の３番目の公式 [1,2] を利用しています。
ローレンツは、簡単化した大気モデル方程式の振る舞いを調べ、 系の時間発展が初期値の違いに鋭敏に反応することを発見しました。 そして、この仕事を契機に、現代的なカオス研究がはじまったとされています。 粒子クラスを使って、ローレンツモデルの振る舞いを解いてみましょう。
ローレンツの方程式系 [1] は
\begin{align} \dot{x} &= \sigma (y-x) \\ \dot{y} &= x (r-x) - y \\ \dot{z} &= xy - b z\\ \end{align}で、３つのパラメーターは $\sigma=10.0$, $r=28.0$, $b=8/3$ です。 ちょっとトリッキーですが、この式の２階微分を考えます。
\begin{align} \ddot{x} &= - \sigma \dot{x} + \sigma \dot{y} \\ \ddot{y} &= - \dot{x}z - x\dot{z} + r \dot{x} - \dot{y} \\ \ddot{z} &= \dot{x}y + x\dot{y} - b \dot{z}\\ \end{align}これを外力ベクトル F として使うと、 粒子クラスでローレンツ方程式を解くことができます。 結果を以下に示します。２つの特異点のまわりで、 系が不思議な振る舞いをしているのがわかります。 これが、ストレンジアトラクターの中で 最も有名なローレンツアトラクターです。
不明な点がありましたら直接ソースファイルを見るか、あるいは作者にメールでご連絡下さい。バグの報告やその他の要望も歓迎します。