☆【プログラム例】シンプソン法による関数の積分 ☆【課題2-1】放物運動 ☆【課題2-2】着陸船のシミュレーション → シンプソン法については自ら確かめよ → 課題2-1または2-2を実行し、それをレポートとして提出せよ標準入出力関数ライブラリの読み込み #include <stdio.h> #include <math.h> 数学関数ライブラリの読み込み；コンパイル時に"-lm"が必要 #define N 刻み数の定義 double 30 被積分関数を与えるサブルーチン（関数）の定義 func_f( double ); main() { double double double int f[N+1]; xa=0.0, xb=3.0; z=0.0, h=0.0, x, s; i; h: 刻み幅 h=(xb-xa) / (double)N; for( i=0; i<=N; i++ ) { x = xa + h*(double)i; f[i] = func_f(x); } for( i=1; i<N; i++ ) 積分点（N+1点）における関数の値；積分区間はＮ個積分範囲の指定：下限; xa、上限; xb z; 台形公式で第２項目からN項目の関数の値の格納場所を初期化 z += 2.0*f[i]; s = (h/2.0) * (f[0]+z+f[N]); i=0が積分の左端（第1項目）、 i=Nが積分の右端（第(N+1)項目）であることに注意全部で(N+1)個の積分点のx座標を生成各積分点での被積分関数の値台形公式の第１項目と第(N+1)項目は1/2の重み第2項目からN項目まで（ループ）は1の重みで関数値を積算 printf("width_h, Answer = %8.4lf %12.8lf\n",h,s); return; 結果の出力（刻み幅と積分値） } double func_f( double x) { return pow(x,4.0) + 2.0*x; } 被積分関数を計算するサブルーチン（関数）被積分関数の形 x 4 + 2x プログラムは、三井田、須田「数値計算法」（森北出版）より放物運動:擬似コード define t,g,h,m 変数、定数の定義 x[0] = 0.0 z[0] = 0.0 初期値の設定 px[0] = 1.0 pz[0] = 2.0 t = 0 begin do-loop (while z>0) 地面より上にある条件で繰り返し t = t+h 時間の更新 x = x+h*(px/m) x座標の更新（オイラーのアルゴリズム） z座標の更新 z = z+h*(pz/m) 運動量の更新（x成分は一定） px= px + 0 pz= pz - h*(m*g) 運動量の更新 print out ; t, x, z 逐次結果を出力 #include <stdio.h> 標準入出力関数ライブラリの読み込み #define 重力加速度（9.80665 m/s2) g main() { double double double double double 9.80665 t=0.0; h=0.001; x=0.0, z=0.0; m=1.0; px=10.0, pz=10.0; printf(" t, do { x, h: 時間の刻み幅（0.001秒）初期座標（原点：地面）物体の重さ（1 kg）初期運動量（初速度のx,z成分がそれぞれ10m/sに相当） z\n"); 結果（テーブル型出力）のヘッダ do whileループ（条件文はループを終了したときにチェックする書式） printf("%6.3lf, %8.4lf, %8.4lf\n", t,x,z); t += h; x += h*(px/m); z += h*(pz/m); px += 0; pz += -h*m*g; }while(z>0); return; } 時間、x,z座標の出力（初期値から）時間の更新 x座標の更新 z座標の更新運動量（px）の更新（実際には更新していない：等速度）運動量（pz）の更新 z座標が正（地面より上）の間do-loopを計算し続ける z [m] 放物運動:計算結果 x [m] 演習課題2-1: （一応の提出期限：１１月２５日授業前）前述の放物運動について、x方向のみ、速さに比例する摩擦力を受けるとき、軌跡はどのように変化するか。（１）まず摩擦がない場合の放物運動を計算・プロットせよ（２）摩擦力を受ける場合の運動方程式を立てよ（３）アルゴリズムのどこを変更すれば良いか？これがわかればプログラムのどこを書き換えればよいかはすぐわかる。（４）適当な摩擦係数を仮定しシミュレーションしてみよう。＄難しい場合は例えば、x方向からの向かい風があるような場合を考えてもよい。実際の空気中の摩擦は大変複雑である。野球のボール程度のスピードの場合には、速度の二乗に比例する抵抗力が働く。ボールが回転している場合にはさらに複雑である（しかし、我々は、カーブやシュートの投げ方を経験的に知っている）。演習課題2-2: （一応の提出期限：１１月２５日授業前）教科書2.1.2節「着陸船のシミュレーション」のプログラムを実行し（lander.c）、様々な条件で着陸のシミュレーションを試してみよ。この際、フリーグラフソフトである、"gnuplot"の使い方を覚えることを強くお勧めする。Win版、Mac版、UNIX版いずれもフリーで配布されており、ウェブを検索すれば、インストール方法・使用方法もすぐにわかるはずだ。