VisualC++ と OpenGL を利用した仮想物理実験室
【1-3-5】等加加速度直線運動の数値解と解析解との比較

文責：遠藤理平（2010年10月 1日）カテゴリ：仮想物理実験室(325)

【1-2-5】等加速度直線運動の数値解と解析解との比較とと同様に、【1-3-2】等加加速度直線運動のシミュレーションでの数値解の結果と、【1-3-4】等加加速度直線運動の解析解の導出の導出で導いた解析解とを比較するためにグラフ化を行います。

【1-3-3】等加速度直線運動のグラフ化と同様に、シミュレーションの数値解の結果をファイルに書き出します。今回は$\Delta t = 1, 0.1, 0.01, 0.01$ の４つの場合について、計算結果と解析解との比較を行ないます。また、加加速度 j =0 [m/s^3]、加速度 a =0 [m/s^2]}、初速度 v_0 = 0 [m/s]、初期位置 x_0 = 0 [m] とします。解析解は

です。 Δt = 1 のときのファイル名を「ball3_1.0.data」、 Δt = 0.1 のときのファイル名を「ball3_0.1.data」、 Δt = 0.01 のときのファイル名を「ball3_0.01.data」 Δt = 0.001 のときのファイル名を「ball3_0.001.data」とそれぞれの結果を出力するファイル名を変更します。

VisualC++ の変更

//////////////////////////////////////////////////////////////////////////
// 変数の定義
//////////////////////////////////////////////////////////////////////////
#define _BITMAP 0//アニメーション作成用ビットマップの保存 0:しない　1:する
ofstream file1( "ball3_0.001.data" );
//--------------------------------------------------------
// 仮想物理実験室変数の定義
//--------------------------------------------------------
double t = 0.0;    //時刻
double dt= 0.001;  //時間刻み
int tn = 0;        //ステップ数
int skipCal = 10;         //描画回数の間引き
int skipBMP = skipCal *10; //BMP書出し回数の間引き数
int skipDat = 100;        //データ書出し回数の間引き回数

// ボールの定義
struct BALL {  //構造体の定義
  double x, y, z;
  double vx, vy, vz;
  double ax, ay, az;
  double jx, jy, jz;
};
BALL ball1;

void SetUp(void){//初期値を設定（main関数の中で実行）
  ball1.x = -50.0;
  ball1.y = 0.0;
  ball1.z = 20.0;
  ball1.vx = 0.0;
  ball1.vy = 0.0;
  ball1.vz = 0.0;
  ball1.ax = 0.0;
  ball1.ay = 0.0;
  ball1.az = 0.0;
  ball1.jx = 1.0;
  ball1.jy = 0.0;
  ball1.jz = 0.0;
}

等加加速度直線運動の数値誤差の計算

ファイル出力の３つ目に誤差を書き出すために、「Calculate()」関数を少し変更します。

//--------------------------------------------------------
// 計算と物体の描画
//--------------------------------------------------------
void Calculate(){
  double delta_x = abs(ball1.x - 1.0/6.0 * ball1.jx * pow(t,3)); //誤差
  if(tn%skipDat == 0 )  file1 << t << " " << ball1.x << " " << delta_x << endl;
  //加速度の算出
  ball1.ax = ball1.ax + ball1.jx * dt;
  ball1.ay = ball1.ay + ball1.jy * dt;
  ball1.az = ball1.az + ball1.jz * dt;
  //速度の算出
  ball1.vx = ball1.vx + ball1.ax * dt;
  ball1.vy = ball1.vy + ball1.ay * dt;
  ball1.vz = ball1.vz + ball1.az * dt;
  //位置の算出
  ball1.x = ball1.x + ball1.vx * dt;
  ball1.y = ball1.y + ball1.vy * dt;
  ball1.z = ball1.z + ball1.vz * dt;
}

等加加速度直線運動の t-x グラフとgnuplot のテンプレート

等加加速度直線運動のグラフ化と同様に、描画用のテンプレートを用意します。

#############################################
#jpg出力の場合
#############################################
set terminal jpeg  enhanced font "Times" 20
set output 'ball3_.jpg'
#############################################
#eps出力の場合
#############################################
#set term postscript eps enhanced color "Times" 20
#set output 'ball3_.eps'

set tics font 'Times,18'
set nokey
plot "ball3_0.001.data" using 1:2 linewidth 5 ,\
     "ball3_0.01.data" using 1:2 linewidth 5 ,\
     "ball3_0.1.data" using 1:2 linewidth 5 ,\
     "ball3_1.0.data" using 1:2 linewidth 5,\
     1.0/6.0*x**3 w l

等加加速度直線運動の数値解と解析解の比較

Δt = 1 の場合の解析解とのずれ、すなわち誤差が目立ちます。しかしながらこのグラフでは、どの程度の誤差があるのかがよくわかりません。式

で定義される、数値解と解析解との差、誤差|δx(t)|をグラフ化します。

等加加速度直線運動の数値誤差|δx(t)|の時間依存グラフとgnuplot のテンプレート

#############################################
#jpg出力の場合
#############################################
set terminal jpeg  enhanced font "Times" 20
set output 'ball3__.jpg'
#############################################
#eps出力の場合
#############################################
#set term postscript eps enhanced color "Times" 20
#set output 'ball3__.eps'

set tics font 'Times,18'
set nokey
plot "ball3_0.001.data" using 1:3 linewidth 5 ,\
     "ball3_0.01.data" using 1:3 linewidth 5 ,\
     "ball3_0.1.data" using 1:3 linewidth 5 ,\
     "ball3_1.0.data" using 1:3 linewidth 5

等加加速度直線運動の数値誤差

誤差は式

で定義されている加速度と加加速度の関係のしかたによって変わります。オイラー法の場合、誤差は$\Delta t$に比例し、等加加速度直線運動では時刻$t$の2乗に比例します。

オイラー法を用いた場合の、等速度直線運動、等加速度直線運動、等加加速度直線運動の誤差をまとめると次のようになることがわかりました。