C言語基礎まとめ
プログラミング言語Ｃ
以降、K & Rと略記します。
常微分方程式は、物理の様々な分野で登場します。 例えば、力学のニュートンの運動方程式は、 物体の運動を決定するための常微分方程式です。 常微分方程式を解くための数値解法が多数開発されています。
今日の演習では、
微分方程式の近似として、差分をあつかう。
しかし、元々の方程式じたい近似ではないのか。
数値計算じたいの諸問題。
以下の参考書の項、洲之内(第１章)、および伊理・藤野を参照
のこと。
洲之内治男『数値計算[新訂版]』(サイエンス社、2002年)。
第4章。
伊理正夫・藤野和建『数値計算の常識』(共立出版、1985年)、
第12章。
放物運動を考察する。下図のようにx軸、y軸をとると、
y軸方向に重力のみが働く物体のニュートンの運動方程式は、
となる。
重力加速度をg=9.8 m/s2とし、質量を m=1 Kgとした方程式を考察する。
この常微分方程式をオイラー法で解くために、１階常微分方程式に書き換えよう。
そのために、
を導入する。
運動方程式は x,y,v,uを用いて、１階の微分方程式に
書きなおすことができる：
オイラー法では、
t=0,dt,2dt,...の時刻におけるx,y,v,uの値
xi,
yi,
ui,
vi
を
の漸化式に従って計算する。
物体は初期時刻 t=0 で原点から速さ 10 m/sで
x 軸から角度 Θ°の方向に投げ出されたとする。
このとき、x(0)=0, y(0)=0,
v(0)=10 cos(πΘ/180),
u(0)=10sin(πΘ/180) である。
以下のソースコード(dynamics1.c)は、
この運動をオイラー法によって解くものである。
#include <stdio.h>
#include <math.h>
#define PI 3.141592
int main() {
double t,dt; /* t: 時刻(秒)、dt: 時間刻み(秒) */
double t1; /* t1: 終了の時刻 */
double x, y, xa, ya; /* x, y: 旧いx, y-座標。xa, ya: 新しいx, y-座標。*/
double v, u, va, ua; /* va, ua, u, v, はそれぞれ新旧の x, y の速度を表す変数 */
double th; /* 投げ上げの角度θ(度) */
printf("投げ上げの角度(度)、 終了の時刻(秒)、\n");
printf("時間刻み(秒)を空白を置いて入力してください。\n");
scanf("%lf %lf %lf", &th, &t1, &dt);
x=0.0; /* 初期位置と初期速度の設定 */
y=0.0;
v=10.0*cos(PI/180*th);
u=10.0*sin(PI/180*th);
for(t=0.0; t<=t1; t= t + dt){
printf("%f %f %f\n", t, x, y);
xa=x+v*dt;
ya=y+u*dt;
va=v;
ua=u-9.8*dt;
x=xa;
y=ya;
v=va;
u=ua;
}
return 0;
}
このプログラムには、実行時に、３つの数値をあたえます。 物体の投げ出される角度Θ 、終了時刻、 時間の刻み幅です。これら三つの数値を半角の空白を 挟んで半角で入力します。
$ g++ dynamics1.c あるいは $ gcc dynamics1.c -lm (gccで数学関数の入ったファイルをあつかう場合、 オプション「-lm」が要ります。 $ ./a.out 投げ上げの角度(度)、 終了の時刻(秒)、 時間刻み(秒)を空白を置いて入力してください。 45 1.5 0.05 0.000000 0.000000 0.000000 0.050000 0.353553 0.353553 ..... ...... ...... ･･･角度45°で 1.5秒まで 0.05秒刻みで時刻 と物体の座標を出力する。 （t、ｘ、ｙ座標が次々と出力される） $ファイルに出力するときは、
$ ./a.out > d0.dat 45 1.5 0.05とすれば、ファイル「d0.dat」に出力される。
$ cat d0.dat として、「Enter」すると、端末に 投げ上げの角度(度)、 終了の時刻(秒)、 時間刻み(秒)を空白を置いて入力してください。 0.000000 0.000000 0.000000 0.050000 0.353553 0.353553 ...... ...... ...... などと出る。
放物運動の解は、解析的に求めることもできる。
これより、軌跡の方程式 y(x) を導出せよ。
オイラー法の計算が正しく行えているかを調べるために、
時間刻みdtを小さくすると、オイラー法の計算結果が解析解の軌跡に近づいていくことを確認しよう。
このため、
時間刻みdtが
0.125,
0.015625,
0.001953125,
0.000244140625
の場合について
角度45°で打ち出された物体の運動を1.5秒まで計算する。
これらの出力データをファイルに出力しよう。
gnuplot を用いて
課題１のy(x)のグラフとともに、
すべての出力データを重ねてプロットすることで、
dtを小さくするとどんな様子でオイラー法が正確になっていくか大まかに把握しよう。
出力データファイル「d0.dat」のデータの2列目をx軸、 3列目をy軸に取って、グラフに表わすには、
gnuplot> plot "d0.dat" using 2:3とすればよいでしょう。
時間刻みdtが
0.125,
0.015625,
0.001953125,
0.000244140625
の出力データファイルから
t=1.5秒のときのx,y座標の値を読み取る。
さらに、
解析的に求められたx,y座標の値からのずれΔを、それぞれのdtについて求めよう。
オイラー法が1次の数値解法であることから
Δはdtに比例することが予想されるが、
実際にそうなっているであろうか？
gnuplotを用いて、
Δをdtの関数としてプロットすることでこれを確かめよ。
ヒント
いろんなやり方があるが、愚直にやるには、
まず、emacsをつかって、dt_xy.datという名前のファイルを開いておく。
それから、
catで0.125, 0.015625, 0.001953125, 0.000244140625 の出力データファイルを出力
して、t=1.5秒のときのx,y座標の値を端末でコピーして、
dat_xy.datにペーストする。
dtの値は手で入力すると、ファイルの内容は、
dtの値 xの値 yの値 dtの値 xの値 yの値 （以下省略）
となる。
Δはdtに比例しているかどうかは、
を確かめれば良い。gnuplotを起動して、
$ gnuplot gnuplot> set logscale xy gnuplot> plot ...
とすると両対数プロットが出来る。
dt_xy.datというファイルの１、２、３行に dt、 x座標、 y座標の値を記録したとしよう。すると、 gnuplotを用いて望みのグラフを作るにはつぎのようにすればよい。
gnuplot> v=10/sqrt(2) gnuplot> g=9.8 gnuplot> t=1.5 gnuplot> xt=v*t gnuplot> yt=v*t-g*t**2/2 gnuplot> print xt, yt 10.6066017177982 -0.418398282201791 gnuplot> set logscale xy gnuplot> plot "dt_xy.dat" using 1:(sqrt(($2-xt)**2+($3-yt)**2)) w lp, 8*x gnuplot>
解説
gnuplotでも簡単な計算ができる。上の
using 1:(sqrt(($2-xt)**2+($3-yt)**2))
のところで、ずれΔを計算している。
課題２と課題３のグラフをGIMPで開く。
ファイルメニューから印刷を選ぶ。
このようにして、印刷したものを提出すること。
クラス(A2, B2)、学籍番号、氏名を手書きで
おおきく記入すること。提出の締め切りは、最終回までです。