// dirty_main.cpp : Numerical Recipes の rk4.c を使って // 2次振動系の微分方程式を解く(「杉戸」問題の正しいアプローチ) // 2002/04/23 by SK // 注意:Numerical Recipesのソースを1行も書き直さないで, // C のソースをC++プログラムとして使い倒し, // さらにディフォルトの浮動小数 float を,私の好きな double として使う。 // このため,以下では汚い技を2つ使うからね! // MS-Visual C++ 使用者へ:注意その1 // メニューの ツール>オプション>ディレクトリのところで // 表示するディレクトリ=「インクルードファイル」の欄に // 次の2つのパスを書き込んでおくこと。 // ヘッダファイル(include)の在り処: // E:\SKON\C\NUMRECIPI\INCLUDE(近藤の場合は) // C のソースと nrutil.c プログラムの在り処: // E:\SKON\C\NUMRECIPI\SRC(近藤の場合は) // 上の2つ,どちらも「インクルードファイル」の欄だからね! // これで,Visual C++ が勝手に見つけてくれる。 // float を double に読み替える。(あまりにも強引,きったねー!その1) #define float double // 次の3行は,ほとんどいつも必要!! #include "nr.h" #include "nrutil.h" #include "nrutil.c" // 使いたい Numerical Recipes の C ソースを取り込む! #include "rk4.c" // このファイル dirty_main.cpp の拡張子が .cpp だから, // 取り込んでしまえばコンパイラには C++ プログラムに見える。 // だまし討ち。(あまりにも強引,きったねー!その2) // MS-Visual C++ 使用者へ:注意その2 // ワークスペースの外部依存関係の中に // nrutil.c や rk4.c // が見つかるが,Source Files へ絶対に移動してはいけない!!! // (もしそうすると,二重定義でコンパイラに叱られる。) // 一方,***.h ファイルは,Header Files へ移動してもいいけど, // 外部依存関係の中に入れておいたほうがよろしい!(Visual C++ が面倒みてくれるから。) // ******* 以上で Numerical Recipes を使う準備ができた!! ******** // ここから,杉戸問題を解こう! つまり,2次の微分方程式を解く // cos, sin の振動をシミュレーションし,ルンゲ・クッタの精度を確認する。 #define PI 3.14159265358979323846 // cos, sin 振動する系の微分方程式の右辺を準備する。周期 = 1.0 sec void func_sin(double t, double y[], double dydt[]) { dydt[1] = - 2.0 * PI * y[2]; // dy1/dt = d/dt(cos) = -omega*sin, frequency = 1 Hz dydt[2] = 2.0 * PI * y[1]; // dy2/dt = d/dt(sin) = +omega*cos, } #include #include // exit() #include // sin(), cos() <== 計算誤差を評価するため! int main(void) { double t=0.0, h=0.01/1.0, tstart=0.0, towari=5.0; FILE *d_2ndA2; double *x, *dxdt; int Norder = 2; x = vector(1,Norder); // rk4() の引数にあわせて,2次元ベクトルを確保する。 dxdt = vector(1,Norder); x[1] = 1.0; x[2] = 0.0; // x[1] = cos(0), x[2]=sin(0) 状態変数の初期値 if( (d_2ndA2 = fopen("d_2ndA2.dat","w")) == NULL) { printf("\nerror: ファイル d_2ndA2.dat が開けません\n"); exit(110); } printf("\ntime = %f x1 = %f x2 = %f :err %g %g", t, x[1], x[2], x[1]-cos(2.0*PI*t), x[2]-sin(2.0*PI*t) ); fprintf(d_2ndA2, "t x1 x2 x1err x2err \n"); fprintf(d_2ndA2, "%f %f %f %g %g \n", t, x[1], x[2], x[1]-cos(2.0*PI*t), x[2]-sin(2.0*PI*t) ); while (t < towari) { // Numerical Recipes の rk4() は次の2行をペアにして使うこと。 func_sin(t, x, dxdt); // rk4 の第1ステップの準備 rk4(x, dxdt, 2, t, h, x, func_sin); // 上の2行はペア(可変刻みの時に備えてこうなっている。) t = t + h; printf("\ntime = %f x1 = %f x2 = %f :err %g %g", t, x[1], x[2], x[1]-cos(2.0*PI*t), x[2]-sin(2.0*PI*t) ); fprintf(d_2ndA2, "%f %f %f %g %g \n", t, x[1], x[2], x[1]-cos(2.0*PI*t), x[2]-sin(2.0*PI*t) ); } printf("\n"); fclose(d_2ndA2); free_vector(dxdt, 1, 2); // vectorを使ったら,後片付けすること!!!! free_vector(x, 1, 2); return 0; // おわり! } // 刻み幅h と 計算誤差の関係! 違うプログラムでやったので異なるかもしれない。 // 刻み幅 h t=10[s] での誤差(x1, x2 の大きいほう) 有効桁数 // 0.01*20 -0.742036 振幅=0.257964 使用不可 0 // 0.01*10 -0.0675592 振幅=0.957799 減衰大きい 1.5 // 0.01*8 -0.0300467 振幅=0.986073 減衰 2 // 0.01*4 -0.00204127 振幅=0.999564 3.5 // 0.01*2 -0.000129831 振幅=0.999986 4 // 0.01 -8.14902e-6 振幅=1.000000 5 // 0.01/2 -5.09854e-7 振幅= 6.5 // 0.01/4 -3.18752E-8 振幅= 7 // 0.01/8 -1.99501E-9 振幅= 8 // 0.01/16 -1.24307E-10 振幅= 9 // 0.01/32 2.08992e-11 振幅= 10 // end of file