微分方程式の数値解法

微分方程式の数値解法
微分方程式の例
毎秒 a [m3/s] の水を注入
a
dx(t)
a
dt
体積 x(t) [m3]
解
高さx(t) [m]
x(t)  x0  at
x
底面積 1 m2
微少時間 Dt の間に注入
される量
aDt
Dx  aDt
傾き a
x0
0
t
微分方程式の例(2)
dx(t)
 bx(t)
dt
bx(t)
流出速度が x(t) に
比例するとする
微少時間 Dt の間に流出
する量 bx(t)Dt
解
x(t)  x0 ebt
x
x0
Dx  bx(t)Dt
0
t
微分方程式の例(3)
a
dx(t)
 a  bx(t)
dt
bx(t)
流入・流出が
両方ある場合
おそらくこんな解になるだろう
x
a
b
x0
0
t
厳密に解ける
dx(t)
 a  bx(t)
dt
厳密解：
(a, b は正の定数)
dx(t)
a

 b  x(t)  
右辺を

dt
b
と変形
a
y(t)  x(t) 
と置くと
b
これは直ぐ解けて
よって
dy(t)
 by(t)
dt
a  bt

bt
y(t)  y(0)e   x0   e

b
a  bt a

x(t)   x0   e 

b
b
オイラー法
たとえば
dx(t)
 a  bx(t) という微分方程式を
dt
t  0 で x  x0 という初期条件で解きたいとする。
プログラムで全ての時刻 t での x(t) を求めるのは無理
飛び飛びの時刻で x(t) を求める
x
x0
0
Dt 2Dt 3Dt
t
オイラー法(2)
dx(t)
 f x(t) の初期値問題を解く
1階の常微分方程式
dt
dx(t)
1 d2x
2
テイラー展開 x(t  Dt)  x(t) 
Dt 
(Dt)
L
2
dt
2! dt
Dt が十分小さければ無視出来るだろう
dx(t)
x(t  Dt)  x(t) 
Dt  x(t)  f x(t)Dt
dt
今の時刻 t での x(t)から、
次の時刻 t  Dt
での
x(t が求まる
 Dt)

時刻 t での速度 f x(t)
のまま時刻
tまで等速運動
 Dt
オイラー法(3)
現在時刻の x(t)が正確にわかっているとして…
x(t) のグラフ
x
dx(t)
x(t) 
Dt
dt
x(t)
dx(t )
接線の傾き
dt
x(t  Dt)
Dt
t
t  Dt
dx(t)
Dt
dt
t
オイラー法で解いてみる
オイラー法
今の例では
よって
x(t  Dt)  x(t)  f x(t)Dt
dx(t)
 a  bx(t) つまり f (x)  a  bx
dt
x(t  Dt)  x(t)  (a  bx(t))Dt
a=1, b=1, x(0)=0 の場合を解いてみる。
オイラー法
厳密解
x(t  Dt)  x(t)  (1 x(t))Dt
x(t)  1  et
#include <stdio.h>
#include <math.h>
int main()
{
double x, t, dt, tmax;
dt = 0.1; tmax = 10.0;
t = 0.0; x = 0.0;
/* 初期値 */
printf("%f\t%f\t%f\n",t,x,1-exp(-t));
while(t < tmax) {
x += (1.0-x)*dt;
/* 次の時刻の x */
t += dt;
printf("%f\t%f\t%f\n",t,x,1-exp(-t));
}
return 0;
}
結果をgnuplot で表示してみる
リダイレクトにより結果を x.txt に保存：
$ ./a.out > x.txt
$ less x.txt
0.000000
0.000000
0.200000
0.200000
..
..
時刻 t
オイラー法
によるx(t)
0.000000
0.181269
..
正確なx(t)
less を終了するには
q を押す
gnuplot でグラフ化してみる：
$ gnuplot
gnuplot> plot 'x.txt' using 1:2, 'x.txt' using 1:3
gnuplot> quit
$
演習
dt を変化させるとt=2.0 での誤差がどのように変
化するかを調べる。
誤差
？
dt
差分
差分
dx(t )
微分係数
を
dt
飛び飛びの時刻での x(t) から近似的に求める
前進差分
dx(t) x(t  Dt)  x(t)

dt
Dt
微分方程式
dx(t)
 f x(t)
dt
x(t  Dt)  x(t)  f (x(t))Dt
オイラー法
x(t  Dt)  x(t)
 f x(t)
Dt
前進差分
テイラー展開
2
dx(t)
1d x
2
x(t  Dt)  x(t) 
Dt 
(Dt)
L
2
dt
2! dt
2
dx(t)
1d x
2
Dt  x(t  Dt)  x(t) 
(Dt)
L
2
dt
2! dt
dx(t) x(t  Dt)  x(t) 1 d x


Dt  L
2
dt
Dt
2! dt
O(Dt) 高々 Dtに比例して
2
dx(t) x(t  Dt)  x(t)

dt
Dt
ゼロに近づく量
各種の差分
前進差分
dx(t) x(t  Dt)  x(t)

 O(Dt)
dt
Dt
後退差分
dx(t) x(t)  x(t  Dt)

 O(Dt)
dt
Dt
中心差分
dx(t) x(t  Dt)  x(t  Dt)

 O (Dt)2
dt
2Dt


必ず精度が高くなるわけではない
中心差分
2
3
dx(t)
1d x
1d x
2
3
x(t  Dt)  x(t) 
Dt 
(Dt)

(Dt)
L
2
3
dt
2! dt
3! dt
3
dx(t)
1 d2x
1
d
x
2
3
x(t  Dt)  x(t) 
Dt 
(Dt)

(Dt)
L
2
3
dt
2! dt
3! dt
上式から下式を引くと
3
dx(t)
2d x
3
x(t  Dt)  x(t  Dt)  2
Dt 
(Dt)
L
3
dt
3! dt
dx(t) x(t  Dt)  x(t  Dt)

 O (Dt)2
dt
2Dt


微分係数とグラフの傾き
x(t)  x(t  Dt)
後退差分
Dt
x
x(t  Dt)  x(t)
前進差分
Dt
中心差分
x(t  Dt)  x(t  Dt)
2Dt
x(t) のグラフ
dx(t )
接線の傾き
dt
x(t  Dt)
x(t)
x(t  Dt)
t  Dt
t
t  Dt
t
運動方程式を解く
運動方程式
d 2 x(t)
m
 F(t)
2
dt
質点の位置 x(t) が解れば質点に働く力 F(t) が解る場合を考える
例えば、バネに繋がれた質点の場合
0
自然の長さ
x(t)
x(t)
F(t)  kx(t)
x
オイラー法で解く
dx(t)
速度 v(t) 
を導入して１階の連立常微分方程式にする：
dt
dx(t)
 v(t)
dt
dv(t) 1
 F(t)
dt
m
前進差分
x(t  Dt)  x(t)
 v(t)
Dt
v(t  Dt)  v(t) 1
 F(t)
Dt
m
現在の x(t), v(t), F(t) から
x(t  Dt)  x(t)  v(t)Dt
次の時刻の x(t  Dt), v(t  Dt)
1
v(t  Dt)  v(t)  F(t)Dt
m
が計算できる
厳密解
d 2 x(t)
 kx(t)
運動方程式 m
2
dt
両辺を m で割って
一般解
d 2 x(t)
2



x(t)
2
dt
但し

x(t)  C1 cos  t  C2 sin  t
初期条件を
t = 0 で x = x0 , v = 0
x(t)  x0 cos  t
とすると
k
m
#include <stdio.h>
#include <math.h>
int main()
{
double k, m, w, x, v, t, dt, tmax, F;
k = 1.0; m = 1.0; w = sqrt(k/m); /* パラメタ */
tmax = 10.0; dt = 0.02;
/* 計算条件 */
t = 0.0; x = 1.0; v = 0.0;
/* 初期条件 */
while( t < tmax ) {
printf("%f¥t%f¥t%f¥n", t, x, cos(w*t));
F = -k*x;
/* 質点に働く力 */
x += v*dt;
/* 次の時刻の x */
v += F/m*dt;
/* 次の時刻の v */
t += dt;
}
return 0;
}
実行してみる
リダイレクトにより結果を x.txt に保存：
$ cc oscillator.c -lm
$ ./a.out > x.txt
$ less x.txt
0.000000
1.000000
0.020000
1.000000
..
..
時刻 t
オイラー法
によるx(t)
1.000000
0.999800
..
less を終了するには
q を押す
正確なx(t)
gnuplot でグラフ化してみる：
$ gnuplot
gnuplot> plot 'x.txt' using 1:2, 'x.txt' using 1:3
gnuplot> quit
$
ベルレ法
dx(t)
1 d 2 x(t)
2
x(t  Dt)  x(t) 
Dt 
(Dt)

2
dt
2! dt
dx(t)
1 d 2 x(t)
2
x(t  Dt)  x(t) 
Dt 
(Dt)

2
dt
2! dt
1 d 3 x(t)
3
(Dt)
L
3
3! dt
1 d 3 x(t)
3
(Dt)
L
3
3! dt
この２式を加えると
d 2 x(t)
2
x(t  Dt)  x(t  Dt)  2x(t) 
(Dt)
L
2
dt
よって
d 2 x(t) x(t  Dt)  2x(t)  x(t  Dt)

2
2
dt
Dt 
F(t)
x(t  Dt)  2x(t)  x(t  Dt) 
(Dt)2
m
次の時刻
現在
一つ前
の時刻
現在
2階微分に対
する差分
ベルレ法
ベルレ法(2)：別の見方
dx(t)
v(t) 
dt
であるから、中心差分により
d 2 x(t) dv(t) v(t  Dt / 2)  v(t  Dt / 2)


2
dt
dt
Dt
再度中心差分を用いると
x(t  Dt)  x(t)
v(t  Dt / 2) 
Dt
x(t)  x(t  Dt)
v(t  Dt / 2) 
Dt
よって
Dt
t
2
t  Dt
Dt
t
2
t
t  Dt
Dt
d 2 x(t) 1  x(t  Dt)  x(t) x(t)  x(t  Dt) 
 


2
dt
Dt 
Dt
Dt
x(t  Dt)  2x(t)  x(t  Dt)
前ページと同じ式

2
Dt 
ベルレ法(3)
F(t)
x(t  Dt)  2x(t)  x(t  Dt) 
(Dt)2
m
例えば
F(t)  kx(t)
●初期条件として x(0) に加え x(-Dt) が必要
例えば x(Dt)  x(0)  v(0)Dt
1 F(0)
2
Dt 
あるいは x(Dt)  x(0)  v(0)Dt 
2 m
●速度が必要な場合、 x(t  Dt) を求めた後
x(t  Dt)  x(t  Dt)
v(t) 
2Dt
で現在の速度を求める。
次の時刻での位置が解って初めて現在の速度が解る
（この不便を解消する方法もあるが、やや煩雑なのでここ
では単純なベルレ法を使う）
#include <stdio.h>
#include <math.h>
int main()
{
double k, m, w, xold, xnow, xnew, t, dt, tmax, F, E;
k = 1.0; m = 1.0; w = sqrt(k/m);
tmax = 10.0; dt = 0.02;
t = 0.0; xnow = 1.0;
xold = xnow;
/* t=0で速度v=0だから、xoldの予測値=xnow */
while( t < tmax ) {
F = -k*xnow;
xnew = 2*xnow - xold + F/m*dt*dt; /* 次の時刻の x */
vnow = (xnew - xold)/(2*dt);
/* 現在の速度 */
E = 0.5*(m*vnow*vnow + k*xnow*xnow); /* エネルギー */
printf("%f¥t%f¥t%f\t%f¥n",t, xnow, cos(w*t), E);
t += dt;
xold = xnow;
xnow = xnew;
}
}
エネルギー保存則
力学的エネルギー
1
1
2
E  mv(t)  kx(t)2
2
2
運動エネルギー
ポテンシャルエネルギー
E は時間に依らない一定値になるはず
x(t)  x0 cos  t の場合
v(t)   x0 sin t
1
1 2
2 2
2
E  m x0 sin  t  kx0 cos2  t
2
2
1 2
1 2
2
2
 kx0 sin  t  cos  t  kx0
2
2

k
m
演習(今回はこれだけ)
オイラー法でエネルギー保存則がどれくらい良く
満たされているかについて
•dt を変化させて誤差がどう変化するか
E
オイラー
？
log誤差
？
t
log dt
演習
エネルギー保存則がどれくらい良く満たされて
いるかについて
•オイラー法とベルレ法で比較
•dt を変化させて誤差がどう変化するか
E
オイラー
？
ベルレ
log誤差
？
t
log dt
放物運動
d2x
m 2 0
dt
d2y
m 2  mg
dt
初期条件
x(t  0)  0,
y(t  0)  0
厳密解
x(t)  vx 0t
1 2
y(t)  vy0 t  gt
2
y
v0
q
0
x
vx (t  0)  vx0  v0 cosq ,
vy (t  0)  vy0  v0 sinq
オイラー法で解く
x(t  Dt)  x(t)  vx (t)Dt
vx (t  Dt)  vx (t)
y(t  Dt)  y(t)  vy (t)Dt
vy (t  Dt)  vy (t)  gDt
#include <stdio.h>
#include <math.h>
int main()
{
double x, vx, y, vy, t, dt;
double m, g, v0, theta, vx0, vy0, deg;
deg = 4.0*atan(1.0)/180.0;
/* 度をラジアンに変換 */
g = 9.8; m = 0.1;
/* 質量と重力加速度 */
v0 = 10.0; theta = 30.0*deg; /* 初速と投げる方向 */
vx0 = v0*cos(theta); vy0 = v0*sin(theta);
dt = 0.01;
t = 0.0; x = 0.0; y = 0.0;
vx = vx0; vy = vy0;
(次ページに続く)
tan

4
1
tan 1 1 

4
/* 初期条件 */
/* 初期条件 */
while(y >= 0.0) {
printf("%f¥t%f¥t%f¥t%f¥t%f¥n",
t, x, y, vx0*t, vy0*t - 0.5*g*t*t);
x += vx*dt;
vx += 0;
y += vy*dt;
vy += -g*dt;
t += dt;
}
return 0;
}
地面に落ちてくるまで計算する
空気抵抗
速度 v で空気中を運動する物体に働く抵抗力 F は、およそ
1
F   CD air Svv
2
となることが知られている。
(速度と逆向きで、速度の2乗に比例する）
但し  air =空気の密度=約 1 kg/m3
C D =抵抗係数(物体の形状に依存)
S =物体の断面積
v  v  vx2  vy2
密度 ball 、半径 a の球の場合 C D  0.5 程度であり
S  a
2
m  ball
air 1
k  0.2
ball a
4 3
a
3
とおくと
F
air 1
 0.2
vv
m
ball a
F
 kvv
m
空気抵抗
 
 
抵抗力 F の x,y 成分を Fx, Fy とすると、運動方程式は
d2x
m 2  Fx
dt
d2y
m 2  mg  Fy
dt
両辺を m で
割って
d2x
 kvvx
2
dt
d2y
 g  kvvy
2
dt
オイラー法を使うとすると、例えば y は
dy(t)
y(t  Dt)  y(t)
左辺を
 vy (t)
 vy (t)
dt
Dt
前進差分
dvy (t)
vy (t  Dt)  vy (t)
で近似
 g  kv(t)vy (t)
 g  kv(t)vy (t)
dt
Dt
よって

y(t  Dt)  y(t)  vy (t)Dt


vy (t  Dt)  vy (t)  g  kv(t)vy (t) Dt
但し
v(t)  vx (t)2  vy (t)2
#include <stdio.h>
#include <math.h>
int main()
{
double x, vx, y, vy, t, dt;
double g, v0, theta, vx0, vy0, deg;
double k, v;
deg = 4.0*atan(1.0)/180.0;
/* 度をラジアンに変換 */
g = 9.8;
/* 重力加速度 */
k = 0.2/(0.01*1000);
/* 0.2/(半径*密度) */
v0 = 10.0; theta = 30.0*deg; /* 初速と投げる方向 */
vx0 = v0*cos(theta); vy0 = v0*sin(theta);
dt = 0.01;
t = 0.0; x = 0.0; y = 0.0;
vx = vx0; vy = vy0;
(次ページに続く)
/* 初期条件 */
/* 初期条件 */
while(y >= 0.0) {
printf("%f¥t%f¥t%f¥t%f¥t%f¥n",
t, x, y, vx0*t, vy0*t - 0.5*g*t*t);
v = sqrt(vx*vx + vy*vy); /* 現在の速さ */
x += vx*dt;
vx += -k*v*vx*dt;
y += vy*dt;
vy += -g*dt - k*v*vy*dt;
t += dt;
}
return 0;
}
dt は十分小さくすること

Download Report