next up previous 次へ: レポートに関して 上へ: 実験1：C++/Octaveによる微分方程式の数値解法 戻る: 例題
実験：野球ボールの軌道計算
前節までで学習したオイラー法を用いて，ピッチャーが投げた野球ボールを打者が打った時のボールの軌道を計算するプログラムをC++言語を用いて作成しよう．
まず，ボールの運動に関する基礎方程式を説明する．運動を記述するには，投射物の位置ベクトル と速度ベクトル
を計算する必要がある．この運動方程式は，
(14)
となる．ただし， は投射物の質量，
は空気抵抗による力，
は重力加速度，
は鉛直方向の単位ベクトルである．実際の空気抵抗の計算は非常に複雑になるので，ここでは次の近似式を用いることとする．
(15)
ここで， は抵抗計数，
は空気の密度，
は投射物の横断面積である．
空気抵抗を無視できる場合，運動方程式は以下の方法で解析的に解くことができる．
(16)
ただし， と
は位置と速度の初期値である．
以下にプログラムの概要をまとめる．
\begin{itembox}[c]{プログラム概要} \begin{itemize} \item ボールの初期位置 $$r_{1... ...�屬鯢充┐垢襦\item ボールの軌道をグラフ表示する． \end{itemize}\end{itembox}
#include "NumMeth.h"void main() {
//* ボールの初期位置及び初期速度を設定する． double y1, speed, theta; double r1[2+1], v1[2+1], r[2+1], v[2+1], accel[2+1]; cout << "高さの初期値(メートル) : "; cin >> y1; r1[1] = 0; r1[2] = y1; // 初期位置ベクトル cout << "初期速度(m/s) : "; cin >> speed; cout << "初期角度(度) : "; cin >> theta; const double pi = 3.141592654; v1[1] = speed*cos(theta*pi/180); // 初期速度 (x) v1[2] = speed*sin(theta*pi/180); // 初期速度 (y) r[1] = r1[1]; r[2] = r1[2]; // 初期位置および初期速度を設定 v[1] = v1[1]; v[2] = v1[2];
//* 物理パラメータを設定(質量, Cd値など) double Cd = 0.35; // 空気抵抗(無次元) double area = 4.3e-3; // 投射物の横断面積(m^2） double grav = 9.81; // 重力加速度(m/s^2) double mass = 0.145; // 投射物の質量(kg) double airFlag, rho; cout << "空気抵抗(あり:1, なし:0) : "; cin >> airFlag; if( airFlag == 0 ) rho = 0; // 空気抵抗なし else rho = 1.2; // 空気の密度(kg/m^3) double air_const = -0.5*Cd*rho*area/mass; // 空気抵抗定数
//* ボールが地面に着くまで, あるいは最大の刻み数になるまでループ
double tau;
cout << "時間刻みτ(秒) : "; cin >> tau;
int iStep, maxStep = 1000; // 最大の刻み数
double *xplot, *yplot, *xNoAir, *yNoAir;
xplot = new double [maxStep + 1];
yplot = new double [maxStep + 1];
xNoAir = new double [maxStep + 1];
yNoAir = new double [maxStep + 1];
for( iStep=1; iStep<=maxStep; iStep++ ) {
//* プロット用に位置(計算値および理論値)を記録する xplot[iStep] = r[1]; // プロット用に軌道を記録 yplot[iStep] = r[2]; double t = ( iStep-1 )*tau; // 現在時刻 xNoAir[iStep] = r1[1] + v1[1]*t; // 位置(x) yNoAir[iStep] = ?????? // 位置(y)
//* ボールの加速度を計算する double normV = sqrt( v[1]*v[1] + v[2]*v[2] ); accel[1] = air_const*normV*v[1]; // 空気抵抗 accel[2] = air_const*normV*v[2]; // 空気抵抗 accel[2] -= grav; // 重力
//* オイラー法を用いて, 新しい位置および速度を計算する
/* ここにオイラー法のアルゴリズムを書く */
//* ボールが地面に着いたら(y < 0)ループを抜ける
/* if文を使って書く */
}
//* 最大到達高さと滞空時間を表示する
cout << "最大到達高さは" << r[1] << "メートル" << endl;
滞空時間の表示をここに書く
//* プロットする変数を出力する
// xplot, yplot xNoAir, yNoAir
ofstream xplotOut("xplot.txt"), yplotOut("yplot.txt"),
xNoAirOut("xNoAir.txt"), yNoAirOut("yNoAir.txt");
int i;
for( i=1; i<=iStep+1; i++ ) {
xplotOut << xplot[i] << endl;
yplotOut << yplot[i] << endl;
}
for( i=1; i<=iStep+1; i++ ) {
xNoAirOut << xNoAir[i] << endl;
yNoAirOut << yNoAir[i] << endl;
}
delete [] xplot, yplot, xNoAir, yNoAir; // メモリを開放
}
問題:3
ボールの軌道のグラフを gnuplot で出力せよ．なお，時間刻みは
の間で5個選ぶこと．補助的に Octave を用いてよもい．
問題:4
時間刻み
問題:5
オイラー法よりも高精度な数値計算アルゴリズムについて調べて よ．余力のある人は，アルゴリズムを実装しその結果をオイラー法 と比較してみよ．
Tomoki MIYAZATO 平成16年6月29日
