微分方程式：情報科学の視点から
Introduction to Diﬀerential Equations:
from the viewpoint of Information Science
亀島鉱二
2014 年 2 月 11 日
ダイナミックな
システム
2
3
電子回路
15kΩ
j
i
e
c
100µF
f
g
h
sw
a
3V
d
• コンデンサ両端の電圧はどのように変化する?
• 抵抗を流れる電流はどのように変化する?
b
4
ある疑問
バルカン半島とイタリア半島に挟まれたアドリア海でのこ
と．第一次世界大戦 (1914–1918) の前は水産資源に恵まれ，
漁業が盛んに行われていた．戦争が始まり，漁民が兵士と
して出征したため，漁業は一休み．戦争が終わり故郷に帰っ
た漁民は，漁業を再開．長い休漁期間の後なのでさぞや魚
が増えたことだろうと期待したが魚はちっとも獲れず．そ
れどころか，人間を襲うフカが増えている始末．豊かなア
ドリア海は，戦争の間にすっかり荒涼とした海に変わっていた．
year
cartilaginous fish(%)
1914
11.9
1919
27.3
1915
21.4
1920
16.0
1916
22.1
1921
15.9
1917
21.2
1922
14.8
1918
36.4
1923
10.7
イタリア人の数学者による説明
year
cartilaginous fish(%)
1914
11.9
1919
27.3
1915
21.4
1920
16.0
1916
22.1
1921
15.9
5
1917
21.2
1922
14.8
食物連鎖のダイナミクス食われるもの（サバ）を取り除く
（漁獲する）と，餌不足で食うもの（フカ）が減少し，その
結果，食われるもの（サバ）の平衡密度は逆に増加する．
1918
36.4
1923
10.7
ダイナミックな食物連鎖
d
人口 = (成長率 − 死亡率) × 人口.
dt
成長率 − 死亡率:
> 0? = 0? < 0?
6
競合成長モデル (LV.c)
7
LV.c[lv]/ LV P GRAPH, LV L GRAPH
8
GPS
c = const.
∆Ts − ∆Te ≃ 40µs
ダイナミックなシステム
9
10
0
自然現象の数学モデル (1)
11
Galilei の実験
H
6
y(t)
6
6v(t)
−v(t) ?
- t
12
力と速度
慣性力 :
重力 :
位置 :
f = ma,
f = −mg,
g = 9.8m/s2
y(t),
dy(t)
速度 :
v(t) =
,
dt
(
)
dv(t)
d
dy(t)
加速度 :
a=
=
,
dt
dt
dt
d2 y(t)
ma = m
= f = −mg,
2
dt
d2 y(t)
微分方程式 :
= −g.
2
dt
(1)
(2)
微分方程式による法則の表現
(
• 微係数
(
• 方程式
dy d2 y
,
dt
dt2
)
を左辺に含む
d2 y
= −g
2
dt
)
(
• t = 0 のときの位置
(
が定まっている
13
)
= y(0)
)
初期条件 .
(
と速度
dy(0)
=
dt
)
微分方程式を解く
微分方程式をみたす関数
14
y(t) はシステムの振舞をあらわす．
(
d2 y(t)
= −g
• 方程式
2
dt
関数 y(t) を解といい，
)
をみたす
• 方程式を満たす関数を求めることを
微分方程式を解く
という．
d2 y
t を省略して方程式
= −g の解を y とかく．
2
dt
微分と積分
d
dt
(∫
)
y(t)dt
∫
15
dy
dt
dt
= y(t).
= y(t).
(1)
(2)
落下物体：微分方程式
2
d y(t)
dt2
= −g,
1 回積分すると
...
2 回積分すると
...
(
)
dy(0)
y(0) = H,
=0
dt
16
d2 y
微分方程式
=
−g
dt2
∫
2
d y(t)
dy(t)
dt =
2
dt
dt
∫
17
∫
= −
gdt
= −gt + C1 ,
(1)
∫
dy(t)
dt = y(t) =
(−gt + C1 ) dt
dt
gt2
= −
+ C1 t + C2 . (2)
2
微分方程式を解くことを “微分方程式を積分する” ともいう．
微分方程式の解：検算
dy(t)
dt
d2 y(t)
dt2
gt2
y(t) = −
+ C1 t + C2 .
2
(
)
2
d
gt
=
−
+ C1 t + C2 = −gt + C1 ,
dt
2
(
)
d
dy(t)
=
dt
dt
)
d (
=
−gt + C1 = −g.
dt
微分方程式の検算はやさしい
得られた解は必ずもとの微分方程式に代入してみる!!
18
d2 y
微分方程式
=
−g
dt2
19
初期条件: y(0) = H, v(0) = 0.
積分定数 (C1 , C2 ) の決定: 測定デ−タに対応する関数
dy(t)
v(t) =
dt
= −gt + C1
(
=⇒ 0
C1 = 0,
gt2
y(t) = −
+ C2
2
C2 = H,
gt2
y(t) = H −
.
2
)
t=0 ,
(
=⇒ H
(1)
)
t = 0 , (2)
微分方程式の解は測定値と対応付けることができる!!
(3)
(
• 代数方程式
微分方程式
)
を解くと
)
b
y=−
が得られ，
a
ay + b = 0
(
条件をみたす数値
(
• 微分方程式
20
)
d2 y
+g =0
2
dt
(
現象を説明する関数
を解くと
)
gt2
y=H−
が得られる．
2
演習問題 – 0
Pisa の斜塔から
質量 m の物体を放した．
• 速度と高さの関係を求めよ．
• 速度と重力加速度の関係を求めよ．
• この現象を記述する微分方程式を解け．
未知関数を結びつける関係
未知関数を含む等式
を導く
を導く
21
22
1
自然現象の数学モデル (2)
冷却
23
シミュレ−ション (DE1.c)
微分方程式:
プログラム:
初期条件:
dy
= −0.5y
dt
dydt = -0.5*y;
y(0) = 1
24
(1)
(2)
(#define Y0)
方程式を定義するには (DE1.c)
#define Y0 1
....................
// Differential Equation ... start
// First Order Equation
....................
dydt = -0.5*y;
....................
// Differential Equation ... end
....................
dy
プログラミング:dydt ⇔
dt
25
冷却：Newton の法則 (1)
26
dTℓ (t)
• 温度の変化速度は熱流の強さに比例:
∝ Q(t)
dt
• 熱流の強さは温度差に比例:
y(t) = Tℓ (t) − Ta :
積分すると
Q (t) ∝ ta − Tℓ (t)
dy(t)
= −ay(t), y(0) = y0 . (1)
dt
...
∫
y(t) = −a
y(t)dt + C,
∫ t
y(s)ds.
y(t) = y0 − a
0
????
(2a)
(2b)
微分の始まり
dy(t)
dt
y=C
:
y=t
:
y = t2
:
y = tn
:
y(t + h) − y(t)
= lim
h→0
h
C −C
0
lim
= lim
= 0,
h→0
h→0 h
h
h
(t + h) − t
lim
= lim
= 1,
h→0
h→0 h
h
(t + h)2 − t2
2th + h2
lim
= lim
= 2t.
h→0
h→0
h
h
·········
dtn
= ntn−1 .
(?)
dt
27
冷却：Newton の法則 (2)
微分方程式:
dy(t)
= −ay(t),
dt
28
y(0) = y0 . (1)
冷却：Newton の法則 (3)
微分方程式:
y = tn ?:
指数関数を微分すると
29
dy(t)
= −ay(t),
y(0) = y0 . (1)
dt
dy(t)
= ntn−1 ̸= −atn .
(2)
dt
...
冷却：Newton の法則 (4)
30
指数関数の性質:
d t
e
dt
= et ,
e0 = 1.
(1)
Newton の法則
dy(t)
dt
= −ay(t),
y(0) = y0 .
(2)
冷却：Newton の法則 (5)
31
Newton の法則
dy(t)
dt
= −ay(t),
y(0) = y0 .
(1)
Newton の法則をみたす関数
一般に
y(t) = e−at ,
(2)
y(t) = Ce−at ,
(3)
Newton の法則と初期条件をみたす関数
y(t)
= y0 e−at ,
y(0) = y0 .
(4)
検算
解？:
左辺:
y(t) = e−at y0 .
dy
d ( −at )
=
e y0
dt
dt
d
−at
=e
·
(−at) · y0
dt
= −ae−at y0
(
)
= −ay(t)
= 右辺!!
32
放射性物質の崩壊
33
放射性物質の崩壊
dm(t)
dt
m(t)
= −am(t),
= Ce−at ,
34
m(0) = M.
(1)
C = m(0) = M.
(2)
35
演習問題 – 1
沸騰している水
を室内に放置した．
• 温度低下速度と室温 (あるいは温度，あるいは温度 −
室温) の関係を求めよ．
宇宙線の照射を受けた上空で放射性炭素
古生物の体内に吸収された．
14
C が生成され
• 化石中の残留 14 C 質量と時間の関係を求めよ．
（質量の変化を表す微分方程式を示せ）
• 上記微分方程式をみたす関数をひとつ示せ．
36
2
変数分離型微分方程式
微分方程式の問題
37
• ガリレオがピサの斜塔から落したボールはニュートン
の運動方程式にしたがって落下する．
• 抛物線は運動方程式をみたす．
• 飲み物の温度は徐々に室温に近付く．
• 放射性物質が崩壊する速度は徐々に 0 に近付く．
• 指数関数 e−at は緩和特性を示す．
• 現象を表す方程式をみたす「関数」を見付けられるか?
• その「関数」は観測データに適合するよう特定できるか?
変数分離型：簡単に解けるタイプ
dy(t)
dt
dy
y
∫
dy
y
log y
y(t)
(
= −ay(t),
dy
= −ay
dt
38
)
(1)
= −adt
= −at + C,
(2)
= −at + C,
= e−at eC .
(3)
初期条件： y(0) を用いて eC を定める．
変数分離型：初期条件
dy
= −ay,
dt
y(0) = y0 .
39
(1)
(2)
初期条件を定めるには. . .
y(t) = e−at eC .
y(0) = e0 eC
eC
= y0 ,
= y0 .
(3)
シミュレーション (DE1.c)
微分方程式:
プログラム:
初期条件:
dy
= −0.5y
dt
dydt = -0.5*y;
y(0) = 1
40
(1)
(2)
(#define Y0)
シミュレーション (DE1.c)
#define Y0 1
....................
// Differential Equation ... start
// First Order Equation
....................
dydt = -0.5*y;
....................
// Differential Equation ... end
....................
dy
プログラミング:dydt ⇔
dt
41
42
変数分離型：別の例
dy(t)
dt
∫
dy
y
dy
y
(
= −aty(t),
dy
= −aty
dt
)
(1)
= −atdt
t2
= −a
+ C,
2
t2
log y = −a
+ C,
2
[
]
2
at
y(t) = y0 exp −
.
2
(2)
(
y(0) = y0
)
(3)
シミュレ−ション (DE1.c)
微分方程式:
プログラム:
初期条件:
dy
= −0.5ty
dt
dydt = -0.5*t*y;
y(0) = y0
43
(1)
(2)
(#define Y0)
方程式を定義するには (DE1.c)
....................
// Differential Equation ... start
// First Order Equation
....................
dydt = -0.5*t*y;
//dydt = -0.5*y;
....................
// Differential Equation ... end
....................
dy
プログラミング:dydt ⇔
dt
44
45
一般解
(
• 微分方程式
(
したもの
dy
= −ay
dt
)
y = Ce−at
)
を満たす関数すべてを表
を一般解といい，
• t = 0 での y の値 y(0) を指定する条件を
「初期条件」という．
• 初期条件を指定することにより現象を特定できる．
初期条件を定めるには (DE1.c)
#define Y0 1
....................
// Differential Equation ... start
// First Order Equation
....................
dydt = -0.5*y;
....................
// Differential Equation ... end
....................
dy
プログラミング:dydt ⇔
dt
46
演習問題 – 2
47
放射性物質が単位時間に崩壊する量はその時点での物質の量
に比例する．
• この関係を微分方程式で表せ．
• 放射性物質の半減期を 5,000 年とするときの比例係数
を求めよ．
• 上記微分方程式を解け．
• 上記解が微分方程式をみたすことを示せ．
48
3
1 階微分方程式：未定係数法
微分方程式の問題
• 飲み物の温度は徐々に室温に近付く．
• 放射性物質が崩壊する速度は徐々に 0 に近付く．
• 指数関数 e−at は緩和特性を示す．
• 現象を表す方程式をみたす「関数」を見付ける．
（変数分離）
• 観測データに適合するよう初期値を指定する．
• 「外部入力」が加えられたらどうなる?
• 「外部入力」が振動していたらどうなる?
49
スロ−スタ−ト
50
沸騰している水を室内に放置した. . .
51
• 時刻 t ≥ 0 のときの水温を y(t), 室温を Ta (t) とする．
• 室内空気から水に流れ込む熱量は，熱抵抗 R を用いて
1
(Ta (t) − y(t)) のように表される:
R
dy(t)
• 水温の上昇速度
は流入熱量に比例するから
dt
dy(t)
C
dt
=
1
(Ta (t) − y(t)) ,
R
が得られる．この C は熱容量とよばれる．
沸騰している水を室内に放置した
微分方程式:
52
1
を用いて
a=
CR
dy
+ ay(t) = Q(t),
dt
Q(t) = aTa (t)
(1)
が得られる．
( )
(
)
すなわち Ta (t) が高低くなると水温は速くゆっくりと
上昇する．
注意:
Ta (t) は一定でなくてもかまわない!!
1 階定係数線形微分方程式
dy
+ ay
dt
= Q(t).
53
(1)
「強制項」をもつ微分方程式の解：一般解 + 特殊解
y(t) =
一般解
特殊解:
y¯(t) + η(t),
d¯
y
+ a¯
y = 0.
dt
(2a)
(2b)
Q(t) の形をみて仮定
d
η(t) + aη(t)
dt
= Q(t).
(3)
検算
解？
：
左辺:
y(t) = y¯(t) + η(t),
d¯
y
+ a¯
y = 0,
dt
dη
+ aη = Q(t).
dt
d
(¯
y + η) = −a¯
y − aη + Q(t)
dt
= −a (¯
y + η) + Q(t)
= −ay + Q(t).
(右辺)
54
未定係数法： Q(t) = A
η(t) =
dη
+ aη = ak
dt
k
=
=
A
.
a
k,
55
(1)
A.
(2)
検算
特殊解？
：
dη
+ aη
dt
A
η=
.
a
A
= a
=A
a
= Q(t).
56
シミュレ−ション (DE1.c)
微分方程式:
プログラム:
初期条件:
外力:
57
dy
= −0.5y + 1
(1)
dt
dydt = -0.5*y + 1;
(2)
y(0) = 0
(#define Y0)
A=1
(3)
温度調節
58
方程式を定義するには (DE1.c)
#define Y0 0
....................
// Differential Equation ... start
// First Order Equation
....................
//dydt = -0.5*y + sin(t);
dydt = -0.5*y + 1;
....................
// Differential Equation ... end
....................
dy
プログラミング:dydt ⇔
dt
59
未定係数法： Q(t) = sin t
η(t)
dη
+ aη
dt
= k cos t + m sin t,
60
(1)
= −k sin t + m cos t
+a(k cos t + m sin t)
= sin t,
am − k
= 1,
(2)
m + ak
= 0.
(3)
検算
特殊解？
：
dη
+ aη
dt
61
η(t) = k cos t + m sin t,
am − k = 1, m + ak = 0.
= (−k sin t + m cos t) + a(k cos t + m sin t)
= (am − k) sin t + (m + ak) cos t
= sin t
= Q(t).
シミュレ−ション (DE1.c)
微分方程式:
プログラム:
初期条件:
外力:
dy
= −0.5y + sin(t)
dt
dydt = -0.5*y + sin(t);
y(0) = 0
(#define
Q(t) = sin(t)
62
(1)
(2)
Y0)
(3)
方程式を定義するには (DE1.c)
#define Y0 0
....................
// Differential Equation ... start
// First Order Equation
....................
dydt = -0.5*y + sin(t);
//dydt = -0.5*y + 1;
....................
// Differential Equation ... end
....................
dy
プログラミング:dydt ⇔
dt
63
気温の変動
64
演習問題 – 3
(1) 微分方程式
を解け．
(2) 微分方程式
を解け．
(3) 微分方程式
を解け．
65
dy
+ ay
dt
= 0.
(1)
dy
+ ay
dt
= 1.
(2)
dy
+ ay
dt
=
sin t.
(3)
66
4
1 階微分方程式：定数変化法
67
定常解
t → ∞ のときの解:
dy
+ ay
dt
dy
=0
dt
= A,
=⇒ y∞
A
=
.
a
便利な書き方:
dy
dt
= a [y∞ − y] .
t → ∞ のとき y が一定でなければ?
(入力が変動)
68
定数変化法のアイデア
dy
+ ay
dt
= Q(t).
(1)
Q(t) = 0 の解 y¯ を求める:
(
y¯ = e−at C,
de−at
+ ae−at = 0
dt
Q(t) の影響により C が変化すると考える:
y
= e−at C(t),
C(t) がみたす微分方程式を考える: ???
)
69
dy
定数変化法
+ ay = Q(t)
dt
(
y
= e−at C(t),
de−at
+ ae−at = 0
dt
)
（左辺）⇒（右辺）
dy
+ ay
dt
=
dC(t)
e
− ae−at C(t) + ae−at C(t)
dt
= Q(t),
−at
C(t) がみたす微分方程式:
dC(t)
dt
=
eat Q(t).
(1)
dC(t)
at
不定積分
= e Q(t)
dt
∫
表記 (1)
表記 (2)
eat Q(t)dt +（積分定数）
∫ t
: C(t) = C(0) +
eas Q(s)ds
0
(
)
: C(t) =
C(0) が指定されたとき便利
C(0) = C0 と指定されたとき
∫
C(t)
t
eas Q(s)ds.
= C0 +
0
70
∫
t
71
as
不定積分
e Q(s)ds
0
表記 (2) の定義より
(∫
d
dt
(∫
...
)
t
eas Q(s)ds
= eat Q(t),
0
)
t+dt
eas Q(s)ds
≃ eat Q(t)dt
t
検算:
dC(t)
dt
(
=
)
d
C0 +
eas Q(s)ds
dt
0
(∫ t
)
d
eas Q(s)ds = eat Q(t).
=
dt
0
∫
t
72
dC(t)
at
定数変化法
= e Q(t)
dt
∫
eat Q(t)dt + C0 ,
[∫
C(t) =
y
= e−at C(t) = e−at
= e−at
∫
∫
(1a)
]
eat Q(t)dt + C0
eat Q(t)dt + e−at C0 .
(1b)
t
eas Q(s)ds + C0 ,
0
∫ t
= e−at
eas Q(s)ds + e−at C0 .
C(t) =
(2a)
y
(2b)
0
73
検算 (1)
解？
：
右辺：
y = e−at
∫
eat Q(t)dt + e−at C.
e−at F (t),
(
)
∫
∫ t
F (t) = eat Q(t)dt + C =
eas Q(s)ds + C .
0
dy
dt
]
d [ −at
=
e F (t)
dt
de−at
−at dF (t)
=
F (t) + e
,
dt
dt
74
検算 (2)
de−at
dt
dF (t)
dt
=
dy
dt
d
dt
= −ae−at ,
)
(∫
eat Q(t)dt + C
= eat Q(t),
= −ae−at F (t) + e−at eat Q(t)
[ −at
]
−a e F (t) + Q(t)
= −ay + Q(t).
定数変化法と初期値: 表記 (2)
dy
+ ay
dt
y
=
= Q(t),
y(0) = y0 .
(
e−at C(t),
∫
75
(1)
)
C(0) = y(0) = y0
t
eas Q(s)ds,
C(t) = y0 +
0
∫
y
t
eas Q(s)ds
= e−at y0 + e−at
0
∫ t
= e−at y0 +
e−a(t−s) Q(s)ds.
0
(2)
76
定数変化法の物理的意味
Q(s) Q(t)
e−a(t−s)
y(t)
t, s
室温の制御
77
演習問題 – 4
78
(1) 微分方程式
dy
+ ay
dt
= Q(t).
を解け．
(2) とくに Q(t) = A のときの解を求めよ．
(3) t = 0 で y = 0 をみたす解を求めよ．
(4) Q(t) = sin(t) のときの解はどうなるか？
(1)
79
5
微分方程式をみたす指数関数
振動を利用して情報を得る
80
振動を利用して情報を得る
81
demo ecalyi
微分方程式の問題
82
• 抛物線は運動方程式をみたす．
• 指数関数 e−at は緩和特性を示す．
• 引っ張ったバネから手を離すと振動する．
• 三角関数 (sin ωt, cos ωt) は振動波形を表す．
• 振動現象を表す方程式をみたす「三角関数」を見付け
られるか?
• その「関数」は観測データに適合するよう特定できるか?
シミュレ−ション (DE2.c)
微分方程式:
プログラム:
初期条件:
d2 y
= −9y
2
dt
d2ydt2 = -9*y;
y(0) = y0
dy(0)
= y˙ 0
dt
83
(1)
(2)
(#define Y0)
(#define DY0)
方程式の定義 (DE2.c)
....................
// Differential Equation ... start
// Second Order Equation
d2ydt2 = -9*y;
// Differential Equation ... end
....................
d2 y
dy
, d2ydt2 ⇔
プログラミング:dydt ⇔
dt
dt2
84
三角関数 A sin ωt + B cos ωt
微分方程式では
sin ωt と cos ωt はセットで扱う!!!
85
初期値を変えるには (DE2.c)
#define Y0 0//1
#define DY0 3//0
....................
// Differential Equation ... start
// Second Order Equation
d2ydt2 = -9*y;
// Differential Equation ... end
....................
dy
d2 y
プログラミング:dydt ⇔
, d2ydt2 ⇔
dt
dt2
86
87
λt
微分方程式と指数関数 y = e
(
dy
+ ay = 0,
dt
λeλt + aeλt = 0
λ2 − (ai)2 = 0
(
=
0.
(1)
y = C · e−at
λ + a = 0 =⇒
d2 y
2
+
a
y
2
dt
)
)
λ2 eλt + a2 eλt = 0
=⇒ y = C1 · eait + C2 · e−ait
(?)
演習問題 – 5
88
Hooke の法則によるとバネの伸び（縮み）に比例した力
が変位の方向と逆の向きに働く．
(i) この関係を微分方程式で表せ．
(ii) 上記微分方程式を充たす関数を一つ示せ．
指数関数 y = eλt
が微分方程式
d2 y
dy
+ 5y
+4
2
dt
dt
= 0
を充たすという．このときの λ を求めよ．
89
6
複素指数関数と微分方程式
90
iθ
Euler: e
iy
cosθ + i sin θ
i sin θ
θ
cos θ
x
複素指数関数と三角関数
91
eiθ
= cos θ + i sin θ,
(1)
e−iθ
= cos θ − i sin θ,
(2)
cos θ
=
sin θ
=
eiθ + e−iθ
,
2
eiθ − e−iθ
.
2i
(3)
(4)
複雑な振動
92
93
Euler: eiωt
iy
+ω
cosωt + i sin ωt
−ω
ωt
cos ωt
i sin ωt
x
微分方程式の問題
94
• 抛物線は運動方程式をみたす．
• 指数関数 e−at は緩和特性を示す．
• 三角関数 (sin ωt, cos ωt) は振動波形を表す．
• 指数関数と三角関数を組合せて振動現象を表す方程式
の「解＝関数」を構成できるか?
シミュレ−ション (DE2.c)
微分方程式:
プログラム:
初期条件:
d2 y
dy
= −3
− 2y,
2
dt
dt
dy
d2 y
= −2
− 2y,
2
dt
dt
d2ydt2 = -2*dydt - 2*y;
d2ydt2 = -3*dydt -2*y;
y(0) = y0
dy(0)
= y˙ 0
dt
95
(1)
(2)
(3)
(4)
(#define Y0)
(#define DY0)
方程式を定義するには (DE2.c)
....................
// Differential Equation ... start
// Second Order Equation
d2ydt2 = -2*dydt - 2*y;
d2ydt2 = -3*dydt -2*y;
//d2ydt2 = -9*y;
// Differential Equation ... end
....................
d2 y
dy
, d2ydt2 ⇔
プログラミング:dydt ⇔
dt
dt2
96
複素指数関数がみたす微分方程式 (1) 97
d λt
e
dt
d2 λt
e
2
dt
d2 iωt
e
2
dt
d2 −iωt
e
2
dt
= λeλt
(1)
= λ2 eλt
(2)
= −ω 2 eiωt
(3)
= −ω 2 e−iωt
(4)
複素指数関数がみたす微分方程式 (2) 98
⇒
d2 λt
d λt
λt
e
+
b
e
+
ce
= 0,
2
dt
dt
(λ2 + bλ + c)eλt = 0.
(λ = σ + iω)
特性方程式:
λ2 + bλ + c = (λ + s1 )(λ + s2 ),
s1 = σ + iω, s2 = σ − iω.
λ = −s1
λ = −s2
y
y = e−s1 t = e−σt e+iωt ,
y = e−s2 t = e−σt e−iωt .
= C1 e−s1 t + C2 e−s2 t .
(1)
(2)
(3)
(4)
(5)
2 階線形微分方程式（減衰解）
解の型:
y = eλt ,
99
(1)
d2
d
y+3
y + 2y
2
dt (
dt
( 2
) λt )
⇐⇒
λ + 3λ + 2 e
((
)
( 2
) λt
)
λ + 3λ + 2 e = 0
λ2 + 3λ + 2 = 0
(λ + 2)(λ + 1) = 0
y = e−2t , y = e−t .
一般解:
y
= C1 e−t + C2 e−2t .
(2)
(3)
(σ+iω)t
複素指数関数 e
[
実部:
虚部:
]
, σ < 0, ω = 0
ℜ e(σ+iω)t = eσt cos(ωt)
[
]
ℑ e(σ+iω)t = eσt sin(ωt)
実部
100
2 階線形微分方程式（振動解）
解の型:
y = eλt
(1)
d2
d
y+2
y + 2y = 0
2
dt
(dt
)
⇐⇒ λ2 + 2λ + 2 = 0
λ1 = 1 + i,
y1 = e−t−it ,
101
λ2 = 1 − i,
y2 = e−t+it .
(
)
λ 1 = λ2
(2)
一般解:
y
= C1 e−t eit + C2 e−t e−it
(3)
= (C1 + C2 )e−t cos t + i(C1 − C2 )e−t sin t
= A1 e−t cos t + A2 e−t sin t
(4)
(σ+iω)t
複素指数関数 e
[
実部:
虚部:
]
, σ < 0, ω > 0
ℜ e(σ+iω)t = eσt cos(ωt)
[
]
ℑ e(σ+iω)t = eσt sin(ωt)
実部
102
複素指数関数と三角関数
103
eiωt
= cos ωt + i sin ωt,
(1)
e−iωt
= cos ωt − i sin ωt,
(2)
cos ωt
=
sin ωt
=
eiωt + e−iωt
,
2
eiωt − e−iωt
.
2i
(3)
(4)
演習問題 – 6
つぎの微分方程式を考える．
d2 y
+ 9y
2
dt
= 0.
• 一般解を求めよ．
• 初期条件
y(π/2) = 2,
dy(π/2)
=6
dt
をみたすように係数を定めよ．
104
105
7
2 階線形微分方程式の一般解
（同次方程式）
衝撃の吸収
106
微分方程式の問題
107
• 抛物線は運動方程式をみたす．
• 指数関数 e−at は緩和特性を示す．
• 三角関数 (sin ωt, cos ωt) は振動波形を表す．
• 指数関数と三角関数を組合せて振動現象を表す方程式
の「解＝関数」を構成できるか?
• 「方程式の係数」を観測データに適合するよう特定で
きるか?
• 「解の初期値」は観測データ（「実数」!!!）に適合する
よう特定できるか?
初期値問題：単振動 (1)
d2
m 2 y
dt
=
d2
2
y
+
ω
y =
2
dt
y = eλt :
108
d
−ky,
y(0) = y0 ,
y(0) = y˙ 0 ,
dt
√
(1)
0,
ω = k/m,
( 2
) λt
2
λ + ω e = 0,
(λ − iω)(λ + iω) = 0,
y = eiωt ,
y = e−iωt .
d2
2
sin(ωt)
+
ω
sin(ωt)
2
dt
d2
2
cos(ωt)
+
ω
cos(ωt)
2
dt
(2)
=
0,
(3)
=
0.
(4)
λt
初期値問題：単振動 (2) y = e
d2
2
y
+
ω
y
2
dt
y
e±iωt
y
= 0,
((
λ2 + ω
)
2
eλt = 0
= C1 eiωt + C2 e−iωt ,
109
)
(1)
= cos(ωt) ± i sin(ωt).
(2)
= A1 cos(ωt) + A2 sin(ωt),
A1 = y0 ,
(3)
A2 ω = y˙ 0
シミュレ−ション (DE2.c)
微分方程式:
プログラム:
初期条件:
d2 y
dy
= −9y − α
2
dt
dt
d2ydt2 = -9*y;
d2ydt2 += -0*y;
y(0) = 1
dy(0)
=0
dt
110
(1)
(2)
(#define Y0)
(#define DY0)
方程式を定義するには (DE2.c)
#define Y0 1 // 0
#define DY0 0 // 1
....................
// Differential Equation ... start
// Second Order Equation
d2ydt2 = -9*y;
d2ydt2 += -0*dydt;
// Differential Equation ... end
....................
d2 y
プログラミング:d2ydt2 ⇔
dt2
111
初期値問題：減衰振動 (1)
d2 y
m 2
dt
dy
= −ky − f
,
√ dt
k
ωn =
,
m
y = eλt :
112
dy(0)
y(0) = y0 ,
= y˙ 0 ,
dt
1
f
ζ=
.
2ωn m
d2 y
dy
2
+
2ζω
+
ω
n
n y = 0.
2
dt
dt
λ2 + 2ζωn λ + ωn2 = 0
(λ + s1 ) (λ + s2 ) y = 0,
λ = s1 , λ = s2 ,
y = C1 e−s1 t , y = C2 e−s2 t .
(1)
(2)
初期値問題：減衰振動 (2)
d2 y
m 2
dt
一般解:
dy
= −ky − f
,
dt
113
dy(0)
y(0) = y0 ,
= y˙ 0 .
dt
y = C1 e−s1 t + C2 e−s2 t ,
(1)
λ2 + 2ζωn λ + ωn2 = (λ + s1 )(λ + s2 ) = 0,
√
si /ωn = −ζ ± ζ 2 − 1,
(2)
√
k
1
f
ωn =
, ζ=
.
m
2ωn m
初期条件:
C1 + C2 = y0 ,
(3)
−s1 C1 − s2 C2 = y˙ 0 .
(4)
シミュレ−ション (DE2.c)
微分方程式:
プログラム:
初期条件:
d2 y
dy
= −2
− 2y,
2
dt
dt
d2ydt2 = -2*dydt - 2*y;
y(0) = 1
dy(0)
=0
dt
114
(1)
(2)
(#define Y0)
(#define DY0)
方程式を定義するには (DE2.c)
#define Y0 1
#define DY0 0
....................
// Differential Equation ... start
// Second Order Equation
d2ydt2 = -2*dydt - 2*y;
// Differential Equation ... end
....................
dy
d2 y
プログラミング:dydt ⇔
, d2ydt2 ⇔
dt
dt2
115
演習問題 – 7
116
パラメ−タ ζ, (ζ > 0), ωn , (ωn > 0), が変化したときの微分
方程式
d2 y
dy
2
+
2ζω
+
ω
n
n y = 0,
2
dt
dt
dy(0)
= 0,
y(0) = 1,
dt
の解の様子を示せ．
117
8
未定係 y 数法
（非同次方程式）
衝撃の吸収
118
乗り心地
119
微分方程式の問題
120
• 指数関数 e−at は緩和特性を示す．
• 三角関数 (sin ωt, cos ωt) は振動波形を表す．
• 指数関数と三角関数を組合せて振動現象を表す方程式
の「解＝関数」を構成できる．
• 「外部入力」が加えられたらどうなる?
• 「外部入力」が変化したらどうなる?
衝撃の吸収 (1)
d2 y
dy
2
+
ω
+
2ζω
n
ny
2
dt
dt
= 0,
y(0) = A,
121
dy(0)
= 0.
dt
衝撃の吸収 (2)
d2 y
dy
2
+
ω
+
2ζω
n
ny
2
dt
dt
= A
122
123
定常入力 (1)
d2 y
dy
2
+
2ζω
+
ω
n
ny = A
2
dt
dt
同次方程式：s2 + 2ζωn s + ωn2 = (s + s1 )(s + s2 ) = 0
d2 y1
+ 2ζωn
2
dt
d2 y2
+ 2ζωn
2
dt
一般解：(y1 , y2 )
dy1
+ ωn2 y1
dt
dy2
+ ωn2 y2
dt
(
y
= C1 y1 + C2 y2 .
y1 = e−s1 t ,
(1)
=
0,
(2)
=
0.
(3)
y2 = e−s2 t
)
124
定常入力 (2)
入力： Q(t) = A
y
= C1 y1 + C2 y2 + η.
(
)
y1 = e−s1 t , y2 = e−s2 t
(1)
特殊解： η = const.
d2 η
dη
2
+
ω
+
2ζω
n
nη
2
dt
dt
ωn2 η
= A.
= A.
(2)
シミュレ−ション (DE2.c)
微分方程式:
プログラム:
初期条件:
外力:
d2 y
dy
= −2
− 2y + 3
2
dt
dt
d2ydt2 = -2*dydt -2*y + 3;
y(0) = 0
dy(0)
=0
dt
A=3
125
(1)
(2)
(#define Y0)
(#define DY0)
(3)
方程式を定義するには (DE2.c)
#define Y0 0
#define DY0 0
....................
// Differential Equation ... start
// Second Order Equation
d2ydt2 = -3*dydt -2*y + 3;
d2ydt2 = -2*dydt -2*y + 3;
// Differential Equation ... end
....................
dy
d2 y
プログラミング:dydt ⇔
, d2ydt2 ⇔
dt
dt2
126
127
乗り心地
d2 y
dy
2
+
2ζω
+
ω
n
ny
2
dt
dt
= A sin αt
未定係数法： Q(t) = A sin αt
η(t)
=
= k cos αt + m sin αt,
128
(1)
d2 η
dη
2
+
2ζω
η
+
ω
n
nη
2
dt
dt
−kα2 cos αt − mα2 sin αt
+2ζωn (mα cos αt − kα sin αt)
+ωn2 (k cos αt + m sin αt) = A sin αt.
sin αt :
cos αt :
−mα2 − 2kζωn α + mωn2 = A,
−kα2 + 2mζωn α + kωn2 = 0.
−2ζωn α · k + (ωn2 − α2 )m = A,
(ωn2 − α2 )k + 2ζωn α · m = 0.
(2)
(3)
シミュレ−ション (DE2.c)
129
プログラム:
d2 y
dy
= −3
− 2y + 3 sin(2t)
(1)
2
dt
dt
d2ydt2 = -3*dydt -2*y + 3*sin(2*t);
初期条件:
(2)
(#define Y0)
微分方程式:
外力:
y(0) = y0
dy(0)
= y˙ 0
dt
Q(x) = 3 sin(2t)
(#define DY0)
(3)
方程式を定義するには (DE2.c)
....................
// Differential Equation ... start
// Second Order Equation
d2ydt2 = -3*dydt -2*y + 3*sin(2*t);
// Differential Equation ... end
....................
dy
d2 y
プログラミング:dydt ⇔
, d2ydt2 ⇔
dt
dt2
130
演習問題 – 8
つぎの微分方程式の一般解を求めよ．
d2 y
dt2
= −y + 3 sin 2t
131
132
9
連立微分方程式
133
微分方程式の問題
• 指数関数 e−at は緩和特性を示す．
• 三角関数 (sin ωt, cos ωt) は振動波形を表す．
• 指数関数は複素数に拡張できる
e−at =⇒ e(σ+iω)t .
• 2 つの微分方程式
dy
+ ay
dt
dy
d2 y
2
+
2ζω
+
ω
n
ny
2
dt
dt
は別物か?
= 0,
= 0,
連立微分方程式への書き換え (1)
d2 y
dy
+a
+ by = Q(t).
2
dt
dt
dy
x1 = y,
x2 =
,
dt
dx1
dy
=
= x2 ,
dt
dt
dx2
d2 y
dy
=
= −a
− by + Q(t)
2
dt
dt
dt
= −bx1 − ax2 + Q(t),


dx1
]
[
x2
 dt 
=
 dx 
−bx1 − ax2 + Q(t)
2
[
] [
]
dt
x2
0
=
+
.
−bx1 − ax2
Q(t)
134
(1)
連立微分方程式への書き換え (2)
d2 y
dy
+a
+ by
2
dt
dt
= Q(t).
135
(1)
[ ]
[
d x1
0
=
−b
dt x2
][ ] [
]
1
x1
0
+
,
−a x2
Q(t)
dy
x1 = y,
x2 =
.
dt [ ]
]
[
]
[
0
1
x1
0
,x =
,B =
,
A=
−b −a
x2
Q(t)
dx
dt
= Ax + B.
(2)
シミュレ−ション (DE.c)
微分方程式:
プログラム:
初期条件:
外力:
dx
= y,
dt
dy
= −2x − 2y + 3
dt
dxdt = y;
dydt = -2*x -2*y + 3;
x(0) = 0
y(0) = 0
A=3
136
(1)
(2)
(3)
(4)
(#define X0)
(#define Y0)
(5)
方程式を定義するには (DE.c)
#define Y0 0
#define X0 0
....................
// Differential Equation ... start
// System
dxdt = y;
dydt = -2*x -2*y + 3;
// Differential Equation ... end
....................
dy
プログラミング:dydt ⇔
dt
137
対角マトリクスをもつシステム (1) 138
dx
dt
=
Ax + B,
マトリクスの指数関数
[ −λ t
]
1
e
0
eAt =
−λ2 t ,
0
e
[
−λ1
A=
0
0
e = I,
[
]
0
.
−λ2
]
At −1
e
(1)
= e−At .
指数関数を解にもつ微分方程式群
dx1
+ λ1 x1 = 0,
dt
dx2
+ λ2 x2 = 0.
dt
(2)
対角マトリクスをもつシステム (2) 139
微分方程式のシステム
[ ]
x1
x=
,
x2


dx1
dx
 dt 
=  dx 
2
dt
dt
[ ]
C1
x0 =
,
C2
[
]
−λ1 x1
=
−λ2 x2
[
−λ1
=
0
= Ax,
0
−λ2
]
[ ]
x1
,
×
x2
x(0) = x0 .
対角マトリクスをもつシステム (3) 140
微分方程式群
dx1
+ λ1 x1 = 0,
dt
初期値 C1 , C2 をもつ解
微分方程式のシステム
dx2
+ λ2 x2 = 0.
dt
x1 = e−λ1 t C1 ,
dx
= Ax
dt
初期値 x0 をもつ解
[ ] [ −λ t ] [ −λ t
x1
e 1 C1
e 1
x =
= −λ2 t
=
x2
e
C2
0
=
eAt x0 .
x2 = e−λ2 t C2 .
0
e−λ2 t
]
[ ]
C1
×
C2
対角マトリクスをもつシステム (4) 141
外力を受けるシステム
定数変化法
x(t) =
dx
=
dt
dC(t)
dt
=
x(t) =
dx
= Ax + B
dt
eAt C(t),
)
d ( At
At
At dC(t)
e C(t) = Ae C(t) + e
dt
dt
At dC(t)
= Ax + e
= Ax + B,
dt
e−At B.
∫
t
eAt x0 +
eA(t−s) B(s)ds.
0
非対角マトリクスをもつシステム
マトリクス A の対角化
つぎの等式をみたす T をみつける．
[
][
k
k12 a11
T −1 AT = 11
k21 k22 a21
[
]
−λ1
0
=
= Λ.
0
−λ2
システムの対角化
dx
dt
d˜
x
dt
][
a12 k11
a22 k21
k12
k22
]−1
A = T ΛT −1
= Ax + B = T ΛT −1 x + B,
˜
= Λ˜
x + B,
142
x
˜ = T −1 x,
˜ = T −1 B.
B
143
dx
システム
= Ax
dt
[
−λ1
A=
0
0
−λ2
]
At
=⇒ e
[ −λ t
e 1
=
0
y = C1 e−λ1 t + C2 e−λ2 t
[
]
[ −(σ+iω)t
e
−σ −ω
A=
=⇒ eAt =
ω −σ
0
0
e−λ2 t
]
.
一般解
一般解
0
−(σ−iω)t
e
]
.
y = (A sin ωt + B cos ωt) e−σt
指数関数と三角関数は
連立微分方程式を介して結び付く!
マトリクスの指数関数
144
微分方程式のシステム
dx
dt
解ベクトル:
[
−λ1
A=
0
[
−σ
A=
ω
= Ax,
[ ]
y0
x(0) = x0 =
.
y˙ 0
(1)
x = eAt x0 , e0 = I
]
]
[ −λ t
0
e 1
0
:
exp [At] =
−λ2 t ,
−λ2
0
e
]
[ −(σ+iω)t
]
−ω
e
0
:
exp [At] =
−σ
0
e−(σ−iω)t
]
[
−σt cos(ωt) − sin(ωt)
.
=e
sin(ωt) cos(ωt)
演習問題 – 9
つぎの微分方程式を 1 階微分方程式に変換せよ．
d2 y
dt2
= −y + 3 sin2 t
145
146
10
微分方程式とシステム
147
電離層モデル
Solar Wind
Cosmic Ray
Solar Radiation
He++
n+ e−
5780K
450Km/s
1.4KW/m2
− e−
e−
e− −
e− − e −
−
e
e
e
e
Kennelly-Heaviside Layer
300KV
+
n+
n-
+
Electrosphere
n+
+
Thunder Cloud
-
-
Rn
+
Earth
2 つの微分方程式
de−
dt
dn+
dt
148
= ae− ,
(1)
= −bn+ .
(2)
爆発? 消滅?
de−
= ae− ,
dt
dn+
= −bn+ .
dt
149
微分方程式のシステム
(Lotka-Volterra)
de−
dt
dn+
dt
150
= α(n+ )e− ,
(1)
= −β(e− )n+ .
(2)
天敵の役割
151
demo lv/ LV P GRAPH, LV L GRAPH
演習問題 – 10
微分方程式
dx
dt
dy
dt
= a (1 − y) x,
= −b (1 − x) y,
の解 (x, y) が閉じた曲線上を移動することを示せ．
152
応用:視野外道路パタ−ンの追跡
153
demo im