C で数値計算しようと言う学生がいて、 Runge-Kutta をすると言うのだけど、 C にはベクトルがないので、なんだか悲しいプログラムになったりする。 以前それで Ruby やったりしたのだけど (結構面白かった)、 探せば行列・ベクトル用の C++ クラス・ライブラリィがあるかなあ?と思って、 Eigen というのを見つけた。
3.2.0.tar.gz を取得した。実際の名前は eigen-eigen-ffa86ffb5570.tar.gz だった。
こんな感じでインストール |
tar xzf eigen-eigen-ffa86ffb5570.tar.gz cd eigen-eigen-ffa86ffb5570 tar cf - Eigen | (cd /usr/include;sudo tar xpf -) |
最初に書いたプログラム |
#include <iostream> #include <math.h> #include <Eigen/Dense> using namespace Eigen; double m, g, Gamma, e; VectorXd f(double t, VectorXd x) { VectorXd y(4); y(0) = x(2); y(1) = x(3); y(2) = - Gamma / m * x(2); y(3) = - g - Gamma / m * x(3); return y; } int main() { int n, N; double tau, Tmax, t,pi; VectorXd x(4),k1(4),k2(4),k3(4),k4(4); pi = 4 * atan(1.0); m = 100; g = 9.8; Gamma = 1.0; e = 1.0; Tmax = 20; N = 1000; tau = Tmax / N; x << 0,0,50*cos(pi*50/180),50*sin(pi*50/180); for (n = 0; n < N; n++) { t = n * tau; k1 = tau * f(t, x); k2 = tau * f(t+tau/2, x+k1/2); k3 = tau * f(t+tau/2, x+k2/2); k4 = tau * f(t+tau, x+k3); x = x + (k1 + 2 * k2 + 2 * k3 + k4) / 6; if (x(1)<0) { x(1) = - x(1); x(3) = - x(3); } std::cout << x(0) << " " << x(1) << std::endl; } } |
最初にアクセスしたときから (入手、インストール、プログラム書いて…) 1時間しないうちに書いたプログラムなので、 とんでもないドジがあるかもしれないけれど、とにかく動いた。 すごい世の中になったものだ。
MPFRC++ Support module なんてあるのか。 これは面白いかな??
桂田 祐史