物理演算におけるベクトル型クラスのススメ
物理演算をコンピュータにやらせようとなると、ベクトルの計算は頻出である。
たとえば天文学とか惑星科学の分野とかでは、よく出てくるのが重力である。
\frac{\mathrm{d}^2 \vec{r}_i}{\mathrm{d} t^2} = - \sum_{j = 0}^{N_\mathrm{body}} \frac{m_j}{r_{ij}^2} \frac{\vec{r}_{ij}}{r_{ij}}
一言に重力と言っても、この業界では実は非常に奥が深い問題なのだが、まぁとりあえず今は差し当たって関係ないので省略しておきます。
重力は2体しかないなら解析的に解けるのだが、N体(N >= 3)になると解析的に解けないので、コンピュータ・シミュレーションに頼る事になる。
さてまぁとにかくなんでもいいんだけれど僕らはベクトル計算をする必要があるわけである。
どうするか？
こんなコードが思いつくだろうか。
for(int i = 0 ; i < Nbody ; ++ i){
for(int j = 0 ; j < Nbody ; ++ j){
double r_ij = 0;
for(int d = 0 ; d < 3 ; ++ d){
r_ij += (r[i][d] - r[j][d]) * (r[i][d] - r[j][d]);
}
r_ij = sqrt(r_ij);
for(int d = 0 ; d < 3 ; ++ d){
a[i][d] += - m[j] * (r[i][d] - r[j][d]) / (r_ij * r_ij * r_ij);
}
}
}
ちょっと待った！こんな書き方もできます！
for(int i = 0 ; i < Nbody ; ++ i){
for(int j = 0 ; j < Nbody ; ++ j){
vec3<double> dr_ij = r[i] - r[j];
double inv_r_ij = 1.0 / abs(dr_ij);
a[i] += - m[j] * dr_ij * (inv_r_ij * inv_r_ij * inv_r_ij);
}
}
なんと！たったの三行！
そう我々はバグを出さないためにこういう書き方をすべきなのだ。
各々次元にはどうせ同じ処理をしてしまうのだ、圧縮してしまえば良い。
そこをfor文で書くと、バグを生む可能性が出るし、for文の分計算コストがかかるかもしれない。
(もっとも、これはコンパイラの最適化で何とかなるかもしれないわけだが…)
さて、では何故こんな事ができるか？
もちろんC++にはvec3とかいう怪しい変数型は存在しない。
僕が作ったのである。
オブジェクト指向(と言っていいのかわからないけど)のC++では、
既存の型を組み合わせる事で、新しい型を作ることができる。
空間ベクトルは基本的に(x, y, z)の三成分なので、3つあればいいであろう。
基本的なレシピは以下の通りである。
template <typename type> struct vec3{
type x , y , z;
vec3(void){
x = y = z = 0;
}
vec3(type init_x , type init_y , type init_z){
x = init_x;
y = init_y;
z = init_z;
}
//operator =
vec3 operator=(const vec3& a){
this->x = a.x;
this->y = a.y;
this->z = a.z;
return *this;
}
vec3 operator=(const type a){
this->x = a;
this->y = a;
this->z = a;
return *this;
}
//operator +vec3
vec3 operator+() {
return *this;
}
//operator -vec3
vec3 operator-() {
return vec3(-x, -y, -z);
}
//operator +=
vec3 operator+=(const vec3& a){
this->x += a.x;
this->y += a.y;
this->z += a.z;
return *this;
}
//operator -=
vec3 operator-=(const vec3& a){
this->x -= a.x;
this->y -= a.y;
this->z -= a.z;
return *this;
}
//operator *=
vec3 operator*=(const type& a){
this->x *= a;
this->y *= a;
this->z *= a;
return *this;
}
//operator /=
vec3 operator/=(const type& a){
this->x /= a;
this->y /= a;
this->z /= a;
return *this;
}
//operator +
inline vec3 operator+(const vec3& a)const{
return vec3(this->x + a.x , this->y + a.y , this->z + a.z);
}
//operator -
inline vec3 operator-(const vec3& a)const{
return vec3(this->x - a.x , this->y - a.y , this->z - a.z);
}
//operator *
inline vec3 operator*(const type& a)const{
return vec3(this->x * a , this->y * a , this->z * a);
}
//operator /
inline vec3 operator/(const type& a)const{
return vec3(this->x / a , this->y / a , this->z / a);
}
~vec3(void){
}
};
template <typename type> inline vec3<type> operator*(type a , const vec3<type>& y){
return vec3<type> ( a * y.x , a * y.y , a * y.z );
}
template <typename type> inline type abs(const vec3<type>& a){
return sqrt(a.x * a.x + a.y * a.y + a.z * a.z);
}
template <typename type> inline type abs2(const vec3<type>& a){
return (a.x * a.x + a.y * a.y + a.z * a.z);
}
template <typename type> inline type inner(const vec3<type> &a , const vec3<type> &b){
return (a.x * b.x + a.y * b.y + a.z * b.z);
}
template <typename type> inline vec3<type> outer(const vec3<type> &a , const vec3<type> &b){
return vec3<type>(a.y * b.z - a.z * b.y , a.z * b.x - a.x * b.z , a.x * b.y - a.y * b.x);
}
さて、何か細かい解説でも入れようかと思っんたんだが、
テンプレートとかコンストラクタとかオーバーロードとか色々あって結構面倒なので、
個別の説明はやめて使い方だけ説明する事に。
基本的には僕らが知っているベクトルとして振る舞います。
以下のコードは3次元のベクトルを生成します。
vec3<double> v1;
各要素へのアクセスは"."を使います。
v1.x = 1.0;
v1.y = 2.0;
v1.z = -1.0;
生成時に初期値が決まっている場合、それを直接渡すことが出来ます。
vec3<double> v2 = vec3<double>(1.0, 3.0, 2.0);
加算や減算が使えます。
vec3<double> v3 = v2 + v1;
スカラー倍や除算が出来ます。
vec3<double> v4 = 2.0 * v3 / 3.0;
absで絶対値が帰ってきます。
double r4 = abs(v4);
innerとouterで内積・外積が帰ってきます。
double inner = abs(v4, v1);
vec3<double> outer = abs(v4, v1);
さて、まぁこんなところかしら。
自分用に作ったので結構癖があると思います。
実用の上では皆様自分で是非作ってみてくださいな。
なお、人に聞いたところでは、
内積は*、外積は^にするのが一般的だそうです。
へぇ。