I (12/11)

情報科学概論 11 回目
2014/12/11
今回の目的
• 数値計算法 –常微分方程式の解法–
1
常微分方程式
• 常微分方程式の問題は，いつも 1 階の微分方程式の組みに置き換えて考えることができる．
d2 y
dy
(例) dx
2 + q(x) dx + r(x) = 0 は以下の連立 1 階微分方程式に書き直せる．
dy
dx
dz
dx
= z(x)
= r(x) − q(x)z(x)
x0 から x1 = x0 + h, x1 から x2 = x1 + h, · · ·, xn から xn+1 = xn + h と少しずつ x を変化させ，それに伴う y の
変化を見る．微分を差分にして解いていくが，この差分化にいろいろな方法がある．
dy
dx = f (x, y) の微分方程式の差分化を例に説明する．
• オイラー (Euler) 法 (最も単純な方法)
xn+1 = xn + h とした時，
yn+1 = yn + hf (xn , yn ) を用いて解を求めていく．h が十分に小さければ良いが，h2 のオーダーの誤差が
累積する．
• 2 次のルンゲクッタ (Runge Gutta) 法 (中点法)
区間の中点の座標を用いてステップを計算する．
同じく xn+1 = xn + h とした時，
yn+1 = yn + hf (xn + 12 h, yn + 12 hf (xn , yn ))
h3 の誤差が累積する．
• 4 次のルンゲクッタ法
同じく xn+1 = xn + h とした時，
k1 = hf (xn , yn ), k2 = hf (xn + 12 h + yn + 12 k1 )
k3 = hf (xn + 21 h + yn + 21 k2 ), k4 = hf (xn + h, yn + k3 ) とすると
yn+1 = yn + k61 + k32 + k33 + k64
h5 の誤差が累積する．
2
初期値問題 –ボールの軌道–
z 方向に一定の重力加速度 −g を受けるボールの運動を考える．ボールの位置 r の運動方程式は
d2 r
= −gˆ
z
dt2
(1)
ˆ は z 方向の単位ベクトルである．
但し，z
初期にある速度で投げ上げた場合のボールの運動を考えよう．初期に vy = 0 であれば以降もずっと vy = 0 で
あるので，運動は x − z 面内に限定できる．するとこの運動は，以下のような 1 階の連立微分方程式とすること
ができる．
dx
= vx
(2)
dt
dz
= vz
(3)
dt
1
z
mg
0
x
Figure 1: 一定重力加速度中のボールの軌道
dvx
=0
dt
dvz
= −g
dt
初期に原点から仰角 30 度で投げ上げたボールの運動を求めるプログラムは，例えば以下のようになる．
#include <stdio.h>
#include <math.h>
main()
{
float dt, t;
float x, z; /* 座標変数の宣言*/
float vx, vz; /* 速度変数の宣言*/
float ax, az; /* 加速度変数の宣言*/
dt = 0.01; /* 時間刻みに値を代入*/
x = 0.0; /* 初期値代入*/
z = 0.0;
vx = 20.0 * cos(2.0*M_PI*30.0/360.0);
vz = 20.0 * sin(2.0*M_PI*30.0/360.0);
printf("%f %f\n",x, z);
ax = 0.0;
az = -9.8;
for(t=dt;t<=2.0;t=t+dt){
x = x + vx * dt;
z = z + vz * dt;
vx = vx + ax * dt;
vz = vz + az * dt;
printf("%f %f\n",x, z); /*2 列のデータとして書き出す*/
}
return 0;
}
2
(4)
(5)
2.1
数値計算結果のアニメーション
ボールの運動を gnuplot を使用して動画表示する方法を紹介する．ボールの位置が 2 列のデータとして ball.dat
というファイルに書き出されているとする (上記のプログラムを実行し，実行結果をリダイレトでファイルに書き
出せば良い)．以下のような gnuplot のコマンドスクリプトを作成し，sample.gp として保存する (ファイル名は
好きなように付けて良い)．
if (exist("i")==0 || i<0) i=0 #変数の初期化
set pointsize 3
pl[0:40][0:10]’ball.dat’every ::i::i
pl[0:40][0:10]’ball.dat’every ::0::i with line
set xl "distance"
set yl "height"
i = i+1
pause 0.01
if(i < 207)reread #file の行数
i=-1
そして gnuplot を立ち上げ
gnuplot> load"sample.gp"
とすればアニメーションが表示できる．スクリプト中 every ::i::i はデータファイルの i 行目のみをプロット
するという意味である．every ::0::i は最初から i 行目までプロットするという意味である．
またファイルの行数 (例では 207 行) は以下のように wc コマンドを使用すれば調べられる．
% wc ball.dat
207 414 3876 ball.dat
表示される結果はそれぞれ行数，単語数，バイト数である．
3
境界値問題
前節の斜方投射の問題では，ある投射速度と角度に対しボールがどのような軌跡を描き，飛距離がどうなるかを
調べた．つまり初期値問題であった．今度は，ある飛距離を出すために必要な初速度や投げ上げ角度を計算する
ことを考える．すなわち，境界値問題として微分方程式を解こうという訳である．
x,z 成分をそのまま残したままでも計算可能であるが，ここでは水平 (x) 方向には速度 vx,0 の等速度運動をす
ることを利用し，
dx = vx,0 dt
(6)
を用いて独立変数を t から x に変換する．こうすると，以下のように解くべき方程式を 4 本から 2 本に減らすこ
とができる．(前節の初期値問題でも，このような取り扱いは可能である．)
vz
dz
=
dx
vx,0
(7)
g
dvz
=−
dx
vx,0
(8)
ここでは簡単のため，投げ上げ角度を 30 度に固定し，xm m の飛距離を出すための初速度を求めてみよう (角度も
変更する場合は，演習課題で取り組む)．まず初期の速度の推測値を v とし，その時の飛距離を x とする．当然こ
の x は一般的には xm とは異なるので，ここから解を改良し求める解に近付けていく．解の改良を δv とすると，
Taylor 展開の 1 次までを考えると (Newton-Raphson 法と同じ考えである)
(x − xm ) + δv
∂x
=0
∂v
(9)
から δv を求めれば良いことが分かる．ここで問題となるのが ∂x
∂v であるが，ここでは以下のように数値的に求め
る．そして，求めるべき初速度を逐次的に求める手法を用いる．
3
z
v0 v2
v1
0
x1
x2
xm
x0
x
Figure 2: 飛距離を決めた際に，与えるべき初速度をシューティング法で求める．
n 回の改良後の初期速度を vn とし，その時の飛距離を xn とすると
∂x
xn − xn−1
=
.
∂v
vn − vn−1
すなわち，次の解の改良値を
vn+1 = vn + δv = vn −
vn − vn−1
(xn − xm )
xn − xn−1
(10)
(11)
とすれば良い．以上のような方法をシューティング法と呼ぶ．例えば，以下のようなプログラムを作成すれば良い．
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define ACC 1.0e-12
int main()
{
double x0; /*飛距離*/
double x1,v1,x2,v2,dv;
char buf[30];
int i=0;
double integ(double); /*double 型の関数 integ を double 型の引数で使用する宣言*/
printf("distance ?");
fgets(buf,sizeof(buf),stdin); /*飛距離の入力*/
if (sscanf(buf,"%lf",&x0) != 1) {
fputs("input error\n", stderr);
exit(-1);
}
v1 = 20.0; /*initial guess*/
x1 = integ(v1);
v2 = v1*1.01;/*initial guess を少しずらす*/
x2 = integ(v2);
do{
i++;
4
dv = - (v2 - v1) / (x2 - x1) * (x2 - x0); /*解の改良*/
v1 = v2; x1 = x2;
v2 = v2 + dv;
x2 = integ(v2);
printf("i=%d,v=%25.16f x=%25.16f\n", i, v2, x2);
}while((x2-x0) >= ACC || (x0-x2) >= ACC);
return 0;
}
double integ(double v)
{
double dx, x, z, zp; /* 座標変数の宣言*/
double vx0, vz; /* 速度変数の宣言*/
double az; /* 加速度変数の宣言*/
if (v <= 0.0) { /*速度が 0 以下の場合はエラーとする*/
fputs("input error\n", stderr);
exit(-1);
}
z = 0.0;
x = 0.0;
vx0 = v * cos(2.0*M_PI*30.0/360.0);
vz = v * sin(2.0*M_PI*30.0/360.0);
dx = 0.01*vx0; /*刻み幅 (左の例では dt=0.01 に対応)*/
az = -9.8;
for(x=dx;z>=0.0;x=x+dx){
zp = z; /*前ステップの z を zp として格納しておく*/
z = z + vz/vx0 * dx; /*微分方程式の数値積分*/
vz = vz + az/vx0 * dx; /*微分方程式の数値積分*/
}
if(zp!=z){/*一般に z<0 となってるので，z=0 の時の x を内挿で求める*/
x = x - dx * z / (z-zp);
}
return x; /*最終的な飛距離を main 関数に戻す*/
}
5

Download Report