常微分方程式を学ぶ— 基礎の学習と数値シミュレーション —桂田 祐史

2022 年 3 月 9 日, 2024 年 2 月 17 日

https://m-katsurada.sakura.ne.jp/labo/text/ode-workbook.pdf
https://m-katsurada.sakura.ne.jp/labo/text/ode-workbook/(PDF を正式版とする。HTML 版はあちこちで組版がおかしくなっているが、プログラム例をコピペするときは HTML 版の方が便利なので、両方公開する次第。)

目次
1

はじめに

4

マルサス

2

Malthus の法則
2.1 Malthus のモデルの微分方程式 . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 一般解 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3 方程式の表す現象を考える . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4 初期値問題の解 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.5 数値解法 — Euler 法と Runge-Kutta 法 . . . . . . . . . . . . . . . . . . . . . .
2.5.1 Euler 法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.5.2 Runge-Kutta 法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4
4
5
5
7
7
8
12

3

logistic 方程式
3.1 Malthus のモデルの修正としての logistic 方程式 . . . . . . . . . . . . . . . . .
3.2 初期値問題の解 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3 解の漸近挙動, 平衡点とその安定性 . . . . . . . . . . . . . . . . . . . . . . . .

15
15
16
16

4

単振動の方程式
4.1 単振動の方程式とは . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 単振動の方程式の解 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3 1 階正規形方程式への帰着 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.4 力学系とは . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.5 Python, Julia で数値計算 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

19
19
19
20
20
21

1

5

Lotka-Volterra の方程式 (捕食者被食者の数理モデル)
5.1 Lotka-Volterra 方程式の勉強のススメ . . . . . . . . . . . . . . . . . . . . . . .
5.2 Lotka-Volterra 方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.3 まずは数値シミュレーション . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.4 解を求めて … まず局所解の存在と一意性 . . . . . . . . . . . . . . . . . . . .
5.5 平衡点とその安定性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.6 解軌道の方程式, 解が R 全体に延長できることと中立安定性の証明 . . . . . .
5.7 解の周期性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.8 解の平均値 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.9 1 つのルーツの紹介: D’Ancona の疑問と Volterra の回答 . . . . . . . . . . . .
5.10 Lotka-Volterra 方程式に続いて — 生物の競争系のモデル . . . . . . . . . . . .
5.10.1 競争系とは . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.10.2 2 種競争系 — Lotka-Volterra competition model . . . . . . . . . . . . .
5.10.3 3 種競争系 — May-Leonard のモデル . . . . . . . . . . . . . . . . . . .
5.10.4 魚種交替 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

22
22
22
23
23
24
25
27
28
28
29
29
30
30
31

6

SIR モデル
6.1 感染症の数理モデルとしての SIR モデル . . . . . . . . . . . . . . . . . . . . .
6.2 数値シミュレーション . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.3 数学的解析 (1) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.4 数学的解析 (2) 無次元化と基本再生産数 . . . . . . . . . . . . . . . . . . . . .

32
32
34
34
37

7

挑戦課題
7.1 三種競争系 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.2 渦糸の力学系 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.3 メトロノームの同期現象 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

40
40
40
40

8

参考書案内

40

A 問題の解答42

B 常微分方程式の初期値問題の数値シミュレーションの手引き
B.1 はじめに . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
B.2 Python によるシミュレーション . . . . . . . . . . . . . . . . . . . . . . . . .
B.2.1 Malthus モデル . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
B.2.2 単振動の方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
B.3 Julia によるシミュレーション . . . . . . . . . . . . . . . . . . . . . . . . . . .
B.3.1 Julia のインストール＆使い方 . . . . . . . . . . . . . . . . . . . . . . .
B.3.2 Malthus モデル . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
B.3.3 単振動の方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

46
46
46
47
52
57
57
57
59

C 定数係数線形同次常微分方程式に対する特性根の方法 (素朴な説明)
C.1 とりあえず要点を学ぶ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
C.2 基本的な定理の証明 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
C.2.1 定義と例 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
C.2.2 特性方程式, 特性根 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
C.2.3 相異なる特性根を持つ場合 . . . . . . . . . . . . . . . . . . . . . . . . .
C.2.4 特性根が重根である場合 . . . . . . . . . . . . . . . . . . . . . . . . . .

62
62
64
64
64
65
66

2

C.2.5 特性根が虚数である場合 . . . . . . . . . . . . . . . . . . . . . . . . . .
C.3 まとめ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6872

D 初期値問題の基礎理論 (かけ足で説明)
D.1 はじめに . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
D.2 解の存在 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
D.3 解の存在範囲 (爆発, 極大延長解) . . . . . . . . . . . . . . . . . . . . . . . . .
D.4 解の一意性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
D.5 まとめ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

72
72
73
74
76
78

E 初期値問題の解の存在と一意性の証明
E.1 Lipschitz 条件を仮定する場合 . . . . . . . . . . . . . . . . . . . . . . . . . . .
E.1.1 Banach 空間 C([a, b]; R) . . . . . . . . . . . . . . . . . . . . . . . . . .
E.1.2 不動点, 縮小写像に関する不動点定理 . . . . . . . . . . . . . . . . . . .
E.1.3 連続かつ Lipschitz 条件を満たす場合の解の存在と一意性 . . . . . . . .
E.2 Lipshitz 条件を仮定しない場合 . . . . . . . . . . . . . . . . . . . . . . . . . .

78
78
79
79
81
82

F 線形微分方程式の常識
F.1 取り扱う方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
F.2 高階微分方程式の 1 階微分方程式への帰着 . . . . . . . . . . . . . . . . . . . .
F.3 解の存在と一意性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
F.4 同次方程式の解空間は線形空間 . . . . . . . . . . . . . . . . . . . . . . . . . .
F.5 線形同次方程式の解の 1 次独立性, Wronskian による判定 . . . . . . . . . . . .
F.5.1 ロンスキアンを用いない方法 . . . . . . . . . . . . . . . . . . . . . . .
F.6 線形同次方程式の解空間の次元 . . . . . . . . . . . . . . . . . . . . . . . . . .
F.7 定数係数高階単独線形常微分方程式 . . . . . . . . . . . . . . . . . . . . . . . .
F.8 線形非同次方程式の解空間の構造 . . . . . . . . . . . . . . . . . . . . . . . . .

82
82
83
83
83
84
87
90
90
93

G 力学系についてのメモ
G.1 力学系とは . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
G.2 平衡点とは . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
G.3 平衡点の安定性、漸近安定性 . . . . . . . . . . . . . . . . . . . . . . . . . . .
G.4 線形安定性解析 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
G.5 復習: 定数係数線形常微分方程式 . . . . . . . . . . . . . . . . . . . . . . . . .

93
93
94
94
95
96

H 本文中の図の描き方
97H.1 図 1 の描き方 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
H.2 図 2 の描き方 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
H.3 matplotlib で日本語タイトルを用いる . . . . . . . . . . . . . . . . . . . . . . . 101
I

微分多項式についての注意
103
I.1 個人的な見解・意見 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
I.2 高橋 [1] の記述を検討する . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104(しばらく工事中の予定)書き換えをする。方針は、定義と定理は本文中に書く。

3

はじめに

1

私自身は偏微分方程式の数値計算法の数理にもっとも興味があり、「それが専門です」とい1
2
うことにしているけれど 、常微分方程式 を学ぶことは重要であり、学部ゼミの卒業研究でも現象のモデルとして良く現れる3 。
現象数理学科の 3 年生は、常微分方程式について必要なことの多く (解ける方程式の解き方、数値シミュレーションの方法など) は、色々な授業で習ってきたと思うが、諸君の先輩達を見る限り、実際に使いこなすにはトレーニングが必要のようだ。そういうトレーニングに便利なように、この文書を用意した。時間の限られている現象数理研究 I で利用することを考えている。
多くのテキストでは、何か理論や方法を説明して、それの適用例として具体的な微分方程式を持ち出して来るが、ここではまず具体的な微分方程式を取り上げ、それについてどういう理論・方法が使えるか、という順に話を進める。
勉強する場合は、以前勉強したときのノートや教科書等の資料を手元において、調べながら読み進めると良い。
色々な用語や定理、コンピューターの利用法が現れる。そのうちのいくつかは知っているか、目にした覚えがあることを期待しているが、知らないものがあっても、心配する必要はな
い。とりあえず頭に入れておいて、将来本格的に使うことになった場合は、勉強しなおそうく
らいのノリで構わない。

マ

ル

サ

ス

2

Malthus の法則

2.1

Malthus のモデルの微分方程式

実変数 t についての未知関数 x(t) について

(2.1)

x′ (t) = ax(t)という微分方程式を考える。ここで a は実定数であるとする。
この方程式は次のようにも書かれる。

(2.2)

dx
= ax.
dt
マ ル サ ス

a > 0 の場合として、人口論の Malthus の法則のモデルがある。
a < 0 の場合として、放射性元素の崩壊のモデルがある。
用語の確認をしておこう。

• 微分方程式に未知関数の n 階導関数が現れるが、n + 1 以上の階数の導関数が現れない
とき、その微分方程式は n 階の微分方程式であるという。
• n 階の微分方程式が、n 階導関数について解かれた形 (x(n) = f (t, x, x′ , · · · , x(n−1) )) をしているとき、その微分方程式は正規形 (normal form) であるという。
1

脱線かもしれないが、大学院で指導するときは「研究テーマはなるべく偏微分方程式の数値計算法にして下さい」ということにしているが、学部卒研では研究テーマは自由としている。
2微分方程式は、独立変数の個数が 1 であるとき常微分方程式 (ordinary diﬀerential equation)、2 以上である
とき偏微分方程式 (partial diﬀerential equation) と呼ばれる。
3
常微分方程式自体がある現象のモデルとなることも多いし、偏微分方程式を扱うために常微分方程式が必要になることも多い。

4

方程式 (2.1) は、1 階正規形の微分方程式 (normal form of ﬁrst order diﬀerential equation)である。
一般に、適当な関数 f = f (x, t) を用いて

x′ (t) = f (x(t), t)(2.3)の形に表される微分方程式が 1 階正規形の微分方程式である。微分方程式についての定理や、数値計算法の多くがこの形の方程式について説明される。
問(2.1) の場合に f は何か。

2.2

一般解

この方程式 (2.1) の解を求めるのは簡単である。

(2.4)

x(t) = Ceat

(C は任意定数)が方程式 (2.1) の一般解である。
問 1. (2.4) が (2.1) の一般解というのは、どういうことか述べよ。 (注: 問の解答は付録 A 節(p. 42〜) にある。)
必修課題 2.2.1 (2.1) は変数分離型の微分方程式であり、有名な解法がある。それに従い (2.4)
を導け。
必修課題 2.2.2 (2.1) は定数係数線形常微分方程式であり、定数係数線形常微分方程式に対しては特性根の方法という有名な解法がある (付録の C 節に説明してある)。特性根の方法を説明し (定理を紹介せよ)、それに従い (2.4) を導け。

2.3

方程式の表す現象を考える(この項は雑談みたいなものです。)
微分方程式についての数学のテキストでは、どういう形の微分方程式がどういう方法で解けて、どういう性質を持つか、数学的に説明する (当然のことである)。
一方、多くの研究で微分方程式が取り扱われる理由は、何かの問題を解決したかったからであることが多い （これも当然である)。数学のテキストを作る場合、一般的に役立つ内容を書こうとするので、この部分の説明が省略されたり、ごく簡単に済まされることが多い。ときには、数学的な準備が十分に整わないうちに問題に挑戦することも珍しくなく、反対にその問題を解決するために数学的なことを発展させる、ということもある。
何かの現象について、微分方程式でモデルを作り、それを解析することで解答を得た、という論文を読んだときに、微分方程式と、その解析結果だけを読んで、論文の問題に意識を向けない学生が時々いる。そういう学生と初めて出会ったとき、非常に驚いたのだが、それが唯一の例というわけでなく、程度の差はあるものの似たようなケースにしばしば出食わすようになった。
指導するときにそういうことを念頭におく必要があると考えるようになった。どうか、微分方程式を学ぶ際に、なぜそれを調べる必要があったかにも目を向けるようになって欲しい。
というわけで次の問題を出す (慣れていなくても、30 分程度で出来るかな？)。マルサスの話は非常に有名なので、資料を見つけるのは難しくないであろう。大事なのはきちんと読み取れるか。

5

調べ物課題 2.3.1 マルサスとはどういう人か (生年、没年、出身地、主な業績など)。マルサスの法則とは何か。マルサスは、著書の中でどのように説明しているか。
余談 2.1 この微分方程式は、解が指数関数になるが、指数関数が解となる微分方程式の典型的な場合であると言える。この文章を書いている 2022 年、COVID-19 (新型コロナウイルス
感染症) の流行が収まらないが、しばしば「指数関数を理解しないと」という言葉が出て来る。
マルサスは人口と食糧のバランスについて警告したが、感染者数と医療サービスのバランスに置き換えてみると、200 年以上経っても指数関数は社会常識になっていないのか、とマルサスに嘆かれるだろうか。
余談 2.2 マルサスの警告を考えると、1972 年に発表された「成長の限界」(by ローマ・クラ
ブ, “The Limits to Growth”) も想起される。
Wikipedia から
The Limits to Growth (LTG) is a 1972 report[2] on the exponential economic
and population growth with a ﬁnite supply of resources, studied by computer
simulation.



日本のウィキペディアには「成長の限界は、化学肥料が安定供給されたことにより克服された」という意味のことが書いてあるが、すんなりと納得して良いものか。
指数関数についての算数的な問題 7 日間 (一週間) で 2 倍に増えるものがあるとする。最初に10 あったとして、何日経つと 10 万 (1 億の 0.1%) 以上になるか求めよ。
答えは略。意外に時間がかかる (どう感じるかは人によるだろうけれど)。感染症は流行初期に感染者数が指数関数的に増加することが知られている。日本 (人口 1.2 億) に 10 人の感染者が侵入してきて、1 週間ごとに感染者が 2 倍になるとしたとき、例えば人口の 0.1% になるまで、3 ヶ月以上かかる、ということである。指数関数が非常に速いスピードで増加することを知っている人は多いであろうが、別の見方をすると、気がつくくらいまで増えるのに (1000
人に 1 人のレベルだと、身近な範囲では感染した人に出会えるかどうか、報道されればわかるだろうけれど)、かなり時間がかかるということである。一方、気づいたときには、対処が難しい (10 万人をすぐにきちんと隔離するのは無理、1 週間ぐずぐずしていると新しい 20 万人の患者がやってくる)。
マルサスの法則は、増加する場合の話であるが、一方で同じ形の方程式が何かが減少する現象のモデルとして使えたりする。有名なのは放射性元素の崩壊であろう。
調べ物課題 2.3.2 放射性元素の崩壊のモデルとするとき、a は何か説明せよ。具体的な場合に数値を紹介せよ (いわゆる半減期から計算できる)。

Cf. 時定数 a > 0 の場合、t が 1/a 増加すると x は e 倍になる。a < 0 の場合、t が 1/|a|
じていすう

増加すると x は 1/e 倍になる。このように値が e 倍または 1/e 倍になる量を時定数4 (timeconstant) と呼ぶ。世間の多くの人にとって (世間の人でなくても？)、e はなじみが薄い (あるいは考えにくい) ので、1/2 倍になる半減期 (half-life) や 2 倍になる倍化時間 (doubling time)が持ち出される。

4

時定数は「ときていすう」とも読むそうだ。

6

2.4

初期値問題の解

一般に微分方程式の解は無数に存在するが、何か条件を加えると解が一つに確定することが多い。代表的なものに (初期条件を加えた) 初期値問題がある。
次の (2.5a), (2.5b) を満たす x = x(t) を求めよ、という問題は初期値問題である。(2.5b) が初期条件である。x0 は初期値, t0 は初期時刻と呼ばれる。

(2.5a)x′ (t) = ax(t),(2.5b)x(t0 ) = x0 .

この初期値問題の解はx(t) = x0 ea(t−t0 )(2.6)

である (「解は…である」で、「解がただ一つ存在する」ということを同時に主張している)。
一般に微分方程式の解のグラフのことを、解曲線 (solution curve) または積分曲線 (integral
curve) と呼ぶ5 。後で出て来る解軌道と混同しないように注意が必要である。
0.
0.2 t4

0.4 t3

0.6 t0.8 t2

1. t
1.2 t1

1.4 t
0.2

0.4

0.6

0.8

1.0

t

1.6 t

図 1: 色々な初期値に対する解曲線 (解のグラフ)

必修課題 2.4.1 図 1 のような解のグラフを描け (方法は何でも OK)。
付録 H.1 にこの図をどう描いたか説明してある。
問 2. (2.6) は初期値問題 (2.5a), (2.5b) の唯一の解であるが、それはどうして分かるか (答えは一通りではなく、色々なやり方があるが、重要なことはどういうやり方にしても何か 1 つ答を持っていることである)。

2.5

数値解法 — Euler 法と Runge-Kutta 法

ここでは 1 階正規系の常微分方程式の初期値問題(2.7)

x′ (t) = f (x(t), t),(2.8)

x(t0 ) = x0

しばしば「積分する」という言葉は、微分方程式を解く、という意味に使われる。多分「解」≒ 「積分」というニュアンスだろう (と個人的に考えている)。
5

7

に対する近似解法である Euler 法と、Runge-Kutta 法を説明する。
このうち実用性があるのは Runge-Kutta 法である。Euler 法は効率が良くないが (普通にやると精度が低く、精度を上げようとすると必要な計算量が爆上がりになる) 、原理が理解しやすいため説明されることが多い。一方、Runge-Kutta 法は原理が分かりにくい。ここでは、Euler 法を効率が良くなるように改良したものである、ということを理解してもらえれば良い。
余談 2.3 Runge-Kutta 法はコーディングがほどほどに簡単な数値解法のうちでは、多くの目的6 に十分な性能を持っているので、従来はお勧めアルゴリズムの定番であった。最近の数値計算環境では、高品質なアルゴリズムが簡単に利用可能になっているため、Runge-Kutta 法
推しという相場も徐々に変わっていくのかもしれない、と私は考えるようになった。

t について、有限の範囲 [t0 , T ] で解くことにする (t0 , T は適当に定めた実数である)。
以下に提示するサンプル・プログラム (euler.c, euler.py, runge-kutta.c) は、Malthusの方程式で a = 1, t0 = 0, x0 = 1 の場合を、t について [t0 , T ] = [0, 1] の範囲で解くことにしている。つまり(2.9)

x′ (t) = x(t) (0 < t < 1),(2.10)

x(0) = 1.

この初期値問題の解は x(t) = et である。
注: 数値計算をする場合に得られるのは近似解であるが、そのときもとの微分方程式の解のことを (誤解のないように) 厳密解 (exact solution) と呼んで区別することがある7 。
これは簡単で特殊な問題であるが、ここで説明する計算手順自体はより一般の問題にそのまま適用できる。
常微分方程式の数値解法の理論的なことは、まずは桂田 [2] を見よ。

2.5.1

Euler 法[t0 , T ] に属する t について解 x(t) を “知る” ため、
t0 < t 1 < t 2 < · · · < t N = Tとなる tn (n = 1, 2, · · · , N ) を取り、x(tn ) を近似的に求める (離散変数法という) ことを目標
とする。
そのため

(2.11)

xn+1 = xn + ∆tn f (xn , tn ),∆tn := tn+1 − tn

により {xn }0≤n≤N を定め、xn を x(tn ) の近似値として採用するのが、(前進) Euler 法 (the(forward) Euler method) である。これはx(tn + ∆tn ) − x(tn )∆tn →0∆tn

x′ (tn ) = limに基づく前進差分近似(2.12)x′ (tn ) ≒x(tn+1 ) − x(tn )∆tn6いわゆる硬い方程式など、古典的 Runge-Kutta 法では手に余る問題もあり、それはそれで重要であるが。
私は学生のとき、厳密解に対して「真の解」という表現を使ったことがあるが、「言わんとすることは分かる
けれど、それは普通に使われる表現ではない。厳密解と言いなさい。」と指導された記憶がある。
78

を用いたものである。つまり、微分方程式 x′ (t) = f (x(t), t) を t = tn で考え、左辺の微分係数を (2.12) の右辺で置き換え、j = n, n + 1 について x(tj ) を xj と書き換えると、(2.11) が導出できる。

C 言語プログラマー向けの説明「常微分方程式の初期値問題を解くプログラムの書き方§ 2 Euler 法, Runge-Kutta 法入門— 1 次元の問題」8 特に「2.3.1 Euler 法の C プログラム例」9 にもう少し詳しい説明がある。
そこから C 言語で書かれたサンプル・プログラムを抜き出す。

euler.c/** euler.c (Euler method for Malthusian model)*/#include <stdio.h>
double a = 1.0;
int main(void){
int i, N;
double t, x, dt;
double f(double, double), x0;
double Tmax;// 初期値x0 = 1.0;// 最終時刻
Tmax = 1.0;// 時間刻み
printf("# N: "); scanf("%d", &N);
dt = Tmax / N;// 初期値
t = 0.0;x = x0;
printf("# t x\n");printf("%f %f\n", t, x);// Euler 法
for (i = 0; i < N; i++) {x += dt * f(x, t);t = (i + 1) * dt;
printf("%f %f\n", t, x);}
return 0;}
double f(double x, double t){return a * x;}コンパイルして実行すると N の値をどうするか尋ねてくるので、10, 100 のような “大きい” 数を入力すると、tj と xj (j = 1, 2, · · · , N ) を表示する。

89https://m-katsurada.sakura.ne.jp/labo/text/intro-ode-simulation/node7.html
https://m-katsurada.sakura.ne.jp/labo/text/intro-ode-simulation/node14.html

9

コンパイル&実行% cc euler.c% ./a.out# N: 100# t x0.000000 1.000000
中略
1.000000 2.704814
%







初期値問題の解は x(t) = e であるから、最後の数は x(1) = e = 2.7182818284590 · · · に近いはずである。N を大きくすると、e に近くなることが分かるはずである (試してみよう)。
t

データファイル作成&解曲線描画

% ./a.out > euler.data
100% gnuplot
ここで gnuplot の起動メッセージが表示される (省略する)。

gnuplot> plot "euler.data" with lp
これで解のグラフ (解曲線) が表示される。以下は、グラフを画像ファイルとして保存する
手順を示す。

gnuplot> set term png
gnuplot> set output "euler.png"
gnuplot> replot
gnuplot> quit%
これで画像ファイル euler.png ができる。


必修課題 2.5.1 上のプログラムを実行して、解のグラフ (解曲線) を描いてみよ。

Python プログラマー向けの説明 上の euler.c を Python で書き換えると、以下のようになる。

10

euler.py# euler.py# C プログラムの Python バージョン# euler.c(http://nalab.mind.meiji.ac.jp/~mk/labo/text/ode-workbook/node10.html)# Euler 法という訳で実用性はない。参考プログラム。
def f(x, t):return a * x
a=1
x0=1
Tmax=1
str=input(’N=’); N=int(str);
dt=Tmax/N
t=0.0; x=x0
print(’%f %f’ % (t,x))
for i in range(N):x += dt * f(x, t)t = (i+1) * dt
print(’%f %f’ % (t,x))

グラフを描くように書き直すのは簡単である。matplotlib を使ってみよう。
euler2.py# euler2.py
import numpy as np
import matplotlib.pyplot as plt
def f(x, t):return a * x
a=1
x0=1
Tmax=1
N=100; dt=Tmax/N;
t=np.linspace(0.0,Tmax,N+1)
x=np.empty(N+1)
x[0]=x0#print(’%f %f’ % (t[0],x[0]))for i in range(N):x[i+1] = x[i] + dt * f(x[i], t[i])#print(’%f %f’ % (t[i+1],x[i+1]))
plt.plot(t,x)
plt.show()N
このプログラムでは、{tn }N
n=0 , {xn }n=0 を配列 (numpy の ndarray) t, x に記憶させて、plot

している。
タイトルや、座標軸が何を表すか等、有益な情報を表示するようにプログラムを書き直すのは読者に任せる。

(前進)Euler 法は精度が低く、数値的安定性についても “平凡な性能” なので、実際に利用されることは少ない。多くのテキストでは、次項で説明する Runge-Kutta 法の利用が勧めら
れている。

11

注意 2.4 (後退 Euler 法) 導関数の前進差分近似の代わりに後退差分近似

x′ (tn ) ≒を採用するとx(tn ) − x(tn−1 )∆tn−1

xn − xn−1= f (xn , tn ).
∆tn−1

これを n = 1, 2, · · · について考えることは、

xn+1 − xn= f (xn+1 , tn+1 )∆tn
を n = 0, 1, 2, · · · について考えることと同じである (単なる書き換え)。分母を払うと(2.13)

xn+1 = xn + ∆tn+1 f (xn+1 , tn+1 ) (n = 0, 1, 2, · · · ).

この漸化式で {xn } を定める方法を後退 Euler 法 (the backward Euler method) と呼ぶ。

xn が既知であるが、まだ xn+1 は既知でないというとき、xn+1 を求めるためには、(2.13)
を xn+1 に関する方程式として解く必要がある。f が簡単な関数の場合を除き、反復法などを行って xn+1 を求めることになる。この点は面倒であるが、後退 Euler 法は、数値的安定性が優れているため、(前進 Euler 法とは異なり) ときどき採用される。— 以上は常微分方程式の場合であるが、偏微分方程式で時間に依存する系のモデルであるものに対しては、後退 Euler
法の採用率はかなり高い。

2.5.2

Runge-Kutta 法

常微分方程式の初期値問題の数値解法のうち、比較的簡単で、高精度な解が得られるため、良く用いられる方法として10 、次の Runge-Kutta 法 (4 次の Runge-Kutta 法, 古典的 Runge-Kutta法とも呼ばれる) がある。これは次の漸化式で {xn } を定めるものである。

1(k1,n + 2k2,n + 2k3,n + k4,n ) ,6(2.14a)xn+1 = xn +(2.14b)k1,n = ∆tn f (xn , tn ),(2.14c)k2,n = ∆tn f (xn + k1,n /2, tn + ∆tn /2),

(2.14d)k3,n = ∆tn f (xn + k2,n /2, tn + ∆tn /2),

(2.14e)k4,n = ∆tn f (xn + k3,n , tn + ∆tn ),1
xn+1 = xn + (k1,n + 2k2,n + 2k3,n + k4,n ) .
6(2.14f)Euler 法にせよ、Runge-Kutta 法にせよ、多くの場合 tn は等分点に取る。つまり [t0 , T ] を
N 等分するとして∆t :=T − t0,N

tn := t0 + n∆t (n = 1, 2, · · · , N )10

長い間定番の方法であったが、最近は常微分方程式を解く部分は自分でプログラムを書かないようになりつつあり、古典的 Runge-Kutta 法が採用されるのは、そろそろ変わる (つまりもっと凝ったアルゴリズムが使われることが増えてくる) のかもしれない。

12で定める (∆tn ≡ ∆t ということ)11 。
以下この節では、特に断りのない限り、解こうとしている区間 [t0 , T ] を N 等分する場合を考えることにする。

C プログラマー向けの説明「常微分方程式の初期値問題を解くプログラムの書き方§ 2 Euler 法, Runge-Kutta 法入門— 1 次元の問題」12 の §2.3.2 Runge-Kutta 法の C プログラム例13 で、C 言語によるプログラ
ム例が見られる。

runge-kutta.c/** runge-kutta.c (Runge-Kutta method for Malthusian model)*/#include <stdio.h>
double a = 1.0;
int main(void){
int i, N;
double t, x, dt, k1, k2, k3, k4;
double f(double, double), x0;
double Tmax;// 初期値x0 = 1.0;// 最終時刻
Tmax = 1.0;// 時間刻み
printf("# N: "); scanf("%d", &N);
dt = Tmax / N;// 初期値
t = 0.0;x = x0;
printf("# t x\n");printf("%f %f\n", t, x);// Runge-Kutta 法
for (i = 0; i < N; i++) {k1 = dt * f(x, t);k2 = dt * f(x + k1/2, t + dt / 2);
k3 = dt * f(x + k2/2, t + dt / 2);
k4 = dt * f(x + k3, t + dt);x += (k1 + 2 * k2 + 2 * k3 + k4) / 6;
t = (i + 1) * dt;
printf("%f %20.14f\n", t, x);}
return 0;}
double f(double x, double t){return a * x;}

解曲線の図を描くための手順は Euler 法のときと同様 (データファイルの名前を変えるくら11

Euler 法も Runge-Kutta 法も Runge-Kutta 型公式と呼ばれる分類に属する (詳しくは専門書を見て下さい)。
Runge-Kutta 型公式は、計算の途中で時間刻み ∆t を変更することが簡単というという長所があり、実際に刻み幅の自動調節がしばしば行われる。
12
https://m-katsurada.sakura.ne.jp/labo/text/intro-ode-simulation/node7.html
13
https://m-katsurada.sakura.ne.jp/labo/text/intro-ode-simulation/node15.html

13い) なので省略する。
必修課題 2.5.2 上のプログラムを実行して、解のグラフ (解曲線) を描いてみよ。
必修課題 2.5.3 a = 1, x0 = 1, t0 = 0, T = 1 とし、初期値問題 (2.5a), (2.5b) を、[t0 , T ] を N
分割して、Euler 法、Runge-Kutta 法で解いたときの計算精度を調べよ。N が増加するにつれて (例えば N = 10, 100, 1000, 10000, · · · で計算してみる)、どのように変化するか (図 2 も見てみよう)。
xN は x(1) = x0 ea(T −t0 ) = 1 · e1(1−0) = e = 2.7182818284 · · · の近似であることを用いる。
以下の 2 つの課題は、プログラムとは直接の関係はなく、紙の上の数学の問題である。
課題 2.5.1 a = 1, x0 = 1, t0 = 0, T = 1 とし、[t0 , T ] を N 分割して、(2.5a), (2.5b) を (前進)
Euler 法、(後退) Euler 法で解いた場合の xN を文字式 (∆t を用いる) で表せ。その結果の妥当性を説明せよ。
(ヒント: 高校の数学 III で学んだ lim (1 + x)1/x = e を思い出すと、N → +∞ のとき xN がx→0

真の解に収束することが納得できるはず。 )課題 2.5.2 a = 1, x0 = 1, t0 = 0, T = 1 の場合に (2.5a), (2.5b) を (前進) Euler 法、Runge2+
Kutta 法で解いた場合の x1 を文字式 (∆t を用いる) で表すと、それぞれ 1+∆t, 1+∆t+ (∆t)
2

4
(∆t)3
+ (∆t)となることを示せ。これから x(t1 ) と x1 の差がそれぞれ O ((∆t)2 ), O ((∆t)5 ) で624あることを説明せよ。

Euler 法、Runge-Kutta 法の次数という重要な概念に触れるために、次の研究課題に取り組むことを勧める。
研究課題 2.5.1 a = 1, x0 = 1, t0 = 0, T = 1 の場合に (2.5a), (2.5b) を (前進) Euler 法、Runge-Kutta 法で解くプログラムを作り (プログラミング言語は何でもよい)、N を 10, 20, 40,· · · に対して xN を求めて、x(T ) − xN が O(N −1 ), O(N −4 ) であることを示せ (横軸 N , 縦軸|x(T ) − xN | で両側対数プロットをすること)。
1"euler.data""runge-kutta.data"

0.01
0.0001
1x10-6
1x10-81x10-10
1x10-12
1x10-1410100

図 2: Euler 法と Runge-Kutta 法の精度の比較
この図 2 をどうやって描いたか、付録 H.2 (p. 98) に説明してある。

14

1000

C プログラマー向けの説明 上の runge-kutta.c の Python バージョンは次のようになる。
runge-kutta.py# runge-kutta.py --- Malthus 方程式を Runge-Kutta 法で解く
# C プログラムの Python バージョン# runge-kutta.c(http://nalab.mind.meiji.ac.jp/~mk/labo/text/ode-workbook/node11.html)# Runge-Kutta 法であるが、応用が効きにくい形。参考プログラム。
def f(x, t):return a * x
a=1
x0=1
Tmax=1
str=input(’N=’); N=int(str);
dt=Tmax/N
t=0.0; x=x0
print(’%f %f’ % (t,x))
for i in range(N):k1 = dt * f(x,t)
k2 = dt * f(x + k1/2, t + dt / 2)
k3 = dt * f(x + k2/2, t + dt / 2)
k4 = dt * f(x + k3, t + dt)
x += (k1 + 2 * k2 + 2 * k3 + k4) / 6
t = (i + 1) * dt
print(’%f %f’ % (t,x))

このプログラムを色々書き換えてみよう。

• 解曲線を描くように。タイトルをつけたり、グリッドを描いたり、描画範囲を 0 ≤ t ≤ 1,
0 ≤ x ≤ 3 にしたり、座標軸を描いたり、軸ラベルをつけてみたり。
• 誤差 x(tN ) − xN = e − xN の値を表示するように。

3logistic 方程式

3.1

Malthus のモデルの修正としての logistic 方程式

マルサスの法則は、個体数の少ない間の増殖は良く説明できるが、個体数が多くなると (環境が悪くなって増殖のスピードが鈍って)、ずれが大きくなる。言い換えるとモデルとしての妥当性が落ちてくる。
ピエール・フランソワ・フェルフルスト (Pierre-François Verhulst, 1804–1849, ベルギー) は、
1838 年に出版した論文の中で、人口予測のために次の微分方程式を提唱した (Wikipedia [3] のHistory 参照)。

(3.1)

dx
= (k − λx)x.
dt

ここで k, λ は正の定数である。
次のようにも表される (両者が同等であることは容易に分かるであろう)。
dN
N= rN 1 −
.
(3.2)
dt
K

15

r は内的自然増加率 (intrinsic rate of natural increase, intrinsic natural growth rate), K は環境収容力 (Carrying capacity) と呼ばれる。
(3.1) や (3.2) を logistic 方程式とよぶ。

3.2

初期値問題の解(3.1) に加えて、初期条件(3.3)

x(0) = cを課した初期値問題の解は次式で与えられる。

(3.4)

x(t) =kc.
λc + (k − λc)e−kt

必修課題 3.2.1 (3.1) は変数分離形の微分方程式であるから、よく知られた方法で解くことができる。それを実行して、(3.4) が正しいことを確認せよ。
しばらく (とりあえず) c は実数の定数としておく (人口の場合、負の数はナンセンスであるが)。


細かい話 (解の存在範囲) — 最初はスルーしても良いですk0 ≤ c ≤ の場合、(3.4) の分母は 0 にならないが、そうでない場合は、分母は t∗ :=
λ
1
λc
− log
のとき (のみ) 0 になる。t∗ の符号を調べると、次のことが分かる。
kλc − k(i) 0 ≤ c ≤kの場合、初期値問題の解の存在範囲は R.
λkの場合、初期値問題の解の存在範囲は (t∗ , +∞). t∗ < 0 であるので [0, +∞)λを含む。

(ii) c >(iii) c < 0 の場合、初期値問題の解の存在範囲は (−∞, t∗ ). t∗ > 0 であるので [0, +∞) を含まない (t = +∞ までは解が延長されない)。
(生物の個体数を考えるときは、c < 0 はナンセンスな初期値なので、(iii) の場合に t の範囲に限界があることはあまり問題にされない。)


必修課題 3.2.2 λ, k を適当に決めて、色々な c に対する解のグラフ (横軸 t, 縦軸 x) を 1 つの図に描け。(注: c < 0, c > λk の場合も含めること。)

3.3

解の漸近挙動, 平衡点とその安定性解の公式 (3.4) を認めれば、次のことが成り立つのは容易に分かる。

(a) c = 0 ならば、任意の t ∈ [0, ∞) に対して x(t) = c = 0.
(b) c =kk
ならば、任意の t ∈ [0, ∞) に対して x(t) = .
λλ16

研究課題 3.3.1 「c = 0 のとき、x(t) := 0 (t ∈ [0, ∞)) とおくと、x は logistic equation (3.1)と初期条件を満たすことは簡単に確かめられる。
」この説明を (a) の証明とするには、微分方程式の初期値問題の解の一意性についての定理が必要である。どういう定理か説明し、それを用いて (a) の証明を完成させよ。
この事実 (a), (b) を、
「0 と
言う。

k
は、(3.1) の平衡点 (equilibrium14 point, 平衡解) である」と
λ

必修課題 3.3.1 (常微分方程式の力学系の) 平衡点の定義を調べて説明せよ。

(付録 G.3 に書いておいたので、それを読んで理解すれば良い。)
必修課題 3.3.2

a が微分方程式 dx= f (x) の平衡点であるとき、
dtx(t) := a (t ∈ [0, ∞))とおくと、これは微分方程式の解であることを示せ。
解の公式 (3.4) を認めれば、次のことも容易に分かる。

(i) c > 0 ならば、 lim x(t) =
t→∞

k
.
λ

(ii) c = 0 ならば、 lim x(t) = 0.
t→∞(iii) c < 0 ならば、 lim x(t) = −∞.
t→t∗

このことから次のことが分かる。

(A) 平衡点kは漸近安定である。
λ(B) 平衡点 0 は不安定である。
必修課題 3.3.3 平衡点に対する「安定」、「漸近安定」、「不安定」という用語の定義を調べて、上の (A), (B) が正しいことを示せ。

(やはり付録 G.3 に書いておいたので、それを読んで理解してから取り組むこと。)

14[ì:kw@líbri@m]17図 3: P. F. Verhulst, Notice sur la loi que la population suit dans son accroissement, Correspondance mathématique et physique, Vol. 10, pp. 113–121 (1838) の p. 115. 最近は、多くの原典がネットで少し探しただけで見つけられる。OCR されていれば、機械翻訳 (DeepL 翻訳, Google 翻訳等) にかけて日本語訳を得るのも簡単である。せっかく便利なものが使えるようになったのだから、こまめに検索・調査することを勧める。(2023/7 追記) macOS Venturaから、画像ファイルからでも、OCR によるテキスト情報のコピー＆ペーストができるように
なった。

18

念のため: DeepL 翻訳 https://www.deepl.com/ja/translator

4単振動の方程式

4.1

単振動の方程式とは実数の定数 ω を用いて

x′′ (t) = −ω 2 x(t)(4.1)

と表される微分方程式を単振動の方程式 (simple harmonic motion equation) と呼ぶ。釣り合
いの位置の近傍での微小な振幅な振動現象の多くが、この方程式に帰着される (バネ振り子の振動 mx′′ (t) = −kx(t), LC 回路 (電気回路) の振動 LI ′′ (t) + C1 I(t) = 0, 波動方程式を変数分離して現れる微分方程式, etc.)。
私見では、単振動の方程式は、理工系の大学生が初年級で学ぶ微分方程式のうちでもっとも重要なものである。

単振動の方程式の解4.2

よく知られているように、ω 6= 0 ならば(4.2)

x(t) = C1 cos ωt + C2 sin ωt (C1 , C2 は任意定数)が (4.1) の一般解である (付録の C「定数係数線形同次常微分方程式に対する特性根の方法」の定理 C.3 を見よ)。
三角関数の合成をすると(4.3)

x(t) = A sin (ωt + ϕ) ,A > 0,

ϕ∈Rの形に変形できる。
問

A, ϕ をどのように定めれば (4.3) が成り立つか。

必修課題 4.2.1 特性根の方法により (4.2) を解き、(4.2) が (4.1) の一般解であることを示せ。
必修課題 4.2.2 任意の実数 x0 , v0 を与えたとき、初期条件(4.4)

x′ (0) = v0

x(0) = x0 ,を課した初期値問題の解を求めよ。それが唯一の解であることはどうして分かるか。
課題 4.2.1 k を正定数とするとき、微分方程式

x′′ (t) + kx′ (t) + ω 2 x(t) = 0を減衰振動の微分方程式と呼ぶ。この微分方程式の一般解を求め、解のグラフ (横軸 t, 縦軸x) を描け。k が小さいと、解は単振動のように符号を変えて振動するが、k が大きいと、解は符号を変えずに 0 に収束する。境目となる k を求めよ。

19

4.3

1 階正規形方程式への帰着(4.1) は 2 階の方程式であるが、ベクトル値関数を導入することで、1 階の方程式に変換できる。
これは一般に適用可能な、必ずマスターすべき操作である。(例えば、Euler 法や Runge-Kutta法のような数値解法を適用するにも必要となる。)!
x2 (t) := x′ (t),x1 (t) := x(t),(4.5)

x(t) :=

x1 (t)x2 (t)

とおくと

(4.6)

x′ (t) =x′1 (t)x′2 (t)!
=x′ (t)(x′ (t))′!
=x′ (t)
x′′ (t)!

!
x2 (t)−ω 2 x(t)==!
x2 (t).
−ω 2 x1 (t)

ここで

!
(4.7)

x2
−ω 2 x1

f (x, t) :=とおくと、(4.1) は次のように表せる。

x′ (t) = f (x(t), t).

(4.8)

1 階方程式になったことに注意しておこう。
多くの微分方程式は、以上と同様の手順で (4.8) の形に変換できる。
必修課題 4.3.1 力学で有名な強制振動 (forced oscillation, forced vibration) の微分方程式

x′′ (t) + ω 2 x(t) = F sin ω ′ t(4.9)

を 1 階正規形微分方程式に変換せよ (F , ω ′ は実数の定数、ω ′ は ω と等しいかもしれないし、等しくないかもしれない)。(これは後の課題 4.5.2 を解くためにも必要である。)
必修課題 4.3.2 高校物理で取り上げられることの多い、(空気抵抗を無視した) ボール投げは、
次の連立微分方程式に従う。

mx′′ (t) = 0,
my ′′ (t) = −mg.
ここで m, g は正定数である。(物理的には、t は時刻, x(t) と y(t) は時刻 t における位置を表す座標、m はボールの質量, g は重力加速度 (SI 単位系では、おおむね 9.80 ms−2 ) を表す。)
x(t) := (x(t), y(t), x′ (t), y ′ (t))⊤ とおき、1 階正規系の微分方程式に変換せよ。

4.4

力学系とは

(4.7) の右辺は t を含んでいないので、f を!
(4.10)

f (x) :=x2
−ω 2 x120と定義して、(4.8) の代わりに

x′ (t) = f (x(t)).

(4.11)とすることも可能である。

じれいけい

このように右辺が f (x(t)) と、t があらわに含まれていない微分方程式を自励系の微分方程式(autonomous diﬀerential equation, 自律系と呼ばれることもある)、あるいは力学系 (dynamicalsystem) の微分方程式と呼ぶ。

(力学系という言葉は、有限次元の常微分方程式以外でも使われる (無限次元力学系, 離散力
学系, …)。その点は注意が必要である。)(4.8) と比べると、(4.11) は一般性が低いわけであるが、実際上現れる (興味を持たれる、とでも解釈してもらえばよい) 多くの場合に微分方程式は (4.11) の形をしている。
自励系微分方程式の解 x = x(t) (t ∈ [t0 , T ]) は、x 空間 (ここでは x1 x2 平面) 内の曲線とみなすことが出来る。それを解の軌道 (orbit) と呼ぶ。正確には、x0 に対して、 x(t) を初期
条件 x(0) = x0 を満たす (4.11) の解として、
O+ (x0 ) = {x(t) | t ∈ [0, ∞)}を x0 を通る正の半軌道、

O− (x0 ) = {x(t) | t ∈ (−∞, 0]}
を x0 を通る負の半軌道、

O(x0 ) = {x(t) | t ∈ R}
を x0 を通る軌道と呼ぶ。
問

解曲線と軌道の違いを説明せよ。

軌道を考えるとき、x1 x2 平面を相空間 (phase space) または相平面 (phase plane) とよぶ15 。

4.5

Python, Julia で数値計算

常微分方程式の初期値問題は、素朴な Runge-Kutta 法のプログラムで十分に解ける場合が多いので、そういう取り組みを勧めるテキストが少なくないが (例えばネットで読める桂田 [4])、現在普及している Python や Julia には、常微分方程式の初期値問題の高品質なソルバーが搭載されていて、可視化をするにも便利なので、それらを利用するのを強く勧めたい (私は宗旨替えをしました)。ここでは Python, Julia のサンプル・コードを示す (付録 B 参照)。
課題 4.5.1 減衰振動の方程式

x′′ (t) + kx′ (t) + ω 2 x(t) = 0,x(0) = x0 ,x′ (0) = v0のシミュレーションをするプログラムを作成し、実際にシミュレーションをするプログラムせよ。
課題 4.5.2 強制振動の方程式

x′′ (t) + ω 2 x(t) = F sin ω ′ t,x(0) = 1,

x′ (0) = 0のシミュレーションをするプログラムを作成し、ω 6= ±ω ′ のとき、ω = ω ′ のときをシミュレー
ションしてみよ。(厳密解の求め方は https://m-katsurada.sakura.ne.jp/lecture/kiso4/
kiso4ode.pdf#page=126 (の問題 5 とその 3 ページ後の解答) に載っている。)
15

数学以外の分野では、「位相空間」とよぶこともあるが、数学では「位相空間」は topological space の訳語として使われるので、phase space を「相空間」とよぶのが普通である。

21

図 4: 強制振動 (ω ′ 6= ±1)

図 5: 強制振動 (ω ′ = 1)

5

Lotka-Volterra の方程式 (捕食者被食者の数理モデル)

5.1

Lotka-Volterra 方程式の勉強のススメ

有名な Lotka-Volterra 方程式を学ぼう。
前節までは、微分方程式の解が比較的簡単な式で表現できた (そういうものを解析解と呼ぶことがある)。しかし多くの微分方程式は、解は存在するが、それを簡単な式で表すことはできない。(むしろ、解が式で表せることの方が稀である。)

Lotka-Volterra 方程式も解を簡単な式で表すことはできないが、解軌道の方程式を得ることができる。そのおかげで十分な解析ができる (そういう意味で少しだけ難しい)。解析解が得られない問題に対するアプローチの仕方を 1 つ学ぶ、という意義がある。この節ではその点を強調した説明をしてみよう。
Lotka-Volterra の方程式は色々なテキストで取り上げられているが、ブラウン [5] を勧めて
おく。

5.2

Lotka-Volterra 方程式

複数の生物種の相互作用には、色々な場合があるが、捕食者-被食者の関係にある 2 種の生物種について16 、Lotka と Volterra は独立に次の方程式によるモデルを提案した。

dx= ax − bxydtdy= −cy + dxy.
dt

(5.1a)(5.1b)(注: 捕食者-被食者モデル (Predator-prey model) には、他にも色々ある。)
x(t), y(t) は、それぞれ時刻 t における被食者の個体数、捕食者の個体数を表す。a, b, c, d は正の定数である。
(5.1a), (5.1b) を Lotka-Volterra の方程式 (Lotka-Volterra equations) とよぶ。

(5.2)

x(t) :=

!
x(t)
,
y(t)

f (x) = f (x, y) :=とおくと、(5.1a), (5.1b) は次のように表される。

(5.3)
16

dx
(t) = f (x(t)).
dt

英語では、捕食者は predator, 被食者は prey という。

22

ax − bxy−cy + dxy!

5.3

まずは数値シミュレーション

この問題は有名であるが、いわゆる求積法 (≒式変形) では解くことができない。しかし、コンピューター・シミュレーションは難しくない。
現象数理学科の学生としては、微分方程式を見た瞬間に「これはシミュレーションできるだ
ろうか、どうすればできるか？」と考えて、必要があれば気軽にシミュレーションするようになってほしい。
Lotka-Volterra の方程式のシミュレーションは、単振動の方程式のシミュレーションの真似をすればできる。演習とする。
(プログラム例はネット上にたくさんある。例えば https://m-katsurada.sakura.ne.jp/
knowhow-2021/node26.html の testodeint2.py とか。自分でプログラムを書けるようになりたい場合は、単振動の方程式のプログラムをたたき台にして、Lotka-Volterra の方程式を解
くプログラムに書き換えてみることをお勧めする。)

図 6: Lotka-Volterra の相空間 (解軌道をいくつか描いてみる)

5.4

解を求めて … まず局所解の存在と一意性(理論的な話に慣れいてない人、特に初めて読む場合は、軽く読み流すだけで良いです。一方、この部分は実はほとんどいつでもワンパターンなので、この機会に少し時間をかけて学ぶ
価値は高いと思います。)
コンピューター・シミュレーションで相当のことが “分かる” けれども、例えば周期解になるらしい、ということが分かっても、本当に周期解だろうか？という疑問に答えるには、数学的な考察が必要になる。
ここでは、まず局所解が存在すること、それがただ一つであること (解の一意性、単独性)をどのように確認するか考えよう。
(初期値問題の局所解とは: 解は関数であり、変数のどういう範囲について定義されているかという定義域がある。解の存在を保証する定理の多くは、初期時刻 t0 から十分短い時間の間は、微分方程式を満たす関数が存在する、という形をしている。そのような定義域が狭い範囲に限定された (短い時間の間の) 解のことを局所解と呼ぶ。)

23

微分方程式 (5.3) の右辺に現れる関数 f (x) =ax − bxy−cy + dxy!
は、各成分が多項式なので、C ∞ 級の関数である。
特に f は連続であるから、局所解が存在することが分かる (付録 D.2 の定理 D.1)。長い時間が経っても解が存在するかどうかは、この段階ではまだ分かっていない (後述する)。
また f は C 1 級であるから、解の一意性が成り立つことが分かる (付録 D.4 の定理 D.9)。一意性と言われてもピンと来ないかもしれないが、後で意外と使うことになる。
(Malthus の法則や、logistic equation などのように、微分方程式の解がこういう式で表される、ということが示せたとき、解の公式が得られただけでなく、解の一意性も (多くの場合は)同時に示せたことになる。そのため、具体的な式で解ける問題ばかり調べていると、一意性にあまり注意を払わずに済んでしまう (「そういうことをする必要はあるんですか？」)。しかし、解を表す式が得られない場合は、解の一意性だけでも分かると、色々と有効である。)

5.5

平衡点とその安定性

微分方程式 (5.3) の平衡点を求めよう。f (x) = 0, つまり

ax − bxy = 0,

−cy + dxy = 0を解くことになる。結果は(5.4)

x = (0, 0), c a
,.
d b


(0, 0) はどちらの種もいない場合で、ある意味ナンセンス (興味がない) であるが、 dc , ab は
2 つの種が共存している場合である。念のため意味を述べておくと、
cdx(t) =  a (t ∈ R)bは (5.3) の解である。被食者は捕食者に食べられる分だけ生まれてきて、捕食者は被食者を食べて生まれる分だけ自然死する、というバランスの取れている状態である。
一般に微分方程式の平衡点の安定性を調べるため、線形安定性解析が使われることが多い。
これは微分方程式 x′ (t) = f (x(t)) の右辺に現れる関数 f (x) の (平衡点 x0 における) ヤコビ
行列 f ′ (x0 ) の固有値の実部の符号を調べて判定する、という方法である。より詳しくは付録
G.4 を見よ。
必修課題 5.5.1 ヤコビ行列 f ′ (0, 0), f ′
調べよ。

c a
,
d b



とその固有値を求めて、それらの実部の符号を

f ′ (0, 0) は正の固有値を持つので、(0, 0) は不安定である。

f ′ dc , ab の固有値の実部は、どちらも 0 である。残念ながら、線形安定性解析では、 dc , abは漸近安定であるとも不安定であるとも結論できない。
注意 5.1 資料によっては、「ヤコビ行列の固有値の実部が 0 であるから中立安定である」、と述べているものがあるが (念のため: 安定だが漸近安定でないとき、中立安定という)、間違いである (そういうものが意外と多いのは困ったものだ)。

24

• 定数係数線形微分方程式 x′ (t) = Ax(t) であっても、A の固有値の実部が 0 以下ということから中立安定性は導かれない。
• 線形化安定解析だけで中立安定であることを結論することはできない。(平衡点が、線形化した微分方程式で中立安定であっても、もとの微分方程式では安定でないかもしれない。)

Lotka-Volterra の方程式の場合は、解軌道の方程式を求める (次項で実行する) ことによって、平衡点の中立安定性が証明できる。

5.6

解軌道の方程式, 解が R 全体に延長できることと中立安定性の証明

まず、実は少し問題含みであるけれど、多くのテキストに載っていて、見通しのよい議論を紹介する。
dyy(dx − c)dy−cy + dxy
= dt ==.
dx
dx
ax − bxyx(a − by)dt
これは変数分離形の方程式であるから、いつもの手順で解くことができる。まずa − bydx − cdy =dxyxとしてから両辺を積分してZ

a − bydy =yZdx − cdx.
x

途中を省略して(5.5)

y a xc= 定数.
eby edx

以上の議論は分かりやすいが、いくつか問題がある。

• 割り算をするので、分母は 0 にならないのか？というツッコミがありうる (変数分離形
微分方程式を解くときに良く出て来る話)。これはそれほど難しくなく解決できる (定理
5.2)。
dy
• dx
を考えるということは、y を x の関数とみなしていることになる。そうみなせる根拠は何だろう。実は、それは局所的にしか成り立たないことで、厳密には正しくないことが分かる (数値シミュレーションで求めた解軌道 (図 6 など) を思い出すと、1 つの xに 2 つの y が対応するのが普通であるから)。まあ、発見的考察をしただけであって、うるさいことを言うな、という意見はありうる。実際、(5.5) という式を見つけてしまえ
y(t)a x(t)cば、任意の解 (x(t), y(t)) に沿って by(t) dx(t) が定数であることが、t で微分することで
e
e
示せる。

まず分母 0 の問題を片付けよう。次の定理が成り立つ。



定理 5.2 (解は第 1 象限を飛び出ない) x(0) > 0 かつ y(0) > 0 ならば、任意の t ∈ R に対して (解が t まで存在する限り) x(t) > 0 かつ y(t) > 0 が成り立つ。




25

発見的考察、というのをやろう。x = 0 を (形式的に) (5.3) に代入すると

dx
= 0,
dt

dy
= −cy.
dt

これを眺めれば、x(t) = 0, y(t) = y(0)e−ct という解を発見できる (代入して (5.3) の解である
ことが、すぐ確認できる。直観的には「被食者がいなければ捕食者は死に絶える」ということである。)。
「y 軸の正の部分は解軌道である」ということが分かった。
同様に y = 0 から y(t) = 0, x(t) = x(0)eat という解を発見できる (「捕食者がいなければ被食者は限りなく増え続ける」ということである)。「x 軸の正の部分は解軌道である」ということが分かった。
最後に、原点は平衡点であるから、原点 1 点からなる集合は解軌道である。
課題 5.6.1 定理 5.2 を証明せよ。(ヒント: 第 1 象限は解軌道で囲まれていることが分かった。
〆は背理法と解の一意性定理を用いる。)
以下、x(0) > 0, y(0) > 0 である解だけを考えることにすると、x(t) > 0, y(t) > 0. ゆえに
x(t) 6= 0, y(t) 6= 0 である。微分方程式から

(5.6)

dx
dy
(−cy + dxy) = (ax − bxy)(−cy + dxy) = (ax − bxy)dtdtである。この式を xy で割り算することに問題はない。

a − by dydx − c dx
=
.
y dt
x dt

(5.7)

ゆえに、解が時刻 t∗ まで存在するならば
Z t∗
Z t∗
a − by dydx − c dx
dt =
dt.
y dt
x dt
0
0
置換積分の公式から

Z y(t∗ )0a − bydy =yZ x(t∗ )0dx − cdx.
x∗

これで上の議論に合流できて (簡単のため x(t ), y(t∗ ) を単に x, y と書くことにすると)(5.8)

y a xc= 定数.
eby edxxc y aは、第 1 象限ではつねに正の値を取り (x, y) = (c/d, a/b) のとき最大
eby+dx値を取る。(x, y) を x 軸や y 軸上の点に近づけると、F (x, y) → 0.
0 < L < F (c/d, a/b) を満たす任意の L に対して、
関数 F (x, y) :=(5.9)

F (x, y) = L

は閉曲線を表す (コンピューターで描いてみることを薦める)。とりあえずこのことは認めることにしよう。
(もちろん図を描けば「明らか」である。証明できた人はレポートしよう。)(5.9) が解軌道を表す方程式である。後で、(x(t), y(t)) は周期関数で、方程式の定める閉曲線の周りをぐるぐる回る (ゆえに解軌道は閉曲線に一致する) ことが分かるが、現時点では解が存在する任意の時刻 t で、(x(t), y(t)) がその曲線の上にあることだけが分かっている。
26

課題 5.6.2 (自力で解くのは難しいかも17 ) 解が R 全体で存在する (任意の時刻 t で x(t), y(t)が定まる) ことを示せ。
(ヒント: 色々なやり方があると思われるが、一つの方針を示す。背理法を用いる。解がある有限の T に対して、[0, T ) まで存在するが、T を超えては延ばせないと仮定する。得られた解軌道の方程式から、任意の解軌道は有界であることが分かる。それから、任意の解軌道の上で f が有界であると分かる。それから極限 lim (x(t), y(t)) が存在することが導かれる。
t→T −0

そうすると解は t = T まで延長され、さらに T を超えて延長できることが導かれる。これは背理法の仮定に矛盾する。ゆえに [0, +∞) で解が存在する。— 桂田 [6] §8.3 に少し書いておいた。)
課題 5.6.3 平衡点

5.7

c a
,
d b



は中立安定であることを証明せよ。

解の周期性(c/d, a/b) 以外の第 1 象限内の点では


dx dy
,
6= (0, 0),
dt dt



dx dy
,
dt dt


> 0.

(c/d, a/b) を通らない任意の解軌道で、速さはある正定数以上である。
証明: 第 1 象限内の、平衡点 (c/d, a/b) 以外の解 (x(t), y(t)) に対して、
L := F (x(0), y(0))
で L を定めると、(x(t), y(t)) はつねに閉曲線 F (x, y) = L 上にあり、閉曲線は R2 の有界閉集合であるから、連続関数である速さは閉曲線上で最小値 (当然正の値) を持つ。
すなわち、ある一定値以上の速さで運動する。曲線の長さは有限だから、いつかは (x(0), y(0))に戻って来る。このことから、(x(t), y(t)) は周期関数であることが証明できる。
課題 5.7.1 T > 0, (x(T ), y(T )) = (x(0), y(0)) ならば、(x(t), y(t)) は周期 T の周期関数であることを示せ。
以上の論法は、応用が効く。
課題 5.7.2 単振り子では、鉛直線から測った角度 θ(t) がd2 θg= − sin θ2dtℓという微分方程式に従う (g は重力加速度、ℓ は振り子の糸の長さ、ともに正の定数)。この微分方程式の (ほとんどすべての) 解が周期関数であることを示せ。
(コメント: 楕円関数を用いると、単振り子の運動の厳密解を表すことが出来るが、楕円関数を知らなくても、Lotka-Volterra 方程式と同じ論法で周期解の存在が示せる。)

17

この課題が解けたら、1 つで S 評価級。

27

5.8

解の平均値

平衡点 (d/c, a/b) 以外の第 1 象限内の任意の解 (x(t), y(t)) は、周期関数であることを示した。その周期を T として
Z
Z
1 T
1 T
x :=
(5.10)
x(t) dt, y :=y(t) dt
T 0
T 0
とおく。いわゆる (時間についての) 平均値を意味する。実は

d
x= ,
c

(5.11)y=abが成り立つ。自然できれい、a, b, c, d だけで定まるのは不思議18 、と感じられる (著者個人の感想だけれど、いかがでしょう)。

(5.11) を証明しよう。まず微分方程式から
Z T ′
Z T
Z T
x (t)
ax(t) − bx(t)y(t)(5.12)dt =dt =(a − by(t))dt.
x(t)
x(t)
0
0
0
一方、周期性より x(0) = x(T ) であるから、左辺について
Z T ′
x (t)
dt = [log x(t)]T0 = log x(T ) − log x(0) = 0.
x(t)0(5.12) に代入して

Z T
0=
0

ゆえに

1
y=
T

x′ (t)
dt =
x(t)

Z T
0

Z T
(a − by(t))dt.
01ay(t) dt =
Tb

Z T
0

a
dt = .
b

同様にして次も得られる

d
x= .
c
(x, y) の平均 (x, y) が平衡点であることは、なるほどと感じる人が多いであろう (平衡点を通る解の平均値が平衡点そのものであることは明らかであるが、そうでない場合にもそうなっている、ということである)。

5.9

1 つのルーツの紹介: D’Ancona の疑問と Volterra の回答
ダ ン コ ナVolterra が、このモデルを考えたのは、そもそも D’Ancona から質問 (相談) されたことがきっかけであった。

D’Ancona の疑問
第 1 次世界大戦中 (1914–1918) に、漁業は小規模にしか行えなかったにもかかわらず、食用に適さない軟骨魚類 (サメやエイなど) が大戦前よりも増え、食用魚は減ってしまったのはなぜか (漁をあまりしなかったのだから増えているのでは？)という D’Ancona の問に対して、Volterra は微分方程式のモデルを作って、その解析に基づいた解答を与えた。これが Lotka-Volterra 方程式の 1 つの由来である。
18

周期解であるから、1 つの解軌道のどこから出発しても同じになるのは明らかであるが、解軌道が違っても、
たとえば平衡点のごく近くをまわっても、平衡点から遠く離れたところをまわっても、平均値がまったく変わらないというのは、それほど明らかではない。

28

これは有名な話で、色々なテキストで紹介されているが、ブラウン [5] の説明が詳しくて分
かりやすく、私のお勧めである。

(数学以外の事実の説明をある程度詳しく書かないと、分かりにくいし、面白くないのではと思われるが、以下に数学的な部分だけ抜書きしておく。[5] の出来の悪い引き写しにしかならない可能性が高いが…)人間が漁業をすることにより、食用魚と軟骨魚を採っている効果を考える。個体数に比例して採ることになると考えて19 、dx= ax − bxy − εxdtdy
= −cy + dxy − εy.
dt(5.13a)(5.13b)に従うと仮定しよう (ここで ε は正の定数)。これは (5.1a), (5.1b) で、a を a − ε, c を c + ε に置き換えたものである。ゆえに §5.8 の議論から、(x(t), y(t)) の平均 (x, y) はc+ε a−ε(x, y) =,dbである。まとめると• 平時 (漁業操業している場合) の平均

c+ε a−ε, bd• 戦争中 (あまり漁業操業していない場合) の平均

c a
,
d b



したがって、戦争中は平時と比べて、食用魚が ( c+ε→ dc ) 減り、軟骨魚類が ( a−ε→ ab ) 増えdbる。ゆえに戦争中は魚類全体に占める軟骨魚類の率が上がるはずである。

1914
11.9%

1915
21.4%

1916
22.1%

1917
21.2%

1918
36.4%

1919
27.3%

1920
16.0%

1921
15.9%

1922
14.8%

1923
10.7%

表 1: イタリアのフューメ港での 1914–1923 年の軟骨魚類の漁獲量率のデータ ([5] から引用)

5.10

Lotka-Volterra 方程式に続いて — 生物の競争系のモデル

今・竹内 [7] には色々なことがまとめられていて大変便利である。

5.10.1競争系とはΩ は Rn の開集合, f : Ω → Rn は C 1 級とする。
(5.14)dx= f (x)dtΩ′ ⊂ Ω で (5.14) が競争系とは(5.15)

i 6= j ⇒∂fi(x) ≤ 0∂xjが成り立つことをいう。
19

本来食用魚が欲しいのに軟骨魚も採ってしまうのは、網を使った漁をしているからだろうか。同じ割合
ε を
c + ε2 a − ε1,
用いているのが気になるが、−εx, −εy をそれぞれ −ε1 x, −ε2 y で置き換えると、(x, y) =とdb
なるだけで、食用魚が減って軟骨魚類が増えるという結論は変わらない。

29

5.10.2(5.16a)(5.16b)2 種競争系 — Lotka-Volterra competition modeldN1= N1 (r1 − β1 N1 − γ12 N2 )dt
dN2
= N2 (r2 − β2 N2 − γ21 N1 )dt

β1 , β2 は種内競争係数、γ12 , γ21 は種間競争係数と呼ぶ。
この力学系の解析については、例えば今・竹内 [7] の 10 章, マレー [8] §3.5 などを見よ。
このモデルに限らず、多くの 2 種の競争系では、例外的な場合を除き 1 種しか存続できないという、競争排除則 (Competitive exclusion principle, Gause’s law of competitive exclusion ともいう [9]) が成り立つと考えられている。
Cf. プランクトンは多種が共存している (Hutchinson [10]) というプランクトンの逆理 (paradox of the plankton)、魚種交替などと比較対象すると良い。

5.10.33 種競争系 — May-Leonard のモデル

May-Leonard は、[11] (1975 年) の中で 3 種の生物が競争するモデルを提唱し、解が興味深い挙動を持つことを示した。

dN1


= N1 (1 − N1 − αN2 − βN3 )


dt



dN2
(5.17)= N2 (1 − βN1 − N2 − αN3 )

dt





 dN3 = N3 (1 − αN1 − βN2 − N3 )dtα, β は正の定数である。
これを Ω := {(N1 , N2 , N3 ) | N1 ≥ 0, N2 ≥ 0, N3 ≥ 0} で考察する。
この力学系の平衡点の安定性、t → ∞ のときの解の漸近挙動について詳しいことが分かっている。シミュレーションをすると楽しめる。
平衡点は 8 つあるが、何種が生き残っているかで分類すると次のようになる。

(a) 0 種平衡点: (0, 0, 0) ∈ Ω.
Jacobi 行列の固有値は 1 (3 重) であり、つねに不安定である。
(b) 1 種平衡点: (1, 0, 0), (0, 1, 0), (0, 0, 1) ∈ Ω
Jacobi 行列の固有値は −1, 1 − a, 1 − b. α > 1 かつ β > 1 ならば漸近安定。α と β のど
ちらか一方が 1 未満ならば不安定。

 
 

β−1
β−1
β−1
α−1
α−1
α−1
(c) 2 種平衡点: αβ−1
, αβ−1, 0 , 0, αβ−1
, αβ−1, αβ−1, αβ−1,0 .
これらが Ω に属するには、α, β が(αβ > 1 ∧ α ≥ 1 ∧ β ≥ 1) ∨ (αβ < 1 ∧ α ≤ 1 ∧ β ≤ 1)を満たすことが必要十分である。このとき、これら 2 種平衡点は不安定である。


1
1
1
(d) 3 種平衡点: α+β+1 , α+β+1 , α+β+1 ∈ Ω.

α + β < 2 ならば漸近安定、α + β > 2 ならば不安定である。α + β = 2 の場合は後述する。
30

まとめ直すと(a) α < 1 かつ β < 11 種平衡点は不安定, 2 種平衡点は Ω に属するが不安定, 3 種平衡点は漸近安定。
(b) α > 1 かつ β > 11 種平衡点は漸近安定, 2 種平衡点は存在するが不安定, 3 種平衡点は不安定。
(c) α > 1 かつ β < 1 かつ α + β > 21 種平衡点は不安定, 2 種平衡点は Ω に属さない, 3 種平衡点は不安定。
安定な平衡点は存在しない。1 種平衡点を巡回する解が存在する。
(d) α > 1 かつ β < 1 かつ α + β < 21 種平衡点は不安定, 2 種平衡点は Ω に属さない, 3 種平衡点は漸近安定。
(e) ((c) の裏返し)
α < 1 かつ β > 1 かつ α + β > 21 種平衡点は不安定, 2 種平衡点は Ω に属さない, 3 種平衡点は不安定。
安定な平衡点は存在しない。1 種平衡点を巡回する解が存在する。
(f) ((d) の裏返し)
α < 1 かつ β > 1 かつ α + β < 21 種平衡点は不安定, 2 種平衡点は Ω に属さない, 3 種平衡点は漸近安定。
(c), (e) の場合に、初期値が平衡点でなければ、3 つの 1 種平衡点 (の近く) を巡回する解が存在する。一見周期解のようにも見えるが、だんだん 1 種平衡点の近くに逗留する時間が長くなる。
上の不等式は等号なしの場合であるが、等号つきの場合にも重要なものがある。特に α+β = 2のとき、N (t) := N1 (t) + N2 (t) + N3 (t), p(t) := N1 (t)N2 (t)N3 (t)とおくと、t → ∞ のとき

p(0)N (0)3が証明できる。N1 + N2 + N3 = 1 と N1 N2 N3 = C の共通部分は閉曲線で、これが limit cycleとなる。
May-Leonard のモデルを修正した研究として Jaramillo-Mrad-Stepien [12] がある。
N (t) → 1,

5.10.4p(t) → C :=魚種交替

魚種交替と呼ばれる現象がある。
うきうお

マイワシ、カタクチイワシ、マサバの 3 種類の浮魚類 (海の表層に住んでいて、プランクト
ン食) の日本近海における個体数は、数十年単位で変動し、沖合ではどれか一種だけが分布を広げている。この優占する魚種が交替する現象を魚種交替 (species replacement) と呼ぶ (河井[13], [14], 川崎 [15])。
魚種交替は、3 すくみの競争関係にある力学系モデル
dN1(5.18a)= c1 + N1 (r1 − a11 N1 − a12 N2 − a13 N3 ),dt
dN2(5.18b)= c2 + N2 (r2 − a21 N1 − a22 N2 − a23 N3 ),dt
dN3(5.18c)= c3 + N3 (r3 − a31 N1 − a32 N2 − a33 N3 )dt

31

で説明できるという仮説が出された (Matsuda-Wada-Takeuchi-Matsuyama [16], Takeuchi-OshimeMatsuda [17], 松田 [18])。
ri は内的自然増加率、aij は種 j が種 i に与える負の影響の強さを表す。それぞれの魚種には、相手に侵されない (競争せずにすむ)「聖域」があり、微小な正の定数である ci は “聖域からの補給” を表す。聖域がないと、低水準期の資源量が減り過ぎて、種の存続が危うくなる。
浮魚３種は沖合の回遊域をめぐっては競争しているが、沿岸域はそれぞれの聖域として残っている。
May-Leonard のモデルでは、1 種平衡点を巡回する解が、時間の経過とともに平衡点に近づくことで、その近傍に長く逗留するようになったが、正の ci を入れることにより、1 種平衡点に近づきすぎなくなる効果が生じ、周期的な解が存在しやすくなる。
このモデルは、May-Leonard のモデルとは色々違っているけれど、重要な違いは聖域からの供給効果を導入しているところと考えられる。そこで、May-Leonard のモデルに聖域からの供給効果だけを加えたモデル

dN1


= c + N1 (1 − N1 − αN2 − βN3 )


dt



dN2
(5.19)= c + N2 (1 − βN1 − N2 − αN3 )

dt





 dN3 = c + N3 (1 − αN1 − βN2 − N3 )dtを解析することが考えられる。

SIR モデル

6(グラフをもっとたくさん含めるべきだったと反省している。)

6.1

感染症の数理モデルとしての SIR モデル

感染症の数理モデルとして一躍有名になった SIR モデルは、Kermack-McKendrick [19] (1927年) で提唱された (その基本的な部分は佐藤 [20] で読解されていて参考になる)。感染症の数
理の定番のテキストである稲葉他 [21] には、COVID-19 で広く知られることとなった西浦博
氏の解説が含まれている。また 2020 年暮れに出版された [21] の増補版には、國谷・稲葉による「COVID-19 の数理モデル解析」という章が加筆されている。
書籍でも論文でもないが、ネット (Newsweek) で公開された
西浦博（北海道大学大学院医学研究院教授）【特別寄稿】「8 割おじさん」の数理モデルとその根拠──西浦博・北大教授
THE NUMBERS BEHIND CORONAVIRUS MODELING
2020 年 6 月 11 日（木）17 時 00 分
https://www.newsweekjapan.jp/stories/world/2020/06/8-39.php
は目を通すことを勧める (有名な「何も対策しなければ 42 万人死亡する」という話の根拠が説明されている)。日本の COVID-19 の実データとの照らし合わせは、神永 [22] に載っている(プログラムが紹介されている)。

SIR モデルでは、全人口を、(1) 感受性者 (susceptible. まだ感染しておらず、免疫を持たないので、感染者との接触により感染する可能性がある)32

(2) 感染者 (infectious. 感染してまだ治っておらず、接触することで感受性者を感染させる可能性がある)(3) 除去者 (removed. 以前は感染していたが、回復したり死亡したり隔離されたりして、もう人に移す可能性がなくなった、あるいは免疫を持つ者— オリジナルのモデルでは回復者(recovered) とされていた)の 3 つに分類する。
各時刻 t におけるそれらの人数を S(t), I(t), R(t) と表すことにする。適当な仮定をおくと(これについては調べ物課題とする)、次の微分方程式を満たすことが分かる。

(6.1a)(6.1b)(6.1c)dS= −βSI,
dt
dI= βSI − γI,
dt
dR
= γI
dt

ここで β, γ は正の定数である。例えば、S, I, R の単位を人、時間の単位を日とすると、β の単−1−1−1
位は 人 日 , γ の単位は 日 である。β のことを (感染) 伝達係数 (transmission coeﬃcient),
γ のことを回復率 (recovery rate) あるいは隔離率という。γ は “平均感染性期間” (D と書かれることが多い) の逆数である。
注意 6.1 (平均感染性期間の意味) D := 1/γ をなぜ平均感染性期間と呼ぶか、私自身は詳し
い説明を読んだことがない。以下は私はこう解釈している、というものである。
既に述べたように。感染者数 I(t) は、微分方程式

dI= βSI − γIdtに従う。ある瞬間から感染が起こらないとすると、例えば完全なロックダウンをして感受性者が隔離されて S = 0 になるとか、完全な感染対策ができて β = 0 になるとか仮定すると

dI
= −γI
dt
という微分方程式になる。初めからこの微分方程式に従うと仮定すると (あるいはある時刻からこの微分方程式に従うとして、その時刻を 0 と置き直せば)

I(t) = I(0)e−γtとなる。単位時間に治る人数 (単位は例えば [人/日]) は、時刻 t で −I ′ (t) = I(0)γe−γt である。時刻 t で ([t, t + ∆t] でというべきか？) 治る確率はこれに比例するであろうが、Z ∞γe−γt dt = 10であるから、実は γe−γt が確率密度関数である。治る時刻 (=感染していた時間) の期待値 (平均) は、t との積の積分でZ ∞1
t · γe−γt dt = .
γ0(これはつまり、指数分布の平均の話ということだろうか。)
33

6.2

数値シミュレーション

実は SIR モデルの厳密解が求まるという話があるが (Harko-Lobo-Mak [23])、あまり使いやすい式ではないそうである。
ここでは数値シミュレーションをしてみる。ほとんど Lotka-Volterra 方程式と同様に扱うことができる。Lotka-Volterra の方程式のプログラムがあれば、それをたたき台にしてプログラムを作ることは簡単で良い練習課題である。
桂田 [4] の付録 G「専用の関数を使ってみる – Python, Julia」20 に Python によるプログラ
ムがある。

図 7: 無次元化 SIR モデルと 4 つの可視化 R0 = 2.5 (桂田 [4])S,I,R の時間変化, 2 次元解軌道, 流行曲線, 3 次元解軌道
神永 [22] の §5.2 にも Python のプログラムが載っている。[22] では、さらに SIR モデルを改良した SEIR モデルを用いて、COVID-19 の日本における第一波のシミュレーションを行
い、実データと比較対照している。

6.3

数学的解析 (1)

数学的解析もしてみよう。やってみると、Lotka-Volterra 方程式の時と同じようなところも多い。

Lotka-Volterra の方程式の場合と同様に、局所解の存在、解の一意性はすぐ分かる (簡単なので省略する)。また解が R 全体に延長できることも同様にして証明ができる (以下で示すよ
うに解軌道の方程式が求まるので、それから解の有界性の議論を経て証明できる)。


定理 6.2 解が 0 を含む区間 I で存在すると仮定すると、

(∀t ∈ I) S(t) + I(t) + R(t) = S(0) + I(0) + R(0).

証明d(S(t) + I(t) + R(t)) = S ′ (t) + I ′ (t) + R′ (t) = −βS(t)I(t) + βS(t)I(t) − γI(t) + γI(t) = 0dtであるから S(t) + I(t) + R(t) は定数である。
20

https://m-katsurada.sakura.ne.jp/labo/text/intro-ode-simulation/node66.html

34

総人口は定数であるから、これを N とするとき、(6.2)

S(t) + I(t) + R(t) = N.

SIR モデルの最初の 2 つの方程式
dS= −βSI,
dt
dI= βSI − γIdt

(6.3a)(6.3b)

は R = R(t) を含まないことに注意しよう。初期条件(6.4)

S(0) = S0 ,I(0) = I0

を与えれば、S(t), I(t) が定まる。すると

R(t) = N − S(t) − I(t)で R(t) が求められる。以上より、S(t), I(t) のみを考えれば十分である。



定理 6.3 (1) (S, I) が、(6.3a)-(6.3b) の平衡点であるためには、I = 0 が必要十分である。

(2) S(0) = 0 であれば S(t) = 0, I(t) = I(0)e−γt . I 軸の正の部分、負の部分はそれぞれ解軌道である。


必修課題 6.3.1 定理 6.3 を証明せよ。


定理 6.4 (SIR モデルの解の性質) 初期値 (S(0), I(0)) が S(0) > 0, I(0) > 0 を満たすと仮定する。

(1) 任意の時刻 t で S(t) > 0, I(t) > 0.
(2) S(t) は狭義の減少関数である。
(3) 解 (S(t), I(t)) は R 全体で存在し、次式が成り立つ。
(6.5)

I(t) − I(0) =

γ
S(t)
log
− S(t) + S(0).
βS(0)(これは解軌道の方程式である。)(4) I(0), S(0) を既知の定数と考えると、I は S の関数とみなせる。そしてdI
γ − βS=.
dS
βS
ゆえに、I は S の関数として、S < βγ の範囲で狭義の増加関数、S > βγ の範囲で狭義の減少関数であり、S = βγ で最大値を取る。


特に、S(0) > βγ ならば、I は S の関数として狭義の増加関数、時刻 t の関数として狭義の減少関数である。S(0) < βγ ならば、t が小さい (S(t) > βγ が成り立つ) うちは、I は増加関数
で、t が大きい (S(t) < βγ が成り立つ) ときは、I は減少関数である。

35

この定理の証明は、ほどほどの難しさであろう。ということで、演習課題にしていたのだけ
れど、誰も解いてくれないので、以下に証明を示す。
証明(1) 第 1 象限の境界は {(S, 0) | S ≥ 0} ∪ {(0, I) | I > 0} である。S 軸上の点はすべて平衡点であり、I 軸の正の部分は解軌道であるから、第 1 象限内の点から出発した解は第 1 象限にとどまる。
もう少し詳しくいうと、第 1 象限の点 (S0 , I0 ) から出発した解 (初期条件 (S(0), I(0)) =(S0 , I0 ) を満たす微分方程式の解) が、第 1 象限の補集合に到達するならば、中間値の定理によって、S 軸上の点 (S ∗ , 0) (ここで S ∗ ≥ 0) か、I 軸上の点 (0, I ∗ ) (ここで I ∗ > 0) を通る。
すなわち、ある t∗ ∈ R が存在して (S(t∗ ), I(t∗ )) = (S ∗ , 0) または (S(t∗ ), I(t∗ )) = (0, I ∗ ).
前者の場合、(S a st, 0) 軸上の点が平衡点であることから、解の一意性に反する。また後者∗については、(S(t), I(t)) = (0, I ∗ e−γ(t−t ) ) が解であることから、やはり解の一意性に反する。ゆえに第 1 象限内の点から出発した解は第 1 象限内にとどまる。

(2) 解軌道は第 1 象限に含まれることがわかった。一方d(S(t) + I(t)) = −γI(t) < 0dtであるから、S(t) + I(t) ≤ S(0) + I(0). S(t), I(t) > 0 に注意するとS(t) ≤ S(0) + I(0),I(t) ≤ S(0) + I(0).

以上から解軌道は有界で、f も解軌道上で有界であることが分かる。従って解は R 全体で定義される。

(3) 微分方程式の解であることから
γ 1I (t) = −βS(t)I(t) + γI(t) = βS(t)I(t) −1 +β S(t)



′



γ 1
= S (t) −1 +β S(t)

′

両辺を t について (一時的に変数を τ に置き換えて)、0 から t まで積分すると

Z t
Z t
γ 1
′
I (τ )dτ =−1 +S ′ (τ )dτ.
β
S(τ
)
0
0
右辺は置換積分の公式から次式に等しい:
Z S(t) 
γ1
dS.
−1 +
βS
S(0)
ゆえに

I(t) − I(0) = S(0) − S(t) +

γ
S(t)
log
.
β
S(0)

(4) S(t) > 0, I(t) > 0 が示されたので、S ′ (t) = −βS(t)I(t) < 0. ゆえに S は狭義単調減少関数であり、逆関数 t = t(S) が存在する。ゆえに I = I(t(S)) は S の関数である。
dI
γ − βSdIβSI − γI
dt
= dS
=
.
=
dS
βSI
βS
dt

36

6.4

数学的解析 (2) 無次元化と基本再生産数

有名な基本再生産数 (basic reproduction number) とその意味 (「基本再生産数が 1 より小さいと感染者数は減少して21 、1 より大きいとしばらくの間は増える (流行する)」…— これ
いきちげんりを閾値原理と呼ぶことがある) を説明する。
方程式をいわゆる無次元化して、簡潔な形に変換してから説明する。

e := S(t) ,S(t)Ne := I(t) ,
I(t)
N
R0 :=

(6.6)

e := R(t) ,
R(t)
N

βN
γ

とおく。この定数 R0 は、除去者の数 R と紛らわしい (特に R() の初期値 R(0) を R0 と書きたくなりそう…) が、混同しないように注意する。(本当は「S と I の 2 つだけを考えれば十分」であるから、早めに R を書かないようにしておくのが良いのかも。) R0 を “R nought”(アール・ノート) と読み、基本再生産数 (the basic reproduction number) と呼ぶ。
このとき次式が成り立つ。

dSee= −γR0 SeI,
dt

dIe
e
= γR0 SeIe − γ I,
dt

e
dR
e
= γ I.
dt

さらに時刻変数もet = γtと変数変換すると

dSee= −R0 SeI,
de
t

dIe
e
= R0 SeIe − I,
de
t

e
dR
e
= I.
de
t

以上をまとめておく。



定理 6.5 S(t), I(t), R(t) が (6.1a), (6.1b), (6.1c) の解であるとき、

et := γt,

S(t)
ee
S(
t) =
,
N

I(t)
ee
I(
t) =
,
N

R(t)
ee
R(
t) =
N

と変数変換すると

dSee= −R0 SeI,
de
t

dIe
e
= R0 SeIe − I,
de
t

e
dR
e
= I.
de
t

ただし R0 は (6.6) で定義する。


もともと S, I, R は個体数 (人口) なので単位は例えば「人」、t は時刻なので単位は例えばeeee「日」であるが、S, I, R, t, そして R0 も無次元量であり、単位を持たない22 。
必修課題 6.4.1 R0 が無次元数であることを確かめよ。

21

このことを「流行は起こらない」といい、一時的にでも感染者数が増加することを「流行が起こる」というようである。
22
e I,
e R
e はどれも、人数を人数で割ったものなので単位なし (次元も無次元) である。
例えば、S,

37

e I,
e R
e は、それぞれ感受性者、感染者、除去者の全体 (全人口) に占める率である。COVID-19
S,−1
1
日 であったと言われている (そう仮定してシミュレーションが行わの第 1 波では、γ = 4.8
れた)。この場合、et は 4.8 日を単位として測った時刻ということになる。
以下では、e を省いてdS
(6.7a)= −R0 SI,dt
dI
(6.7b)= R0 SI − I,dt
dR
(6.7c)= I.
dtと表す。これを無次元化 SIR 方程式と呼ぶことにする。
単位を替えただけなので、元の問題と本質的な違いがある訳ではない。次の 2 つの定理は、二つ目の (5), (6) 以外は前の定理で証明ずみとして良い。


定理 6.6 (1) S(t) + I(t) + R(t) = 1.

(2) (S, I) が (6.7a), (6.7b) の平衡点であるためには、I = 0 が必要十分である。
(3) SI 平面の I 軸の正の部分、負の部分はいずれも (6.7a), (6.7b) の解軌道である。






定理 6.7 S(0) > 0, I(0) > 0 と仮定する。

(1) 初期値 (S(0), I(0)) が S(0) > 0, I(0) > 0 を満たすと仮定する。
(2) S(t) は狭義の減少関数である。
(3) 解 (S(t), I(t)) は R 全体で存在し、次式が成り立つ。
(6.8)

I(t) − I(0) =

S(t)
1
log
− S(t) + S(0).
R0S(0)(4) I(0), S(0) を既知の定数と考えると、I は S の関数とみなせて
dI
R0 S − 1
=
.
dS
R0 S
I は S の関数として、S < R10 の範囲で狭義の増加関数、S > R10 の範囲で狭義の減少関数であり、S = R10 で最大値を取る。
(5) (a) R0 ≤ 1 ならば、任意の初期値に対して、I(t) は減少関数である。
(b) R0 > 1 の場合、初期値が R10 < S(0) (< 1) を満たすならば、時刻が小さいうちはI(t) は増加するが、S(tm ) = R10 を満たす tm より t が大きい範囲では I(t) は減少関数である。
(6) 任意の初期値に対して lim I(t) = 0.
t→∞



証明



(1)–(4) については、定理 6.4 の (1)–(4) を書き直すだけである。

(5) (a) S は減少関数であり S(0) ≤ 1 であるから、任意の t ∈ R に対して R0 S(t) − 1 ≤R0 · 1 − 1 = R0 − 1 ≤ 0. I(t) > 0 であるから I ′ (t) = (R0 S(t) − 1)I(t) ≤ 0. ゆえに Iは減少関数である。
38

1< S(0) が成り立つとき、R0 S(0) − 1 > 0. 連続性より、ある ε > 0R0が存在して、R0 S(t) − 1 > 0 (|t| < ε). このとき I ′ (t) = (R0 S(t) − 1)I(t) > 0 であるから (−ε, ε) で I は増加関数である。
1を満たす tm の存在を認める。) S は狭義減少関数なので(とりあえず S(tm ) =R0t > tm に対して、S(t) < S(tm ) = R10 . ゆえに R0 S(t) − 1 < 1 − 1 = 0. ゆえに
I ′ (t) = (R0 S(t) − 1)I(t) ≤ 0 であるから、I は減少関数である。

(b) R0 > 1 かつ

(c)
o
n
S
= 0 ∩ {(S, I) | I ≥ 0} 上に(6) (S(t), I(t)) は曲線 (S, I) I − I(0) = S(0) − S + R10 log S(0)載っている。ゆえに、ある δ1 > 0 が存在して、すべての t ∈ R に対して S(t) ≥ δ1 .

S(t) は単調減少で、S(t) > 0 (t ∈ R) であるから、 lim S(t) が存在する。
t→∞I(t) も最初のうちは単調増加という場合もあるが、あるところから先は単調減少に切り替
わり、やはり I(t) > 0 (t ∈ R) であるから、 lim I(t) が存在する。
t→∞

実は lim I(t) = 0 が成り立つことを背理法で証明しよう。 lim I(t) 6= 0 と仮定すると、あ
t→∞

t→∞

る正数 δ2 が存在して、すべての t ∈ R に対して I(t) ≥ δ2 . これからk(S ′ (t), I ′ (t))k ≥ |βS(t)I(t)| = β |S(t)| |I(t)| ≥ βδ1 δ2 > 0.
これから解軌道の長さが ∞ であることが導かれ、矛盾が生じる。ゆえに lim I(t) = 0.
t→∞

メモ: 十分時間が経過すると S(t) < βγ となることの証明: 背理法を用いる。S(t) ≥ βγ (t > 0)と仮定すると、I ′ (t) ≥ (β βγ − γ)I(t) ≥ 0. ゆえに I(t) ≥ I(0). これから S ′ (t) ≤ −β · βγ I(0) =−γI(t). ゆえに S(t) ≤ S(0)e−γt . これは S(t) ≥ βγ と矛盾する。
問 3. 感受性者ばかりの集団にごく少数の感染者が時刻 t = 0 で侵入した場合、S(0) = 1 − ε,I(0) = ε,R(0) = 0,

I =1−S+

1
log S
R0

0 < ε  1.

この場合の解軌道の方程式はとみなせる。I 軸との共有点の S 座標は方程式1−S+

1log S = 0R0を満たす。この方程式は 2 つの実数解を持ち、その 1 つは 1 である。R0 > 1 の場合に、1 以外の実数解を求めるプログラムを作成し、それを用いて、流行が終息したとき (I = 0) の 1 − Sを求めてグラフにせよ。
以下の (6.9) という方程式を掲げている本も多いので、問にしておく。
問 4. 人口のほぼ全体が感受性者である集団に、時刻 0 で少数の感染者が侵入した場合、流行が収束するまでに感染した人の割合を z (感染しないで感受性者のままでいる割合が 1 − z) と
すると、

(6.9)

e−R0 z = 1 − zが成り立つことを示せ。

39

図 8: 1.1 ≤ R0 ≤ 5 のとき、流行が終息するまでにどれだけ感染するか
独白 6.8 R0 を与えたとき、z を求めると考えると方程式を解く話になるが、逆に z を与えて対応する R0 を求めるのは、

R0 = −(6.10)

log(1 − z)zという単純な計算だ。右辺の関数のグラフを描くのは簡単で、図 8(横軸 R0 , 縦軸 z) の代わり
になる。

挑戦課題

7

(準備中)

7.1

三種競争系

単に「Lotka-Voltarra の方程式」というと、

7.2

渦糸の力学系

7.3

メトロノームの同期現象

8

参考書案内

常微分方程式は非常に豊富な内容があるため、ほとんどのテキストは、説明する内容を絞っている。その絞り方は千差万別である。そのため、自分の目的にあったテキストを書名だけか
ら見つけ出すのは意外と手間がかかる。自分で探す場合は、図書館か大型書店に行って色々見比べることを勧める。

40

最近は電子図書が普及して来て、図書館でオンラインで読めるものもあるので、必要に応じて辞書代わりに使うことも可能になってきた。
さすがに、簡単な微分方程式の初等的な解法は、どんな本を読んでも大して変わらないだろう (「常識だ」という理由で省略されるおそれはあるけれど)。ネットでちょっと確認するという目的には、例えば桂田 [6] がある。明治大学学生の場合は、丸善 eBook Library23 で「微分方程式」を検索すると、何冊か本がヒットして閲覧可能である。
(「桂田研資料室」24 には、ゼミ向け (一般公開していない) の資料も置いてある。)
以下、読んで学ぶ、という観点から何冊か紹介する。
いわゆる教科書的な本として、例えば笠原 [24] を勧めておく。定理とその証明がきちんと書かれているが、力学系の相図もたくさん掲載されている。著者は多くの良質なテキストを著している人で、微分方程式については笠原 [25] という読み物 (こういう本は珍しい) もある。
最近知った竹之内 [26] (新しい本という訳ではない) も、似た感じの本である。こちらは安価な文庫版があり、持ち運びに便利という利点がある。
ブラウン [27], [5] (二巻本) は、微分方程式の基本的な事項を分かりやすく解説してあるだ
けでなく、有名な現象のモデルである微分方程式を数多く取り上げている。また、初等的な数値解法についても述べられている (プログラミング言語の選択が今風ではないけれど)。基礎と応用のバランスがよいと思う。私自身はやったことはないが、数学科の学部のゼミで採用している先生を見たことは何度かある。ゼミで読みたいという場合は相談に応じる。

[27], [5] は分量が多いので通読に向かないかもしれない。もう少し薄い本を 1 冊きちんと読み切りたいという人には、俣野 [28] を推奨する。力学系の記述も多く、現象数理学科の学生に向いていると思う。
以下は、特定のテーマについて調べ物をする場合に参考になるものを紹介する。
解の存在や一意性など、基礎的なことは多くのテキストで言及されているが、私自身が学んだ経験から言うと、説明が分かりやすくないとか (読み手の実力不足だと言われると辛いよね)、実際にそれを適用しようとすると使いにくいとか、意外と苦労させられる。今回手元にある本を見直して、高橋 [1], [29] は良さそうと感じた。どちらも丸善 eBook に入っていて、明治大学図書館ではオンラインで読むことができる。それから、私が学生の時に参考書として紹介されることの多かったコディントン・レヴィンソン [30] は、豊富なページ数 (二巻本) の中で重要なことをゆったりと解説してある本で、あらためてさすがだと感じた。
常微分方程式の数値解析のテキストとしては、三井 [31], [32], 三井・小藤・齊藤 [33], さらに専門書であるハイラー・ネルセット・ヴァンナー [34], [35], Hairer-Lubich-Wanner [36] をあげておく。数値解法の公式をどうやって導くかという原理的な話については、筆者には一松 [37]が分かりやすかった (微積分の本に書いてあるのは意外だけれど面白い)。
プログラムが載っている本については、お勧めが難しいが (プログラミング言語の好みとかあるだろうから)、神永 [22], 小川・宮路 [38] を推しておく (どちらも Python プログラムが載っている)。
入門段階を終えて、これらの本では不足する場合は「自分の目的に応じて探して下さい」となる。以下は参考まで。
比較的新しい坂井 [39] は非常に内容が豊富であり (その点については比較する本が思いつかないほどである)、それが大きな長所である。欠点はそれが 300 ページ程度に詰め込まれてい23
24

https://elib.maruzen.co.jp/elib/
https://m-katsurada.sakura.ne.jp/labo/library/ode/41るため、やや言葉足らずで読みにくいかもしれない、ということ。ネット評では誤植の多さも指摘されている (まだ出て間もないし…)。ある程度地力のついている人向けかもしれない。
(私自身は、最初はイマイチと思ったが、何度かめくっているうちに、「へー、そうなんだ」が多く、評価が上がりつつある。) これも丸善 eBook に入っているので、オンラインで読むことが可能である。辞書的な利用が可能かもしれない。

A

問題の解答

解答 1. (解答) 次の 2 つのことが成り立つことを意味している。なお K は R または C のど
ちらかを表す (考察の対象を実数値関数に限るならば K = R)。

(a) 任意の C ∈ K に対して、x(t) := Ceat で定義した x(t) は方程式 (2.1) を満たす。
(b) 方程式 (2.1) を満たす任意の x(t) に対して、ある C ∈ K が存在して、x(t) = Ceat が成り立つ。
以上の解答の要点として、マルサスの微分方程式は線形常微分方程式であることがあげられる。

(解説) 一般解という言葉は、常微分方程式のテキストには必ず現れる言葉であるが、「どの微分方程式に対しても適用できる」という意味で一般的に定義しようとすると、意外と難しい言葉である。最低限、次の 2 点を理解すべきである (線形とそうでない場合に分けて考える、という 1 つのことなのかもしれない)。
(a) (線形常微分方程式の場合はわかりやすい定義ができる) 線形同次常微分方程式の場合、解全体の集合 (解空間と呼ばれることがある) は線形空間である、という定理が成り立つ。例
えば
y (n) + a1 (x)y (n−1) + · · · + an−1 (x)y ′ + an (x)y = 0(♥)という微分方程式に対して、ある関数 φ1 , · · · , φn が存在して、(♠)(♥) の解全体の集合 = {C1 φ1 + · · · + Cn φn | C1 , · · · , Cn ∈ C}が成り立つ。この事実を(♥) の一般解は
y = C1 φ1 + · · · + Cn φn(C1 , · · · , Cn は任意定数)

である。




と表現する。厳密な意味は (♠) という集合の等式であるから、次の 2 つが成り立つ、ということになる。

(i) 任意の C1 , · · · , Cn に対して、y = C1 φ1 + · · · + Cn φn とおくと、y は (♥) を満たす。
(ii) (♥) を満たす任意の y に対して、ある C1 , · · · , Cn が存在して、y = C1 φ1 + · · · + Cn φnが成り立つ。

42

(ここまでを正確に理解できれば、「一般解というのは、微分方程式のすべての解を 1 つの式で表すものである」、ということもできるだろう。もっとも知らない人にそう言って正確に伝わるかどうか。定義ではなく、1 つの説明にすぎないと考えるべきだろう。)
以上は単独 n 階の同次微分方程式の場合であるが、1 階連立の同次微分方程式

y ′ = A(x)y(y = y(x) がベクトル値関数で、A(x) が行列値関数)の場合も同様である。また、線形非同次常微分方程式

y (n) + a1 (x)y (n−1) + · · · + an−1 (x)y ′ + an (x)y = b(x) あるいは y ′ = A(x)y + b(x)については、有名な
非同次方程式の一般解 = 同次方程式の一般解 + 非同次方程式の任意の特解という公式が成り立つ (知っている人も多いだろう、万一知らなくてもとりあえず大きな問題ないだろう、ということで、詳しい説明は略する)。この場合も、常微分方程式のすべての解を 1 つの式で表現できる。

(b) 線形でない常微分方程式の場合でも、多くの場合に、解全体の集合の「大部分」をひとまとめに式で表すことができる。そのとき，そのような形で表された解を一般解と呼ぶ。その一般解にまとめることが出来なかった、例外的な解を特異解と呼ぶ。一般解は方程式の階数と同じ個数の独立なパラメーターを含むことが多い (その点については、線形微分方
程式の場合と似ている)。そのパラメーターのことを「任意定数」と呼ぶ。
「大部分」という言葉があるため、(この場合の) 一般解は厳密に定義されているとはいえない (これは線形の場合と大きく異なるところで、注意が必要である)。
特異解の具体的な例を 1 つあげておく。クレローの微分方程式の典型例として有名な
y = xy ′ − (y ′ )2 は
y = Cx − C 2 (C は任意定数)という一般解を持つ (確認: y = Cx−C 2 とするとき、xy ′ −(y ′ )2 = x·C −C 2 = Cx−C 2 = y)が、これ以外に (これら直線解の包絡線になっている)

y=

x2
4

という特異解が存在する (念のため: 特異解は他にもある)。
解答 2. 式変形で方程式を解いて公式が得られる場合、しばしば一意性 (それ以外に解がないこと) の証明にもなっているが、変数分離形微分方程式の解法では、割り算するところがあら
われ、分母が 0 となるケースを除外して計算を進めるので、一意性の完全な証明にはならない。
証明 1x(t) を (2.5a), (2.5b) の解と仮定すると
d −ate x(t) = −ae−at · x(t) + e−at · x′ (t) = −ae−at x(t) + e−at · ax(t) = 0.
dt

ゆえにある定数 C が存在して e−at x(t) = C. これから x(t) = Ce−at であるが、初期条件 (2.5b)に代入すると、C = x0 . ゆえに x(t) = x0 e−at .

43

証明 2 x(t) = x0 e−at が初期値問題の解であることは、代入することで確認できる。
f (x) := ax とおくと、微分方程式は x′ (t) = f (x). また f は Lipschitz 条件を満たす (あるいは f は C 1 級という理由でもよい)。ゆえに初期値問題の解の一意性が成り立つ。
証明 3

(準備中)

解答 3. Mathematica が利用できるならば、関数 FindRoot[] を使えば簡単に済む。


Mathematica



i[s_, R0_] := 1 - s + 1/R0 Log[s]lasts[R0_] := s /. FindRoot[i2[s, R0] == 0, {s, 0.1}]
Plot[1 - lasts[R0], {R0, 1, 5}, PlotRange -> All]

他のプログラミング言語では、例えば Newton 法を使うとして、初期値の選び方に試行錯誤するかもしれない。

44



Python プログラム# test-newton3.py
import numpy as np
import matplotlib.pyplot as plt
def newton(f, df, x0, eps, args):x = x0for i in range(10):dx = f(x,args=args) / df(x,args=args)
x = x - dx
# print(’ Δ x=%e, x=%g, f(x)=%e’ % (dx, x, f(x)))
if abs(dx) < eps:
return x
print(’newton: 収束しませんでした。初期値=%g, 修正量=%e’ % (x0,dx))
return x# 解きたい方程式 f(x)=0 の f()
def I(S,args=(2.5,)):R0=args[0]return 1.0-S+np.log(S)/R0# f() の導関数
def dI(S,args=(2.5,)):return -1.0 + 1.0/(R0*S)# 感染せずに残る感受性者の割合を求める
n=40
R0s=np.linspace(1.1,5.0,n+1)
left=np.empty_like(R0s)for i in range(n+1):
R0=R0s[i]if i==0:left[i]=newton(I,dI,0.8,1e-15,args=(R0,))
else:
left[i]=newton(I,dI,0.5*left[i-1],1e-15,args=(R0,))
plt.plot(R0s,1.0-left)
plt.title(’ 基本再生産数と最終的に感染した人の割合’)
plt.grid()
plt.xlabel(’R0’)#plt.ylabel(’Those who remained uninfected’)plt.ylabel(’Those who have infected’)plt.savefig(’test-newton3.png’)plt.show()#for i in range(n+1):# print(left[i])


解答 4.

I(t) − I(0) =
であるから

0−0=

S(t)
1
log
− S(t) + S(0)
R0
S(0)

1−z
1
log
− (1 − z) + 1.
R0
1

整理すると

−R0 z = log(1 − z).
これから e−R0 z = 1 − z.

45



B

常微分方程式の初期値問題の数値シミュレーションの手引き

B.1

はじめに

常微分方程式の初期値問題は、よくある問題で、色々なプログラミング言語 (処理系) でシ
ミュレーションできる。
現象数理学科で 3 年ゼミを受講している人は、すでに常微分方程式の初期値問題の (数値)シミュレーションをしたことがあるはずである。そのときにマスターして気に入っている方法があるならば、それを使って構わない。
私はずっと以前は、Euler 法や Runge-Kutta 法を用いた C 言語のプログラムを書く方法を学生に勧めていたが (やり方はネットで見られる桂田 [4] で説明してある)、その後は Java やRuby を使ったりして、最近のゼミでは Mathematica や Python を使って卒業研究をした学生が多い。
偏微分方程式を解くことを考えると、ゼミの中で Julia の利用を増やしたいと考えているが、常微分方程式ならば Python で十分というか、もしかすると今一番のお勧めかもしれない。
ここでは、Python と Julia での例を見せることにする。
どちらの場合も、最初に解こうとしている微分方程式を一階正規系、つまり(B.1a)(B.1b)

dx(t) = f (x(t), t)dtx(t0 ) = x0の形に変換する。
問題としては、Malthus モデル、単振動の方程式、Lorenz 方程式を取り上げる。これは良く知られた問題であること、相空間の次元がそれぞれ 1, 2, 3 であり、ある種の網羅性がある
ことによる。

B.2

Python によるシミュレーション

Python の処理系は普通はインタープリターなので25 、十分な計算速度が得られるか心配になるが26 、実際には Fortran や C で書かれた高速なライブラリィを呼び出すことで、十分な効率を達成できる。
現象数理学科 Mac では Python が使えるようになっているはずである (何らかの理由でなく27などを参考にしてインストー
なってしまった人は、「Python のインストールの仕方まとめ」
ルすること。(b) Anaconda の導入がとりあえずのお勧め。)。
VS Code や Atom のようなテキスト・エディターでプログラムを書いてもいいし、Jupyter
ノートブックで作業してもよい。
常微分方程式のライブラリィは、定評のある数値計算ライブラリィである scipy に含まれる
scipy.integrate を利用するのが良いと思う。scipy.integrate の公式ドキュメントは25
多くのプログラミング言語処理系は、コンパイラーとインタープリターに分類される。その言語で書かれたプログラムを実行可能な機械語プログラムに翻訳することをコンパイルといい、翻訳するソフトウェアのことをコンパイラーと呼ぶ。また、その言語で書かれたプログラムを解釈して実行するソフトウェアをインタープリ
ター (普通は通訳という意味) と呼ぶ。C 言語や C++ 言語の処理系は、コンパイラーであるのが普通である。
一方、Python 処理系はインタープリターである。Java は、Java 言語で書かれたプログラムを中間言語にコンパイルした上で、中間言語を実行できるプログラム (JVM — Java Virtual Machine) で実行するので、コンパイラーとインタープリターの中間であると言える。
26
その理由から、私は Julia や Crystal に注目している。これらは JIT (just in time) compilter を用いることで実行効率を高めている。
27
https://m-katsurada.sakura.ne.jp/knowhow-2021/node70.html

46Integration and ODEs (scipy.integrate)
https://docs.scipy.org/doc/scipy/reference/integrate.html
で読める。
少し前に神永 [22] という本が出版された。色々参考になることが書いてある。

B.2.1

Malthus モデル

Malthus モデルの初期値問題
dx(t) = ax(t),dtx(0) = x0

(B.2)
(B.3)

を解いてみよう。ここで説明するやり方は、多くの 1 次元問題に適用できる (例えば logistic
方程式が解ける)。

scipy の odeint(func, y0, t [,args,Dfun,col deriv, · · · ]]) という関数を使って初期値問題(B.4a)(B.4b)dy= f (y, t) (t ∈ [t0 , T ])dty(t0 ) = y0が解ける。
念のため: ODE は ordinary diﬀerential equation (常微分方程式) の頭文字を取って作った
略称である。

odeint() の利用 odeint には有名な ODEPACK28 というライブラリィが使われているそうである。ネット上に非常に詳細な ODEPACK の情報がある。
• func は f に対応する。
• y0 は初期値 y0 に対応する。
• t は “時刻” であるが、詳しくは [t0 , T ] を分割したときの分点 t = {tn }N
n=0 (を要素に持29
つ ndarray ) である:
t0 < t1 < · · · < tN = T.
t=np.linspace(t0,T,N+1) とすれば得られる。
詳しいことは scipy.integrate.odeint — ScyPy 1.11.3 Manual を見よ。
微分方程式 (B.1a) の右辺に現れる f を計算する関数 func を作れば、簡単に解が求まる。
次のプログラムでは、malthus() というのが f を計算する関数である。f (x, t) = ax という簡単な関数なので、実質 2 行だけである。

28

https://www.netlib.org/odepack/ndarray は Python の Numpy で定義された用語である。https://docs.scipy.org/doc/numpy-1.13.0/
reference/arrays.ndarray.html
29

47





malthus1.py# malthus1.py --- dx/dt=a*x, x(0)=x0 を解いて解曲線を描く
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
def malthus(x, t, a): # f(x,t, パラメーター)
return a * x
a=1
t0=0; T=1; n=10t=np.linspace(t0,T,n+1)
x0=1
sol=odeint(malthus, x0, t, args=(a,))
plt.plot(t,sol)# 解曲線描画

plt.ylim(0, 3)# 縦軸の範囲指定# 解を求める# タイトル、凡例の位置、横軸と縦軸の説明、グリッド
plt.title(’Malthus: dx/dt=ax, x(0)=x0; a=’+str(a)+’, x0=’+str(x0))
plt.xlabel(’t’)
plt.ylabel(’x’)
plt.grid()
plt.show()Malthus: dx/dt=ax, x(0)=x0; a=1, x0=1

3.0
2.5

x

2.0
1.5
1.0
0.5
0.0

0.0

0.2

0.4

t

0.6

0.8

1.0

図 9: malthus1.py で 1 つの解の解曲線を描いてみる
次の 2 行で、区間 [0, 1] の 10 等分点を計算して、変数 t に記憶している。



t0=0; T=1; n=10t=np.linespace(t0,T,n+1)




48

sol=odeint(malthus, x0, t, arg=(a,)) で解を計算している。 , arg=(a,) というのは見慣れない書き方かもしれないが、a をオプションのパラメーターとして渡す仕組みを使っている (パラメーターが 1 個のときは最後に , をつける必要があるらしい)。
この後、解をどのように出力するかは、色々なやり方がある。ここでは matplotlib.pyplotの plot() を用いて解のグラフ (解曲線) を描いている (グラフ描きは、微分方程式とは直接関係ないので、多くの Python の入門書に説明が書いてある。)。
では次に複数の解の様子を同時に眺めてみる。初期値を変えて解を求めて、その解曲線を描き足していく。
malthus2.py# malthus2.py --- dx/dt=a*x, x(0)=x0 を複数の初期値について解いて解曲線を描く
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
def malthus(x, t, a):return a * x
a=1
t0=0; T=1; n=10t=np.linspace(t0,T,n+1)
for i in range(0,10):x0=i*0.2
sol=odeint(malthus, x0, t, args=(a,))
plt.plot(t,sol)# 解曲線描画
plt.plot(t,sol, label=’x0=’+’{:6.1f}’.format(x0))
plt.ylim(0, 5)# 縦軸の範囲指定# タイトル、凡例の位置、横軸と縦軸の説明、グリッド
plt.title(’Malthus: dx/dt=ax, x(0)=x0; a=’+str(a))
plt.legend(loc=’upper left’)
plt.xlabel(’t’)
plt.ylabel(’x’)
plt.grid()
plt.show()

当たり前のことであるが、解の一意性が成り立つので、解曲線は互いに交わらない。

solve ivp() の利用 scipy には、ODEPACK 由来でない、より新しい関数群も用意されている。例えば solve ivp(func, t span, y0, [,method, t eval,· · · ]]) を使って初期値問題が解ける。
• func は f を計算する関数であるが、引数並びに注意が必要である。solve ivp() のとき
は func(t, x[, param]) という順であり、odeint() のときの func(x, t [, param])とは違っている (t と x が逆)。
• t span は、[t0, T] のような、初期時刻と最終時刻からなるリストである。
区間をどのように分割するかは、ユーザーが指定するのではなく、自動的に決められる。
そのかわり、以下のサンプル・プログラムでは、rtol=, atol= で許容される誤差を指定

49

Malthus: dx/dt=ax, x(0)=x0; a=1

5

x0= 0.0x0= 0.2x0= 0.4x0= 0.6x0= 0.8x0= 1.0x0= 1.2x0= 1.4x0= 1.6x0= 1.8

4

x

3
2
1
0

0.0

0.2

0.4

t

0.6

0.8

1.0

図 10: malthus2.py で複数の初期値に対する解の解曲線を描いてみる(コピーするとモノクロになって分かりにくいけれど…)する。内部で、いわゆる刻み幅の自動調節が行われる。
詳しいことは scipy.integrate.solve ivp — ScyPy v1.11.3 Manual を見よ。

50





malthus3.py# malthus3.py --- dx/dt=a*x, x(0)=x0 を複数の初期値について解いて解曲線を描く
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
def malthus(t, x, a): # 順番に注意 f(t,x, パラメーター)
return a * x
a=1
t0=0; T=1; n=10t=np.linspace(t0,T,n+1)
for i in range(0,10):x0=[i*0.2] # スカラーでない
sol=solve_ivp(malthus, [t0,T], x0, args=(a,),
dense_output=True, rtol=1e-6,atol=1e-8)x=sol.sol(t)# 解曲線描画
plt.plot(t,x.T, label=’x0=’+’{:6.1f}’.format(x0[0]))plt.ylim(0, 5)# 縦軸の範囲指定# タイトル、凡例の位置、横軸と縦軸の説明、グリッド
plt.title(’Malthus: dx/dt=ax, x(0)=x0; a=’+str(a))
plt.legend(loc=’upper left’)
plt.xlabel(’t’)
plt.ylabel(’x’)
plt.grid()
plt.show().BMUIVTEYEUBYY	YBY
Y
Y
Y
Y
Y
Y
Y
Y
Y



Y















U





図 11: malthus3.py で複数の初期値に対する解の解曲線を描いてみる

51

B.2.2単振動の方程式

単振動の方程式の初期値問題(B.5a)

x′′ (t) = −ω 2 x(t),(B.5b)

x(0) = x0 ,x′ (0) = v0

を解いてみよう。これは 2 階方程式であるが、既に説明した方法によって

(B.6)

dx
(t) = f (x(t), t),dtx(0) = x0という 1 階正規形の微分方程式の初期値問題に変換される。ただし
!
x
(t)
1
, f (x, t) = f (x1 , x2 , t) :=v(t) := x′ (t), x1 (t) := x(t), x2 (t) := v(t), x(t) :=x2 (t)!
x2−ω 2 x1

である。
これは 2 次元の微分方程式である点が、Malthus モデルとは違っていて、計算の仕方、可視化の仕方に少し違いがある。それを実例で示そう。
このやり方を応用して、Lotka-Volterra 方程式を解くことも出来る。


shm1.py — scipy.integrate odeint() で解く# shm1.py --- simple harmonic motion#
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt# X=[X[0],X[1]]=(x,x’), f(X)=(X[1],-ω^2 X[0]]def shm(x, t, omega2):
return [x[1], -omega2 * x[0]]omega = 1; omega2 = omega * omegat0=0; T=20; n=100; t=np.linspace(t0, T, n+1)x0=1; v0=0; X0=[x0,v0]
X=odeint(shm, X0, t, args=(omega2,))
plt.title("単振動: x’’=-ω^2 x, ω="+str(omega)+", (x(0),x’(0))=("+str(x0)+’,’+str(v0)+’)’)
plt.plot(t,X[:,0],’b’, label=’x’)
plt.plot(t,X[:,1],’g’, label=’v’)plt.legend(loc=’best’)plt.xlabel(’t’)
plt.ylabel(’x and dx/dt’)
plt.grid()
plt.show()

X1 (t) = x(t) と X2 (t) = x′ (t) のグラフを描いたわけであるが、解曲線 (ベクトル値関数
X(t) = (X1 (t), X2 (t)) のグラフ) は、R3 内の曲線である。それを描くにはどうすればよいかは後述する。
(このプログラムでは、タイトルに日本語を使っている。現時点では、自分で設定しないと日本語が表示されずに文字化けして、警告が表示される。設定法については、付録 H.3 を見よ。設定が面倒な場合は、プログラム内の “単振動” を削ればよい。)
初期値を色々変えた複数の解を求め、それらの解軌道 (解の像) を描いてみよう。

52,

: x''=- ^2 x, =1, (x(0),x'(0))=(1,0)
1.00
0.75

x and dx/dt

0.50
0.25
x
v

0.00
0.25
0.50
0.75
1.00
0.0

2.5

5.0

7.5

10.0
t

12.5

15.0

17.5

20.0

図 12: 単振動: x(t) と v(t) = dx/dt(t) の時間変化

shm2.py# shm2.py --- simple harmonic motion#
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt# X=[X[0],X[1]]=(x,x’), f(X)=(X[1],-ω^2 X[0]]def shm(x, t, omega2):
return [x[1], -omega2 * x[0]]omega = 1; omega2 = omega * omegat0=0; T=20; n=100; t=np.linspace(t0, T, n+1)
for i in range(0,30):x0=0.1*i; v0=0; X0=[x0,v0]
X=odeint(shm, X0, t, args=(omega2,))
plt.plot(X[:,0],X[:,1])
plt.title("単振動: x’’=-ω^2 x, ω="+str(omega)+"の相図")
plt.legend(loc=’best’)plt.xlim(-2,2)
plt.ylim(-2,2)
plt.xlabel(’x’)
plt.ylabel(’dx/dt’)
plt.grid()
plt.show()

解のグラフと解軌道の両方を描くプログラム例も紹介しておく。こちらは (半分気分転換で)scipy.integrate の solve ivp() を使う。
まず、X1 , X2 = v のグラフと、解軌道を描く。
shm graph orbit.py# shm_graph_orbit.py --- simple harmonic motion53

: x''=- ^2 x, =1

2.0
1.5
1.0

dx/dt

0.5
0.0
0.5
1.0
1.5
2.0

2.0

1.5

1.0

0.5

0.0
x

0.5

1.0

1.5

図 13: 単振動: x(t) と v(t) = dx/dt(t) の相図# coding: utf-8#
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt# x=(x1, x2), f(x)=(x2,-ω^2 x_1)
def shm(t, x, omega):
return [x[1], -omega*omega*x[0]]omega=2.5
x0=[1.0,0.0]
Tmax=20.0
sol=solve_ivp(shm, [0.0, Tmax], x0, args=(omega,),dense_output=True, rtol=1e-10,atol=1e-10)
n=2000t=np.linspace(0.0, Tmax, n+1)
x = sol.sol(t)plt.figure(figsize=(15,5))
plt.subplot(121)plt.plot(t,x.T)
plt.xlabel(’t’)
plt.legend([’x’, ’v=dx/dt’], shadow=True)

54

2.0

plt.title(’simple harmonic motion (t-x,v)’)
plt.subplot(122)
plt.xlim(-1.0,1.0)
plt.plot(x[0], x[1], "-")
plt.xlabel("x")
plt.ylabel("v")
plt.title("orbit in xv-plane")plt.show()

図 14: 単振動の解のグラフと解軌道
いわゆる解曲線というのは t 7→ (x(t), v(t)) のグラフのことだから、3 次元空間の中の曲線である。これも描いてみよう。


shm graph orbit v2.py# shm_graph_orbit_v2.py --- simple harmonic motion# coding: utf-8#
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt# x=(x1, x2), f(x)=(x2,-ω^2 x_1)
def shm(t, x, omega):
return [x[1], -omega*omega*x[0]]omega=2.5
x0=[1.0,0.0]
Tmax=20.0
sol=solve_ivp(shm, [0.0, Tmax], x0, args=(omega,),dense_output=True, rtol=1e-10,atol=1e-10)
n=2000t=np.linspace(0.0, Tmax, n+1)
x = sol.sol(t)55fig=plt.figure(figsize=(15,5))
plt.subplot(131)plt.plot(t,x.T)
plt.xlabel(’t’)
plt.ylabel(’x,v’)plt.legend([’x’, ’v=dx/dt’], shadow=True)
plt.title(’simple harmonic motion (t-x,v), ω=’+str(omega))plt.subplot(132)
plt.xlim(-1.0,1.0)
plt.plot(x[0], x[1], "-")
plt.xlabel("x")
plt.ylabel("v")
plt.title("orbit in xv-plane")# SIR 空間での軌道# 3DAxes を追加
ax = fig.add_subplot(133, projection=’3d’)# 軸ラベルを設定
ax.set_xlabel("$t$")
ax.set_ylabel("$x$")
ax.set_zlabel("$v$")#曲線を描画
ax.plot(t, x[0], x[1])plt.title("solution curve")
plt.show()

図 15: 単振動の解のグラフと解軌道

56

B.3

Julia によるシミュレーション

B.3.1

Julia のインストール＆使い方

Julia についてはネットに情報があふれている。「Julia メモ 1.2 インストール」30 には、この文書のプログラムで利用するパッケージ (Printf, DifferentialEquations など) の準備の仕方も説明してある。
とりあえず DifferentialEquations.jl の利用を勧めておく。
Python のときは、テキスト・エディターでプログラムを書いて、それを python なんとか.py のように実行するやり方で十分使えたが、Julia の場合は、あまり勧められない。
私は、以前は、テキスト・エディターでプログラムを書いて、それを REPL (julia コマンド)に include() して利用するやり方を使っていたが、最近は Jupyter を使うのが便利と感じている。Jupyter は Python でおなじみの人が多いだろうか。一石二鳥かもしれない。

B.3.2

Malthus モデル

malthus1.jl

# malthus1.jl --- dx/dt=a*x, x(0)=x0 を解いて解曲線を描く
using DifferentialEquations
using Plotsf(u,a,t) = a*u# x, パラメーター,ta=1.0
t0=0.0; T=1.0; tspan = (t0, T)
x0=1.0prob = ODEProblem(f,x0,tspan,a)sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8);
plot(sol,linewidth=5,title="Malthus: dx/dt=ax; a=1.0, x0=1.0",
xaxis="t", yaxis="x(t)", label="numerical solution") # legend=false
plot!(sol.t, t->x0*exp(a*t),lw=3,ls=:dash,label="exact solution!")savefig("malthus1_julia.png")
savefig("malthus1_julia.pdf")


30



https://m-katsurada.sakura.ne.jp/labo/text/julia-memo/node4.html

57

図 16: Malthus モデル: 1 つの解の解曲線を描く
malthus2.jl# malthus2.jl --- dx/dt=a*x, x(0)=x0 を複数の初期値について解いて解曲線を描く
using Printf
using DifferentialEquations
using Plotsf(u,a,t) = a*u# x, パラメーター,ta=1.0
t0=0.0; T=1.0; tspan = (t0, T)
for i=0:10
x0=i*0.2
prob = ODEProblem(f,x0,tspan,a)sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8);
if i==0p=plot(sol, title="Malthus: dx/dt=ax, x(0)=x0; a=1.0",
xaxis="t", yaxis="x(t)",label=@sprintf("x0=%4.1f", x0), legend=:topleft)else
p=plot!(sol, label=@sprintf("x0=%4.1f", x0))
end
end
# savefig() しないときは次が必要かも
display(current())# 図の保存savefig("malthus2_julia.png")
savefig("malthus2_julia.pdf")




58

図 17: Malthus モデル: 複数の解の解曲線を描く

B.3.3

単振動の方程式

shm1.jl# shm1.jl --- simple harmonic motion#
using DifferentialEquations
using Plots
using Printffunction shm!(y,x,p,t)omega2,=py[1] = x[2]y[2] = - omega2 * x[1]endomega = 1; omega2 = omega * omega; p=(omega2,)t0=0; T=20; tspan=(t0, T)
x0=1; v0=0; X0=[x0,v0]
prob = ODEProblem(shm!,X0,tspan,p)sol = solve(prob, reltol=1e-8, abstol=1e-8);
plot(sol, xlim=(t0, T), ylim=(-1.2, 1.2),
title="simple harmonic motion: dx/dt=y, dy/dt=-ω^2 x, x(0)=1, y(0)=0",
xlabel="t", ylabel="x and y", size=(800,500))
savefig("shm1_julia.png")
savefig("shm1_julia.pdf")




59

図 18: 単振動の方程式: x(t), y(t) のグラフを描く


shm1.jl# shm2.jl --- simple harmonic motion#
using DifferentialEquations
using Plots
using Printffunction shm!(y,x,p,t)omega2,=py[1] = x[2]y[2] = - omega2 * x[1]endfunction draworbits(omega)omega2 = omega * omega; p=(omega2,)t0=0; T=2*pi/omega; tspan=(t0, T)
for i=0:30
x0=0.1*i; v0=0; X0=[x0,v0]
prob = ODEProblem(shm!,X0,tspan,p)sol = solve(prob, reltol=1e-4);
if i==0
plot(sol, vars=(1,2), xlim=(-2, 2), ylim=(-2, 2), aspect_ratio=1,
title="simple harmonic motion: dx/dt=y, dy/dt=-ω^2 x",
xlabel="t", ylabel="x and y", legend=:none)
else
global fig=plot!(sol, vars=(1,2), xlim=(-2,2), ylim=(-2,2),legend=:none)end
end
return fig
end
if abspath(PROGRAM_FILE) == @__FILE__
print("こうやって実行すると時間がかかります。")
draworbits(1.0)
savefig("shm2_julia.png")
savefig("shm2_julia.pdf")
else
println("draworbits(omega)")println("savefig(\"なんとか.png\")")
end


60



図 19: 単振動の方程式: x(t), y(t) のグラフを描く

61

C

定数係数線形同次常微分方程式に対する特性根の方法 (素朴な説明)

ここでは定数係数線形同次常微分方程式に対する特性根の方法を説明する。(非同次方程式の解をどのように求めるかは重要であるが、その説明はさぼる。例えば桂田 [6] などを参照せよ。)
基本的事項なので授業で学んだはずだが、便利なように簡単にまとめて収録しておく、ということである。
この節では、証明はなるべく初等的になるように努めた。そのため詳しい人にはまどろっ
こしく感じられるかもしれない。(初期値問題の解の一意性定理から、n 階線形微分方程式の解空間が n 次元線形空間であることが容易に証明できて、それを使うと色々と見通しがよく
なる。)

C.1

とりあえず要点を学ぶ

p, q を複素数の定数とする。(実数だけで話が完結する場合もあるが、煩雑になるので、まずは p, q が複素数であり、解も複素数値の関数を認めた場合について述べる。)(C.1)

x′′ (t) + px′ (t) + qx(t) = 0.

2 次方程式(C.2)

λ2 + pλ + q = 0

を微分方程式 (C.1) の特性方程式と呼ぶ。特性方程式の根 (解) のことを特性根と呼ぶ。


定理 C.1 (特性根が相異なる場合) p, q を複素数の定数とする。λ2 + pλ + q = 0 が相異なる 2 根 α, β を持つならば(C.3)

x(t) = C1 eαt + C2 eβt(C1 , C2 は任意定数)

は (C.1) の一般解である。すなわち、(a) 任意の定数 C1 , C2 に対して、(C.4) で定まる x は (C.1) の解である。
(b) (C.1) の任意の解は、適当な定数 C1 , C2 を用いて、x(t) = C1 eαt + C2 eβt と一意的に
表される。




カタい言い回しをすると、「(C.1) の解空間 (解全体の集合) は、2 次元の線形空間で、基底(解の基本系) として eαt , eβt が取れる」となる。特に(C.1) の解空間 = C1 eαt + C2 eβt C1 , C2 ∈ C
である。

62





定理 C.2 (特性根が重根となる場合) p, q を複素数の定数とする。λ2 + pλ + q = 0 が重
根 α を持つならば(C.4)

x(t) = C1 eαt + C2 teαt

(C1 , C2 は任意定数)

は (C.1) の一般解である。すなわち、(a) 任意の定数 C1 , C2 に対して、(C.4) で定まる x は (C.1) の解である。
(b) (C.1) の任意の解は、適当な定数 C1 , C2 を用いて、x(t) = C1 eαt + C2 teαt と一意的に
表される。




p, q が実数の場合、特性根 α, β が虚数ならば、それらは互いに共役複素である。すなわち、
実数 a, b (b 6= 0) が存在して
α = a + ib, β = a − ib.
このとき (C.4) が一般解であるのはすでに述べた通りだが、(C.5)

x(t) = C1 eat cos bt + C2 eat sin bt (C1 , C2 は任意定数)

も (C.1) の一般解である。特に C1 , C2 として実数のみを選ぶと、(C.1) の実数値関数の解す
べてを表すことができる。


定理 C.3 (実係数で、特性根が虚数の場合) p, q を実数の定数とする。λ2 + pλ + q = 0 が
a + ib (a, b ∈ R; b 6= 0) を根 α に持つならば(C.6)

x(t) = C1 eat cos bt + C2 eat sin bt (C1 , C2 は任意定数)

は (C.1) の一般解である。
特に(a) 任意の実数の定数 C1 , C2 に対して、(C.6) で定まる x は (C.1) の実数値関数の解である。
(b) (C.1) の任意の実数値関数の解は、適当な実数の定数 C1 , C2 を用いて、x(t) =C1 eat cos bt + C2 eat sin bt と一意的に表される。



以上の話は、一般の自然数 n に対する n 階方程式

(C.7)

x(n) + p1 x(n−1) + · · · + pn−1 x′ + pn x = 0 (p1 , · · · , pn−1 , pn は定数)の場合まで一般化される。
例えば n = 1 の場合、すなわち

x′ + px = 0 (p は定数)

(C.8)

の場合は、特性方程式は λ + p = 0 であるから、特性根は λ = −p であり、x(t) = Ce−pt(C は任意定数).

は (C.8) の一般解である。
例えば n = 3 の場合、すなわち(C.9)

x′′′ + px′′ + qx′ + rx = 0 (p, q, r は定数)の場合、特性方程式は λ3 + pλ2 + qλ + r = 0. その根を α, β, γ とする。

63

(i) α, β, γ がすべて相異なるならば
x(t) = C1 eαt + C2 eβt + C3 eγt(C1 , C2 , C3 は任意定数)

は (C.9) の一般解である。

(ii) α = β 6= γ (二重根の場合) ならば
x(t) = C1 eαt + C2 teαt + C3 eγt(C1 , C2 , C3 は任意定数)

は (C.9) の一般解である。

(iii) α = β = γ (三重根の場合) ならば
x(t) = C1 eαt + C2 teαt + C3 t2 eαt(C1 , C2 , C3 は任意定数)

は (C.9) の一般解である。

(iv) α = a + ib, β = a − ib (a, b ∈ R, b 6= 0) ならば
x(t) = C1 eat cos bt + C2 eat sin bt + C3 eγt(C1 , C2 , C3 は任意定数)

は (C.9) の一般解である。

C.2

基本的な定理の証明(やはり定理の証明もつけておこう。とりあえず他の機会に書いた証明をこちらに持ってきた。しばらくは工事中。)

C.2.1

定義と例

定数 p, q と、区間 I 上で定義された関数 f が与えられたとき、

y ′′ + py ′ + qy = f (x)(C.10)

を定数係数 2 階線形常微分方程式とよぶ。特に f (x) ≡ 0 の場合の

y ′′ + py ′ + qy = 0(C.11)

を同次方程式、一般の場合の (C.10) を非同次方程式という。

C.2.2

特性方程式, 特性根

微分方程式 y ′′ + py ′ + qy = 0 に対して、2 次方程式

λ2 + pλ + q = 0
を微分方程式の特性方程式 (charasteristic equation), 特性方程式の根を特性根 (characteristicroot) とよぶ。
例 C.4 y ′′ − 3y ′ + 2y = 0 の特性方程式は λ2 − 3λ + 2 = 0 で、特性根は λ = 1, 2.

64





補題 C.5 α が λ2 + pλ + q = 0 の解ならば、y = eαx は (C.11) の解である。


証明



y=e

αx

ならば、

y ′ = αeαx ,y ′′ = α2 eαxであるから、y ′′ + py ′ + qy = (α2 + pα + q)eαx = 0 · eαx = 0.


補題 C.6 (重ね合せの原理 (principle of superposition)) y1 , y2 が (C.11) の解ならば、(C.12)

y = Ay1 + By2(A, B は任意定数)

も (C.11) の解である。


証明



y = Ay1 + By2 より
y ′ = Ay1′ + By2′ ,y ′′ = Ay1′′ + By2′′

であるから

y ′′ + py ′ + qy = A(y ′′ + py ′ + qy) + B(y ′′ + py ′ + qy) = A · 0 + B · 0 = 0.
上の二つの補題から、λ2 + pλ + q = 0 の 2 根を α, β とするとき、

y = Aeαx + Beβx

(A, B は任意定数)

は (C.11) の解であることが分かる。
例 C.7 y = Aex + Be2x (A, B は任意定数) は、y ′′ − 3y ′ + 2y = 0 の解である。
実は α 6= β の場合、(C.11) の解は (C.12) 以外にないことが示せる。一方、α = β の場合にはもう一工夫必要である。順番に考察していこう。

相異なる特性根を持つ場合

C.2.3




命題 C.8 (定数係数 2 階線形同次方程式 (1) 相異なる特性根を持つ場合) 特 性 方 程 式
λ2 + pλ + q = 0 が相異なる 2 根 α, β をもつとき、(C.13)

y = Aeαx + Beβx

(A, B は任意定数)

は (C.11) の一般解である。すなわち、(a) 任意の定数 A, B に対して、(C.13) で定まる y は (C.11) の解である。
(b) (C.11) の任意の解は、適当な定数 A, B を用いて、y = Aeαx + Beβx と一意的に表さ
れる。


証明



(a) は済んでいる。(b) については、次の二点を示せばよい。
65

[任意の解は eαx , eβx の 1 次結合で書ける] y が y ′′ + py ′ + qy = 0 の解だとする。
y ′ − βy
y1 :=
,
α−β

y ′ − αy
y2 :=
β−α

とおくと

y1 + y2 =11[(y ′ − βy) − (y ′ − αy)] =(α − β)y = y.
α−β
α−β

また

1
1
[(y ′′ − βy ′ ) − α(y ′ − βy)] =[y ′′ − (α + β)y ′ + αβy]
α−β
α−β
1
=
(y ′′ + py ′ + qy) = 0.
α−βy1′ − αy1 =これから、ある定数 C1 が存在して y1 = C1 eαx となることが分かる。同様にしてある定数 C2が存在して y2 = C2 eβx となることが分かる。ゆえにy = y1 + y2 = C1 eαx + C2 eβx .
[eαx , eβx の 1 次独立性] 定数 C1 , C2 について
C1 eαx + C2 eβx = 0(C.14)
が成り立ったとする。微分してC1 αeαx + C2 βeβx = 0.
この式から (C.14) の α 倍を引くと

C2 (β − α)eβx = 0.
仮定より α 6= β であるから C2 = 0. (C.14) に代入してC1 eαx = 0.
これから C1 = 0. C1 = C2 = 0 が示せたので、eαx , eβx は 1 次独立である。
例 C.9 y ′′ − 3y ′ + 2y = 0 の特性方程式は λ2 − 3λ + 2 = 0 で、特性根は λ = 1, 2. したがって

y = Aex + Be2x

(A, B は任意定数)が一般解である。

C.2.4

特性根が重根である場合



補題 C.10 λ2 + pλ + q = 0 が重根 α をもつとき、y = xeαx は (C.11) の解である。


証明



y = xeαxとすると、

y ′ = eαx + αxeαx ,y ′′ = αeαx + αeαx + α2 xeαx = α2 xeαx + 2αeαx66であるから、
y ′′ + py ′ + qy = (α2 x + 2α + pαx + p + qx)eαx = (α2 + pα + q)x + (p + 2α) eαxとなるが、α は λ2 + pλ + q = 0 の重根であるから、

α=

−p +

p

−p +
=
2

判別式

2

√

α2 + pα + q = 0,
0

=−

p
2

ゆえに

p + 2α = 0が成り立ち、y ′′ + py ′ + qy = [0 · x + 0] eαx = 0.



命題 C.11 (定数係数 2 階線形同次方程式 (2) 特性根が重根の場合) 特性方程式 λ2 +pλ+q = 0 が重根 α をもつとき、(C.15)

y = Aeαx + Bxeαx(A, B は任意定数)

は (C.11) の一般解である。すなわち、(a) 任意の定数 A, B に対して、(C.15) で定まる y は (C.11) の解である。
(b) (C.11) の任意の解は、適当な定数 A, B を用いて、y = Aeαx + Bxeαx と一意的に表
される。




証明

(a) は済んでいる。(b) については、次の二点を示せばよい。

[任意の解は eαx , xeαx の 1 次結合で書ける] y が y ′′ + py ′ + qy = 0 の解だとする。y = eαx uとおくと、
y ′ = αeαx + eαx u′ = eαx (u′ + αu),y ′′ = α2 eαx + αeαx u′ + eαx u′′ = eαx (u′′ + αu′ + α2 u)であるから

y ′′ + py ′ + qy = y ′′ − 2αy ′ + α2 y = eαx (u′′ + αu′ + α2 u) − 2αeαx (u′ + αu) + α2 eαx u= eαx u′′ .
y ′′ + py ′ + qy = 0 であったから eαx u′′ = 0. これから u′′ = 0. ゆえに定数 C1 , C2 が存在して
u = C1 + C2 x. ゆえに
y = eαx u = C1 eαx + C2 xeαx .
[eαx , xeαx の 1 次独立性] 定数 C1 , C2 について
C1 eαx + C2 xeαx = 0(C.16)
が成り立ったとする。微分してC1 αeαx + C2 eαx + C2 αxeαx = 0.
この式から (C.16) の α 倍を引くと

C2 eαx = 0.
ゆえに C2 = 0. (C.16) に代入してC1 eαx = 0.
これから C1 = 0. C1 = C2 = 0 が示せたので、eαx , xeαx は 1 次独立である。

67

例 C.12 y ′′ − 2y ′ + y = 0 の特性方程式は λ2 − 2λ + 1 = 0 で、特性根は 1 (重根). ゆえに一
般解は
y = Aex + Bxex (A, B は任意定数).

C.2.5

特性根が虚数である場合

微分方程式(C.17)y ′′ + y = 0の特性方程式は λ2 + 1 = 0 で、特性根は λ = ±i (i は虚数単位) である。そこで定理 C.8 を機械的に適用すると、一般解は(C.18)

y = Aeix + Be−ix(A, B は任意定数)となるが、eix , e−ix は一体何であろうか？
実は指数関数 ex は、複素変数に一般化され、その一般化された指数関数に対しても

d λx
e = λeλxdx(λ は複素数の定数)などの性質は保たれるので、(C.18) は確かに微分方程式 (C.17) の解を与えるのである。

68

要約: 複素変数の指数関数



指数関数は複素変数まで拡張できる。その定義には色々な方法があるが、どれを採用しても結果は一致する。ここでは z = x + iy (x, y は実数, i は虚数単位) に対して

ez = ex+iy := ex (cos y + i sin y)と定義する。

• 実変数に関する指数関数の拡張になっている。
• 指数法則 ez+w = ez ew が成立する。
• 特に Euler の公式
eiy = cos y + i sin y(C.19)

が成り立つ。y = π とすると eiπ = cos π + i sin π = −1 + i · 0 = −1 より有名な

eiπ + 1 = 0が得られる。(C.19) で y の代りに −y とした

e−iy = cos(−y) + i sin(−y) = cos y − i sin y(C.20)と (C.19) を連立方程式とみて、

cos y =

1 iy
e + e−iy ,2sin y =

1 iy
e − e−iy2iを得る。

• λ が複素数であってもd λz 
= λeλz .
edz• 任意の複素数 z に対して
z

e =

∞
X
zn
n=0

n!

が成立する。

• |ex+iy | = ex .

p, q が実定数の場合、λ + pλ + q = 0 が虚根を持てば、それは互いに複素共役である。ゆえに
λ = a ± ib (a, b ∈ R; b 6= 0)
2

と書ける。

Ae(a+ib)x + Be(a−ib)x = Aeax (cos bx + i sin bx) + Beax (cos bx − i sin bx)= (A + B)eax cos bx + i(A − B)eax sin bx.
C1 = A + B, C2 = i(A − B) とおくと、(C.21)

y = C1 eax cos bx + C2 eax sin bx.
69

また A, B は任意定数であることから、C1 , C2 も任意定数である31 。C1 , C2 を実数の範囲のみで動かせば、(C.21) は任意の実数値関数の解を表す。

補足 特性根が虚数である場合を定理の形にまとめておかなかったが、そうしておくべきで
あった。




命題 C.13 p, q を実定数とする。微分方程式

y ′′ + py ′ + qy = 0(C.11)の特性方程式 λ2 + pλ + q = 0 が虚数解 λ = a ± ib (a, b ∈ R, b 6= 0) を持つとき、

y = Aeax cos(bx) + Beax sin(bx)(C.22)

は (C.11) の一般解である。すなわち(a) 任意の定数 A, B に対して、y = Aeax cos(bx) + Beax sin(bx) で定まる y は (C.11) の解である。
(b) (C.11) の任意の解は、適当な定数 A, B を用いて、y = Aeax cos(bx) + Beax sin(bx) と一意的に表される。


この文書では (元々の対象者が大学 1 年生であることから)、証明はなるべく初等的にする方針で説明しているが、ここではまず線形空間の議論に慣れている人向けの証明をして、それから初等的な (だがやや面倒な) 証明を示すことにする。
n 階線形同次微分方程式に対し、

y = C1 φ1 (x) + · · · + Cn φn (x) (C1 , · · · , Cn は任意定数)が一般解であるとは、線形代数の言葉を用いると、φ1 (x), · · · , φn (x) が微分方程式の解空間(解全体の集合のなす線形空間) の基底であることである。言い換えると解全体の集合 = Span {φ1 (x), · · · , φn (x)}が成り立ち、φ1 , · · · , φn が 1 次独立であることである。
証明 (線形空間の議論に慣れている人向け)

a + ib と a − ib は相異なる特性根であるから

y = Ae(a+ib)x + Be(a−ib)x(A, B は任意定数)

は (C.11) の一般解である (命題 C.8 による)。これは

(C.11) の解空間 = Span e(a+ib)x , e(a−ib)xであり、e(a+ib)x , e(a−ib)x が 1 次独立であることを意味する。
一方
Ae(a+ib)x + Be(a−ib)x = (A + B)eax cos(bx) + i(A − B)eax sin(bx)


C1
これは
C
1
  2
i
C1
−2
2
.
i
1
C2
2
2
31




=

1
i

1
−i



A
B



 
A
=
において係数行列が正則であることから分かる。実際、
B70であるからSpan e(a+ib)x , e(a−ib)x ⊂ Span {eax cos(bx), eax sin(bx)} .

上に述べたことから左辺は 2 次元であり、右辺は 2 次元以下であるので、実は左辺と右辺は等しく、右辺は 2 次元すなわち eax cos(bx), eax sin(bx) は 1 次独立である。

(C.11) の解空間 = Span {eax cos(bx), eax sin(bx)}であるので、

y = Aeax cos(bx) + Beax sin(bx) (A, B は任意定数)
は (C.11) の一般解である。 (証明終)初等的な証明

まず、任意の定数 A, B に対して

y = Aeax cos(bx) + Beax sin(bx)で定義された y が、微分方程式 (C.11) の解であることを示そう。

y ′′ + py ′ + qy =eax A(a2 − b2 + ap + q) + B(bp + 2ab) cos(bx)+ eax B(a2 − b2 + ap + q) − A(bp + 2ab) sin(bx)λ = a ± ib が λ2 + pλ + q = 0 を満たすことから
a2 − b2 + ap + q = 0 ∧ 2ab + bp = 0が得られる。ゆえに y ′′ + py ′ + qy = 0 が成り立つ。

[任意の解は eax cos(bx), eax sin(bx) の 1 次結合で書ける] 命題 C.8 より、(C.11) の任意の解y に対して、ある定数 A, B が存在して y = Ae(a+ib)x + Be(a−ib)x が成り立つ。
Ae(a+ib)x + Be(a−ib)x = (A + B)eax cos(bx) + i(A − B)eax sin(bx)であるから、y は eax cos(bx), eax sin(bx) の 1 次結合で書ける。

[eax cos(bx), eax sin(bx) の 1 独立性] 定数 A, B について
Aeax cos(bx) + Beax sin(bx) = 0
が成り立ったとする。x = 0 を代入して

A = 0.
また、微分してから x = 0 を代入してAa + Bb = 0.
仮定 b 6= 0 に注意すると A = B = 0 が示せる。ゆえに eax cos(bx), eax sin(bx) は 1 次独立である。(証明終)
線形代数を応用すると、計算があまり必要のない証明が得られ、見通しが良くなる。例えば、n 階線形同次方程式の解空間が n 次元線形空間であることが一般的に示される。それが分かれば、1 時独立な解が n 個求まれば、それが基底であることがすぐ結論でき、一般解が得られる。(補足終)√−1 ± 3i
′′
′
2
. ゆえに
例 C.14 y + y + 1 = 0 の特性方程式は λ + λ + 1 = 0 で、特性根は λ =
2
一般解は
√
√
3x
3x
−x/2
−x/2
y = Aecos+ Be
sin(A, B は任意定数).
2
2

71

C.3

まとめ

y ′′ + py ′ + qy = 0 (p, q は定数) の一般解は次のように求まる。
(i) λ2 + pλ + q = 0 が相異なる 2 根 α, β を持つならば、
y = Aeαx + Beβx

(A, B は任意定数).

(ii) λ2 + pλ + q = 0 が重根 α を持つならば、
y = Aeαx + Bxeαx(A, B は任意定数).

(iii) p, q が実定数で、λ2 + pλ + q = 0 が虚根 a + ib (a, b ∈ R; b 6= 0) を持つならば、
y = C1 eax cos bx + C2 eax sin bx (C1 , C2 は任意定数).
このように 2 階方程式の一般解は、適当な二つの関数 φ1 , φ2 を用いて、

y = C1 φ1 (x) + C2 φ2 (x) (C1 , C2 は任意定数)と書ける場合がしばしばあるが、このとき φ1 , φ2 を基本解 (fundamental solution) とよぶ。
つまり

1. 特性方程式が相異なる 2 根 α, β を持つ場合、eαx , eβx は基本解
2. 特性方程式が相異なる重根 α を持つ場合、eαx , xeαx は基本解
3. 特性方程式が互いに複素共役である虚根 a ± ib を持つ場合、eax cos bx, eax sin bx は基本解

D

初期値問題の基礎理論 (かけ足で説明)

数学的議論が好きな人は、一度時間をかけて学ぶことを勧めたい。ここでは常識的なことをかけ足で説明する。

D.1

はじめに

常微分方程式という方法論を使い始めてみてすぐに気がつくのは、問題の解を既知の関数で表すことができない場合がある、という事実である。このことは常微分方程式のもっとも簡単な場合と言える「原始関数を求める問題」ですでにそうであったことから、容易に了解できるであろう。
この困難を解決する一つの方法は、特定の問題の解を表現できるような新しい関数を導入することである。これは一定の成果をおさめたが、一方で解が存在すること自体を保証する方
法32 が求められたのも当然のことであろう。

32

例えば、高校数学の段階でも、方程式の解を具体的に式で表示できなくても、グラフによる考察 (曲線と曲線が交わる — 厳密には中間値の定理を根拠とする) から解の存在の証明せよ、あるいは解の個数を調べよ、という問題があった。

72

D.2

解の存在

常微分方程式の初期値問題の解の存在を保証するには、次の定理が便利である (証明は省略する。例えば高橋 [1] の定理 IV.A.4 や、坂井 [39] の定理 1.7、あるいはコディントン・レヴィ
ンソン [30] を見よ。)。


定理 D.1 (f が連続ならば解は存在する, Peano の定理) f (t, x) が連続関数ならば、

x′ = f (t, x),x(t0 ) = x0に対して、十分小さな範囲では解が存在する。つまり、ある正数 δ と、区間 (t0 − δ, t0 + δ)で定義された関数 u で、

u′ (t) = f (t, u(t)) (t ∈ (t0 − δ, t0 + δ)),u(t0 ) = x0

を満たすものが存在する。


注意 D.2 しばしば「この微分方程式は現実に起っている現象を表しているのだから、解が存在するのは当たり前である」と考える人がいるが、それは誤解である。例えば物理現象を「記述する」微分方程式であるとしても、微分方程式は現象そのものではなく、せいぜい近似としか言えないものである。実際に現象が起こることはせいぜい「状況証拠」であって、数学的な証明にはならない。
余談 D.3 Peano の定理を、Cauchy-Peano の定理、Perron の定理と呼ぶ人もいる。一方、微分方程式の初期値問題の解の存在については、x についての Lipschitz 条件を仮定する定理 (Picardの定理と呼ばれる) も有名である (縮小写像の不動点定理を用いて証明できる)。


定理 D.4 (Lipschitz 条件を満たす場合の解の存在と一意性, Picard の定理)
Rn+1 = R × Rn の閉領域

E = Eτ,r,ξ,ρ := {(t, x) ∈ R × Rn | |t − τ | ≤ r, kx − ξk ≤ ρ}
で定義された連続関数 f : E → Rn が、x についての Lipschitz 条件(∃L > 0)(∀(t, x), (t, y) ∈ E) kf (t, x) − f (t, y)k ≤ L kx − ykを満たすならば、常微分方程式の初期値問題

dx= f (t, x),
dtx(τ ) = ξの解が|t − τ | ≤ r′ := min {r, ρ/M } ,

M :=

max
(t,x)∈Eτ,r,ξ,ρ

kf (t, x)kの範囲で存在し、一意的である。




定まった名前がない、と言う印象がある。高橋 [1] では、証明の方法を「ピカールの逐次近似法」と呼んでいるが、定理自身は「コーシーの定理」と呼んでいる。

73

証明はそのうちここに載せたいと考えているが、いつになるか分からないので、高橋 [1] を推奨しておく。

D.3

解の存在範囲 (爆発, 極大延長解)

前項の定理で、δ という範囲を制限するものが出て来てしまったが、これは仕方がないことである。
例 D.5 常微分方程式

dx
= x2
dt
は変数分離形なので容易に解ける。これにx(0) = a (a は正定数)という初期条件をつけた初期値問題の解はx(t) =

1
1/a − t(t ∈ (−∞, 1/a)).

lim x(t) = +∞ となっている。このように解の大きさが無限大に発散することを解の爆発、
t→1/a−0

そのときの変数の値 (ここでは 1/a) のことを爆発時刻と呼ぶ。t ≥ 1/a の範囲まで解を延長することは不可能であることに注意しよう。
注: 実は有名な logistic 方程式でも解の爆発はおこる (本文中に解の存在範囲が限定されることを説明した)。実際の現象とは関わりがない場合 (初期値が負であったり、環境収容力より大きかったり) なので、普通は問題とされないが。
爆発が起らないための十分条件としては、次のリプシッツ条件が有名である33 。


x に関する Lipschitz 条件
定数 L が存在して|f (t, x1 ) − f (t, x2 )| ≤ L|x1 − x2 |.



特に f (t, x) が x の 1 次関数 (f (t, x) = A(t)x + b(t) という形) である場合 (言い換えると線形微分方程式の場合) は、ごく緩い仮定 (例えば A(t) が有界) の下でリプシッツ条件が成り立つ(したがって解の爆発は起らない)。
爆発が起らない場合は、解は方程式が意味を持つ (f の定義域をはみ出ない) 範囲で存在することが知られている。これについて、以下なるべく簡潔な説明を試みる。

D を R × Rn の領域または閉領域とする34 。f : D → Rn は連続で、微分方程式(D.1a)dx= f (t, x)
dt

に初期条件

(D.1b)x(t0 ) = x0を課した初期値問題について考える。
33

Rudolf Otto Sigismund Lipschitz (1832–1903, Königsberg に生まれ、Bonn にて没する) にちなむ。
領域とは連結開集合のこと。また閉領域とは、ある領域の閉包になっている集合のこと (D が閉領域とは、ある領域 Ω が存在して、D = Ω が成り立つこと) 。
3474以下、初期値問題 (D.1a), (F.2) の解とは、その定義域が t0 を左端とする区間 I で、初期条件 x(t0 ) = x0 と、I 全体で微分方程式を満たすような関数 x(t) のことをいう。また、I のことを解 x(t) の定義区間と呼ぶことにする。
初期時刻 t0 の十分近く (ある正の数 δ に対して、t ∈ [t0 − δ, t0 + δ] を満たす範囲) では解(局所解) が存在する、という定理があるわけだが、それだけでは満足できない。出来る限り広い範囲での解の存在を保証してほしい。
ここではまず素朴な検討をしてみる。新たな初期値問題

dy= f (t, x)
dt
に初期条件

y(t0 + δ) = x(t0 + δ)を考えることで、もう少し先まで ([t0 , t0 + δ + δ ′ ] まで) 解を延ばすことができる。つまり yが [t0 + δ − δ ′ , t0 + δ + δ ′ ] で定義されているとき、()
x(t) (t ∈ [t0 , t0 + δ])e(t) :=
xy(t) (t ∈ (t0 + δ, t0 + δ + δ ′ ])e は、初期値問題 (D.1a), (F.2) の解である。
定めた x
ちょっと考えると、これを続けることで、どこまでも解を延ばすことが出来そうだが、有限の限界 T があるかもしれない (t0 + δ + δ ′ + δ ′′ + · · · = T < +∞ — 正の数を足し続けても、限りなく増えるわけではない)。
(D.1a), (F.2) の解のうち、真に大きな定義区間が存在しないような (それ以上延ばせない)解のことを、(D.1a), (F.2) の極大延長解と呼ぶことにする。例えば
(D.2)
(D.3)

dx
= x2 ,
dtx(0) = a1(−∞ < t < 1/a) は極大延長解である。
1/a − t微分方程式の初期値問題の解の一意性が成り立つ場合 (例えば x について局所 Lipschitz 条
件が成り立つ) には、極大延長解を構成するという方針で、極大延長解の存在が証明できる (例えば高野 [40] の §5.3)。この場合は最大延長解と呼ぶ方が適切かもしれない。
簡単のため、以下の議論では解の一意性が成り立つことを仮定する。
いくつかのテキストで、次のように説明されている。

の場合、x(t) =「極大延長解 x(t) (t0 ≤ t < T あるいは t0 ≤ t ≤ T ) は、t → T のとき x(t) が有界でないか、
f の定義域 D の境界に近づく」
これを用いると、解 x(t) が有界であることを示せば、解が D の境界に近づくことが分かる。
これは満足すべき状況に思えるかもしれない。
しかし D の境界に近づくとは、正確にはどういう意味であろうか。正直なことを言うと、私にはよく分からない (t → T のとき、x(t) は有界であるが、振動し、x(t) は境界から離れたり近づいたりすることがありそうに思われる)。定義を明記してあるテキストを見た覚えがな
い (単に私が不勉強なだけかもしれないが…)。定義を書いていないということは、証明もきちんとは書かれていないことを意味する。

75

仕方がないので、私自身は、次のことを使っている。
極大延長解のグラフ上で f が有界ならば、xT := lim x(t) が存在して、(T, xT ) は D の境界
t→T∁点である (すなわち (∀ε > 0) B((T, xT ); ε) は D と D の両方と交わる)。
証明には、次の二つの定理を用いる。


定理 D.6 D ⊂ R × Rn , f : D → Rn は連続、φ : I = (α, β) → Rn は連続, graph φ ⊂ D,dφ(t) = f (t, φ(t)) (t ∈ I).
dt
さらに f が graph φ 上で有界ならば、lim φ(t) が存在する。
t→β

もちろん、f が D 全体で有界である場合も、同じ結論が成り立つ。






定理 D.7 x が (t0 , t1 ] で連続で、(t0 , t1 ) で微分可能で、(t1 , x1 ) ∈ D,dx(t) = f (t, x(t)) (t ∈ (t0 , t1 )),dt
x(t1 ) = x1を満たすならば、x は (t0 , t1 ] で微分方程式を満たす。すなわち x は t1 でも微分可能であり、
dx(t1 ) = f (t1 , x(t1 ))
dt
を満たす。




Z t

(この 2 つの定理は、関数が初期値問題の解であることと、積分方程式 x(t) = x0 +f (s, x(s)) dst0の解であることが同値であること (これは常識的) を用いれば証明できる。難しくないが証明は省略する。例えば桂田 [41] を見よ。これは、コディントン・レヴィンソン [30] を参考にして考えたことをメモしたものなので、直接 [30] を見る方が良いかもしれない。)
極大延長解の場合に、もしも (T, xT ) が D の内点であれば、初期条件 x(T ) = xT の初期値問
題を考えることで、解が T を超えて延長できることが導かれ、矛盾が生じる。ゆえに (T, xT )は D の境界点である。

D.4

解の一意性

常微分方程式の初期値問題の解の存在が分かったとして、つぎに気になるのは、解がただ一つに限るかということである。
初期値問題
x′ = |x|2/3 , x(0) = 0
を考えよう。すぐにわかる x ≡ 0 以外に、
x(t) =(t − C1 )30(t < C1 )x(t) =(t ≥ C1 ),0(t − C2 )3
 (t − C1 )3 (t < C1 )(t < C2 )0(C1 ≤ t ≤ C2 )x(t) =(t ≥ C2 ),
3(t − C2 ) (t > C2 )なども解である (ただし C1 は 0 以下の、C2 は 0 以上の、ともに任意の定数とする)。
そこで解の一意性を保証する条件が知りたくなるが、次の定理が非常に有名である。(しかし、応用上はその次の定理 D.9 が便利なので、もしも 1 つだけ覚えるならば、定理 D.9 を覚えることを勧める。)

76





補題 D.8 (リプシッツ条件をみたす場合の一意性) 連続関数 f (t, x) が x に関するリプ
シッツ条件|f (t, x1 ) − f (t, x2 )| ≤ L|x1 − x2 | (L は定数)を満たすとき、常微分方程式の初期値問題

dx= f (t, x) (t ∈ [a, b]),dtx(a) = x0の解 x = φ1 (t) (t ∈ [a, b1 ]), x = φ2 (t) (t ∈ [a, b2 ]) に対して、

φ1 (t) = φ2 (t) (t ∈ [a, b∗ ], b∗ := min{b1 , b2 })が成り立つ。
証明

ψ(t) := φ1 (t) − φ2 (t) (t ∈ [a, b∗ ]) とおく。
Z t
φj (t) = x0 +f (s, φj (s)) ds (j = 1, 2)
a

より

Z t
[f (s, φ1 (s)) − f (s, φ2 (s))] ds

ψ(t) =aゆえにZ t|ψ(t)| ≤|f (s, φ1 (s)) − f (s, φ2 (s))| ds
Z t
Z t
≤
L |φ1 (s) − φ2 (s)| ds = L|ψ(s)|ds.
aaa

ここで M := max |ψ(t)| とおくと、
t∈[a,b∗ ]|ψ(t)| ≤ LM (t − a),Z t(t − a)2|ψ(t)| ≤ L
LM (s − a) ds = L2 M
,
2
a
Z t
(s − a)2(t − a)3|ψ(t)| ≤ L
L2 M
ds = L3 M
,···
2
3!
a
以下帰納的に[L(t − a)]n[L(b∗ − a)]n|ψ(t)| ≤ M≤M(n ∈ N).
n!
n!
これから n → ∞ の極限を考えると、ψ(t) ≡ 0 が分かる。ゆえに φ1 (t) = φ2 (t) (t ∈ [a, b∗ ]).
常微分方程式の初期値問題の場合、一意性が成り立つというのは解が枝分かれをしないことであるから、一意性を保証するには Lipschitz 条件は局所的なもので十分であり (つまり L は全体で統一的に取れなくても構わない)、例えば f が C 1 級であればよいことが分かる。すな
わち次が成立する (応用上とても便利なのでぜひ覚えておこう)。

77





定理 D.9 (C 1 級ならば一意性が成立) f が C 1 級の関数であるとき、常微分方程式の初期値問題
dx= f (t, x) (t ∈ [a, b]), x(a) = x0
dtの解 x = φ1 (t) (t ∈ [a, b1 ]), x = φ2 (t) (t ∈ [a, b2 ]) に対して、

φ1 (t) = φ2 (t) (t ∈ [a, b∗ ], b∗ := min{b1 , b2 })が成り立つ。


上の例の f (t, x) = |x|では、x = 0 で f は微分可能でなく (従って C 級でもなく)、
Lipschitz 条件も 0 のところで崩れていることに注意しよう。
2/3

D.5

1

まとめ

1 階正規形常微分方程式の初期値問題
x′ = f (t, x),f (t0 ) = x0

については、(1) f が連続でありさえすれば、(t0 , x0 ) の十分近くで解は存在する (局所解が存在する)。
(2) f が C 1 級であれば、解は一意である。
(3) f が C 1 級であっても爆発という現象がありうる。x に関するリプシッツ条件|f (t, x1 ) − f (t, x2 )| ≤ L|x1 − x2 | (L は定数)が成り立てば爆発は起こらない。特に有界な係数を持つ線形方程式では爆発は起こらない。

E

初期値問題の解の存在と一意性の証明(しばらく工事中)
やはり証明は必要と考えるようになったので、以下とりあえずスケッチだけ書いておく。(そのうちにきちんと完成させたいが、授業をする訳ではないので、カタツムリの速さであろう。)

E.1

Lipschitz 条件を仮定する場合比較的多くのテキストに載っている。たとえば明治大学の学生の場合は、オンラインで高橋[1] が読める。
以下では、関数解析を学ぶと紹介されることの多い、縮小写像に関する不動点定理を用いる証明を紹介する。

78

E.1.1Banach 空間 C([a, b]; R)a, b ∈ R, a < b とするとき、I = [a, b] とおく。[a, b] 上の実数値連続関数の全体を C([a, b]; R)と表す。
f, g ∈ C([a, b]; R) について、和 f + g : I → R を(f + g)(x) := f (x) + g(x) (x ∈ I)
で、また f ∈ C([a, b]; R), c ∈ R について積 cf : I → R を(cf )(c) := cf (x) (x ∈ I)で定めると、f + g, cf ∈ C([a, b]; R) が成り立ち、C([a, b]; R) は R 上の線形空間になる。またf ∈ C([a, b]; R) に対して
kf k := max |f (x)|x∈[a,b]とおくと、k·k は C([a, b]; R) 上のノルムとなる。また C([a, b]; R) はそのノルムにより Banach
空間となる。Banach 空間は距離空間として完備である。

E.1.2

不動点, 縮小写像に関する不動点定理

完備距離空間については Banach の不動点定理という重要な定理がある。
(X, d) を距離空間とする。
写像 T : X 3 x 7→ T x ∈ X について、x ∈ X が T の不動点であるとは、

Tx = xが成り立つことをいう。
写像 T : X 3 x 7→ T x ∈ X が縮小写像であるとは、(∃k ∈ [0, 1))(∀x, y ∈ X) d(T x, T y) ≤ kd(x, y)が成り立つことをいう。



定理 E.1 (Banach の不動点定理, 縮小写像に関する不動点定理) (X, d) を完備距離空間、
T : X → X が縮小写像ならば、T の不動点 x が存在し、それは一意的である。また、xは以下の逐次近似法により得られる： 任意の x0 ∈ X に対して

xn+1 := T xn

(n ∈ N)により X 内の点列 {xn }n∈N を定めると、n → ∞ のとき {xn } は T の不動点に収束する。



証明

任意の x0 ∈ X を選び、xn+1 = T xn

(n = 0, 1, · · · )により {xn }n≥0 を定める。
Step 1 任意の n ∈ N ∪ {0} に対して

(E.1)

d(xn , xn+1 ) ≤ k n d(x0 , x1 )が成り立つ。実際、

79

(i) n = 0 のとき、(E.1) の左辺 = (E.1) の右辺であるから、(E.1) が成立する。
(ii) n = j のとき (E.1) 成り立つ、すなわち d(xj , xj+1 ) ≤ k j d(x0 , x1 ) が成り立つと仮定す
るとd(xj+1 , xj+2 ) = d(T xj , T xj+1 ) ≤ kd(xj , xj+1 ) ≤ k · k j d(x0 , x1 ) = k j+1 d(x0 , x1 ).
すなわち n = j + 1 のときも (E.1) が成り立つ。
数学的帰納法により、(E.1) は任意の n ≥ 0 について成り立つ。
Step 2 {xn }n≥0 は X における Cauchy 列である。実際、n, m ∈ N, n > m とするとき

d(xm , xn ) ≤ d(xm , xm+1 ) + d(xm+1 , xm+2 ) + · · · + d(xn−1 , xn )
=

n−m−1
X

d(xm+j , xm+j+1 )

j=0

≤
≤

n−m−1
X

k

m+j

m

d(x0 , x1 ) = k d(x0 , x1 )

j=0
m

n−m−1
X
j=0

1 − k n−mk = k d(x0 , x1 )
1−k
j

m

k
d(x0 , x1 ).
1−km > n のときはd(xm , xn ) ≤knd(x0 , x1 ).
1−k

一般に、任意の n, m ∈ N に対してd(xm , xn ) ≤k min{n,m}d(x0 , x1 )1−kが成り立つ。ゆえに {xn }n≥0 は Cauchy 列である。
Step 3 X は完備であるから、 ある x∞ ∈ X が存在して、lim xn = x∞ .

n→∞

これから

0 ≤ d(x∞ , T x∞ ) = lim d(xn , xn+1 ) = 0.
n→∞

ゆえに d(x∞ , T x∞ ) = 0. ゆえに T x∞ = x∞ . これは x∞ が T の不動点であることを意味する。
さて、T : X → X, n ∈ N ∪ {0} に対して、T n : X → X を以下のように帰納的に定める。

T 0 x = x (x ∈ X),
T i x = T (T i−1 x) (i ∈ N).
縮小写像の定理の系として得られる次の定理は非常に役立つ。



定理 E.2 (冪が縮小写像ならば不動点を持つ) X を K 上の Banach 空間、T : X → X とする。ある自然数 m に対して、T m が縮小写像ならば、T の不動点 x が存在し、それは一意的である。また、x は以下の逐次近似法により得られる： 任意の x0 ∈ X に対して

xn+1 := T xn

(n ∈ N)により X 内の点列 {xn }n∈N を定めると、n → ∞ のとき {xn } は T の不動点に収束する。



手品のような証明という印象がある。

80



証明

T m は縮小写像なので、T m の不動点 x∗ が存在する。すなわちT m (x∗ ) = x∗ .

このとき

T m (T x∗ ) = T m+1 x∗ = T (T m x∗ ) = T x∗
であるから、T x∗ は T m の不動点である。T m の不動点の一意性から

T x ∗ = x∗である。これは x∗ が T の不動点であることを表わす。T の不動点は T m の不動点でもあるから、T m の不動点の一意性から T の不動点の一意性が導かれる。
任意の x0 ∈ X を取って、
xn+1 = T xn
で点列 {xn }n∈N を定める。各 r ∈ {0, 1, · · · , m − 1} に対して、この点列の部分列 {xmj+r }j∈Nを考えると、これは xr を初期値として、漸化式

xm(j+1)+r = T m (xmj+r )で作った点列である。T m が縮小写像であることから、これは T m の唯一の不動点 x∗ に収束する:
lim xmj+r = x∗ .
j→∞

任意の r に対してこれが成り立つことから、

lim xn = x∗ .

n→∞

E.1.3

連続かつ Lipschitz 条件を満たす場合の解の存在と一意性

D := [α, α + a] × [β − b, β + b], f : D → R は連続であるとする。また、f は x について
Lipschitz 条件を満たすとする。
(∃L ∈ R)(∀(t, x) ∈ D)(∀(t, x′ ) ∈ D)

|f (t, x) − f (t, y)| ≤ L |x − y| .

このとき、次の常微分方程式の初期値問題を考える。

(E.2)

x′ (t) = f (t, x(t)) (t ∈ [α, α + a])(E.3)x(α) = βM := max |f (t, x)| ,t∈[α,α+a]bδ := min a,M,

X := {φ ∈ C([α, α + δ]) | |φ(t) − β| ≤ M (t − α).}とおくと、X は C([α, α + δ]; R) の閉集合である。ゆえに X は完備距離空間である。
写像 T : X → X を
Z t(T x)(t) := β +f (s, x(s)) ds (x ∈ X, t ∈ [α, α + a])α81で定義すると、

Z t|T x(t) − T y(t)| = |[f (s, x(s)) − f (s, y(s))] ds| ≤ L|t − α| max |x(s) − y(s)|,sα|T 2 x(t) − T 2 y(t)| ≤ L max |(T x)(s) − (T y)(s)||t − α| ≤s

(L|t − α|)2
max |x(s) − y(s)|,s2

以下帰納的に|T m x(t) − T m y(t)| ≤(L|t − α|)m
(Lr′ )mmax |x(s) − y(s)|. ≤max |x(s) − y(s)|.
s
s
m!
m!

ゆえに

kT m x − T m yk ≤(Lr′ )mkx − yk .
m!

十分大きな m に対して、(Lr′ )m<1m!
となり、T m : X → X は縮小写像になっている。ゆえに定理 E.2 により、T は X で一意的な不動点 x を持つ:
T x = x.
すなわち

Z t
f (s, x(s)) ds (t ∈ [α, α + δ]).

x(t) = β +αE.2Lipshitz 条件を仮定しない場合

(準備中)

F

線形微分方程式の常識

F.1

取り扱う方程式

I は R の区間、A : I → Rn×n , f : I → Rn は連続とするとき、1 階線形微分方程式(F.1)dx(t) = A(t)x(t) + f (t) (t ∈ I)
dt

を考える。
任意の t0 ∈ I, ξ = (a1 , · · · , an )⊤ ∈ Rn に対して、(F.1) と(F.2)

x(t0 ) = ξを満たす x = x(t) を求めよ、という問題を初期値問題という。a を初期値、(F.2) を初期条
件と呼ぶ。

82

F.2

高階微分方程式の 1 階微分方程式への帰着

1 次元 n 階線形微分方程式
x(n) (t) + a1 (t)x(n−1) (t) + · · · + an−1 (t)x′ (t) + an (t)x(t) = f (t) (t ∈ I)

(F.3)
と初期条件

(F.4)

x(t0 ) = ξ1 ,x′ (t0 ) = ξ2 ,· · · x(n−1) (t0 ) = ξnを考えよう。
ここで aj (j = 1, · · · , n) と f は I で連続な実数値関数、ξ = (ξ1 , · · · , ξn ) ∈ Rn とする。


x(t) := 




x(t)
x′ (t)
..
.



,


x(n−1) (t)A(t) := 




0
0
..
.

1
0
..
.



0
1
...

0
0
−an (t) −an−1 (t) · · ·

...
...

0
0
1
−a2 (t) −a1 (t)




,





0
 .. 
 . f (t) := 

 0 
f (t)


とおくと、(F.3) は n 次元の 1 階線形微分方程式 (F.1) に、また (F.4) は (F.2) に帰着される。

F.3

解の存在と一意性

任意の初期値 ξ に対して、初期値問題 (F.1), (F.2) の一意的な解 x ∈ C 1 (I) が存在する。
実際
F (t, x) := A(t)x + f (t) (t ∈ I, x ∈ Rn )とするとき、I に含まれる任意のコンパクト区間 K に対して

kF (t, x1 ) − F (t, x2 )k ≤ L kx1 − x2 k ,L := max kA(t)kt∈K

ゆえに F |K×Rn は、連続かつ x について Lipschitz 条件を満たす。ゆえに t ∈ K で解が存在し、一意的である。K の任意性から、t ∈ I で解が存在し、一意的である。

F.4

同次方程式の解空間は線形空間

特に f = 0 の場合、

dx(t) = A(t)x(t) (t ∈ I)
dt

(F.5)

を線形同次方程式と呼ぶ。
この微分方程式の解全体の集合 (以下では解空間と呼ぶ)

V := {φ | φ′ (t) = A(t)φ(t) (t ∈ I)}
は、通常の関数の和・スカラー積の定義のもとで線形空間をなす。実際 φ1 , φ2 ∈ V とすると
き、任意の c1 , c2 ∈ R に対して φ := c1 φ1 + c2 φ2 とおくと

d(c1 φ1 (t) + c2 φ2 (t)) = c1 φ′1 (t) + c2 φ′2 (t) = c1 A(t)φ1 (t) + c2 A(t)φ2 (t)dt= A(t) (c1 φ1 (t) + c2 φ2 (t)) = A(t)φ(t)φ′ (t) =であるから、φ ∈ X.

83

注意 F.1 (線形空間の定義, 線形部分空間) 集合 X と、写像 X × X 3 (f, g) 7→ f + g ∈ X,
R × X 3 (c, f ) 7→ cf ∈ X が定義されていて、以下の条件 (1)–(8) を満たすとき、X は R 上の線形空間 (linear space, ベクトル空間, vector space) という。

(L1) (f + g) + h = f + (g + h)(L2) (∃0 ∈ X) (∀f ∈ X) f + 0 = 0 + f = f(この条件を満たす 0 は一意的に定まる。)(L3) (∀f ∈ X)(∃g ∈ X) f + g = g + f = 0(f に対して g は一意的に定まり、普通 −f と表される。)(L4) (∀f, g ∈ X) f + g = g + f(L5) (∀c ∈ R) (∀f, g ∈ X) c(f + g) = cf + cg
(L6) (∀c1 , c2 ∈ R) (∀f ∈ X) (c1 + c2 )f = c1 f + c2 f
(L7) (∀c1 , c2 ∈ R) (∀f ∈ X) (c1 c2 )f = c1 (c2 f )(L8) (∀f ∈ X) 1f = f
一般に、空でない集合 Ω が与えられたとき、Ω 上の実数値関数全体の集合 {f | f : Ω → R}
を X とおき、f, g ∈ X に対して f + g ∈ X を(f + g)(x) = f (x) + g(x) (x ∈ Ω)
で定め、c ∈ R, f ∈ X に対して cf ∈ X を(cf )(x) = cf (x) (x ∈ Ω)で定めると、X は R 上の線形空間となる。
線形空間 X の部分集合 V が (加法・スカラー乗法は X のそれを使って)(F.6a)

f, g ∈ V ⇒ f + g ∈ V,

(F.6b)c ∈ R, f ∈ V ⇒ cf ∈ V

を満たすとき、V は線形空間となる。この V を X の線形部分空間という。
この節では、線形微分方程式の解空間がしばしば登場するが、それらは上の X の線形部分空間であり、(F.6a), (F.6b) という 2 つの条件を確認するだけでよくなる (一々上の (L1)–(L8)を確認する必要はない)。

F.5

線形同次方程式の解の 1 次独立性, Wronskian による判定(F.5) の解空間の次元は n であることを証明することが大事な目標であるが、まず 1 次独立性 (線形独立性) の定義を想い出そう。
線形空間 X の要素 φ1 , · · · , φn が 1 次独立であるとは(F.7)(∀c1 , · · · , cn ∈ R)

[c1 φ1 + · · · + cn φn = 0 ⇒ c1 = · · · = cn = 0]が成り立つことをいう。

84

φ1 , · · · , φn が線形同次微分方程式 x′ (t) = A(t)x(t) (t ∈ I) の解である場合は、
c 1 φ1 + · · · + c n φn = 0
は(∀t ∈ I) c1 φ1 (t) + · · · + cn φn (t) = 0(F.8)

を意味する。このことは容易に理解できるが、この条件 (無限個の t について等式が成り立つ)を確認するのは、一見してそれほど簡単でないように思われるかもしれない。実は、この後紹介する定理 F.3 によると、線形同次微分方程式 x′ (t) = A(t)x(t) の n 個の解 φ1 , · · · , φn に対しては、任意に選んだ 1 つの t0 ∈ I に対して

c1 φ1 (t0 ) + · · · + cn φn (t0 ) = 0 ⇒ c1 = · · · = cn = 0(φ1 , · · · , φ1 が Rn のベクトルとして 1 次独立)が成り立つことを確認すれば良い。


定義 F.2 (Wronskian) I は R の区間、A : I → Rn×n とする。線形同次微分方程式 x′ (t) =A(t)x(t) の n 個の解 φ1 , · · · , φn に対して

Φ(t) := (φ1 (t) φ2 (t) · · · φn (t)) ,W (t) = W [φ1 , · · · , φn ](t) := det Φ(t)ロ ン ス キ ーとおき、この W を、微分方程式 x′ (t) = A(t)x(t) の解 φ1 , · · · , φn のWronski 行列式あるいは Wronskian と呼ぶ。







定理 F.3 (Wronskian の満たす微分方程式) I は R の区間、A : I → Rn×n とする。線形
同次微分方程式 x′ (t) = A(t)x(t) の n 個の解 φ1 , · · · , φn に対して、W (t) を φ1 , · · · , φnの Wronskian とすると

W ′ (t) = tr A(t)W (t)

(F.9)
def.

が成り立つ (tr(aij ) =nXajj )。ゆえに次式が成り立つ。

j=1



Z t
(F.10)
tr A(s) ds

W (t) = W (t0 ) exp
t0(t ∈ I).
証明(F.9) を証明しよう。行列 Φ(t) の第 (i, j) 成分を φji (t) と表すことにする (i 番目の解nX′であるベクトル値関数 φi の第 j 成分である— φji (t) =ajk φki (1 ≤ j ≤ n) が成り立つこ
k=1

85

とに注意する)。

φ11 φ12 · · ·d φ21 φ22 · · ·W ′ (t) =..
dt ...
.
φ11 φ12 · · ·

φ1n
φ2n
..
.
φ1n

φ11
φ′11 φ′12 · · · φ′1nφ′21
φ21 φ22 · · · φ2n
+
= .
..
..
..
..
.
.
.
φn1
φn1 φn2 · · · φnnnnXXa1k φk1
a1k φk2 · · ·
k=1

=

k=1

φ21
..
.

φ22
..
.

···

φn1

φn2

···

+ ··· +

n
X

φ11
φ21
..
.

k=1

+ ··· +

φ11
φ21
..
.

φn2 · · ·nX
a1k φkn

φ11 φ12 · · ·φ1n′φ21 φ22 · · ·φ2n..
.. + · · · ..
.
.
.
′φ1n φ′1n · · ·
φnn
n
X

k=1

φ2n
..
.

+ k=1

φnn

n
X

ank φk1

a11 φ11 a11 φ12 ·
φ21
φ22
···
=
..
..
.
.
φn1
φn2 · · ·φ12 · · ·φ′22 · · ·
..
.

φ11

ank φk2 · · ·

k=1

n
X

φ12

φ′nn
···

a2k φk2 · · ·

..
.
φn2

n
X
k=1

···

φ1n
a2k φkn

..
.
φnn

φ1n
φ2n
..
.
ank φkn

k=1

φ11
φ12
···
a11 φ1na22 φ21 a22 φ22 · · ·
φ2n
+
..
..
..
.
.
.
φn1
φn2 · · ·
φnn
φ12
φ22
..
.

k=1

..
.
φn1

···
···

φ12
φ22
..
.

a2k φk1

n
X

φ1n
φ2n
..
.

···
···

φ1n
a22 φ2n
..
.
φnn

φ1n
φ2n
..
.

ann φn1 ann φn2 · · ·ann φnn

φ11 φ12 · · ·φ21 φ22 · · ·= (a11 + a22 + · · · + ann ) .
..
..
.
φn1 φn2 · · ·

φ1n
φ2n
..
.
φnn

= tr A(t) · W (t).
(F.9) を解けば (F.10) が得られる。
この定理から W (t) について、(a) 任意の t に対して W (t) = 0 であるか(b) 任意の t に対して W (t) 6= 0 であるかどちらか一方のみが成り立つが成り立つことが分かる。ゆえに、ある t0 ∈ I に対して、φ1 (t0 ),· · · , φn (t0 ) が Rn の 1 次独立なベクトルならば、W (t0 ) 6= 0.
86

(a) が成り立たないので、(b) が成り立つ。すなわち任意の t ∈ I に対して
W (t) 6= 0.
これは φ1 (t), · · · , φn (t) が Rn の 1 次独立なベクトルであることを意味する。

F.5.1

ロンスキアンを用いない方法(しばらく工事中)
ロンスキアンを用いた議論は分かりやすいが、ロンスキアンが登場しないテキストも結構ある。実は解の一意性定理だけで十分であるようだ。
定理を述べて証明する前に、二、三注意しておく。
φ1 , · · · , φn が 1 次独立という条件は
c1 φ1 + · · · + cn φn = 0 ⇒ c1 = · · · = cn = 0の形で述べられることが多いが、以下では対偶(c1 , · · · , cn ) 6= 0 ⇒nXci φi 6= 0i=1の形で用いる。
関数として 1 次独立・1 次従属という場合、= 0 と 6= 0 は定数関数 0 に等しい、等しくないということであるから、nX

ci φi = 0 ⇔ (∀t ∈ I)

n
X

i=1

i=1

n
X

n
X

ci φi 6= 0 ⇔ (∃t0 ∈ I)i=1

ci φi (t) = 0,ci φi (t0 ) 6= 0

i=1

である。




定理 F.4 A : I → Rn×n は連続、φ1 , · · · , φn は微分方程式 x′ (t) = A(t)x(t) の解とする。
このとき次の 3 つの条件は互いに同値である。

(i) (∃t0 ∈ I) {φi (t0 )} は (Kn のベクトルとして) 1 次独立nXnci φi (t0 ) 6= 0 ということである。
この条件は (∃t0 ∈ I)(∀c ∈ K \ {0})i=1(ii) (∀t ∈ I) {φi (t)} は (Kn のベクトルとして) 1 次独立nXnci φi (t) 6= 0 ということである。
この条件は (∀t ∈ I)(∀c ∈ K \ {0})i=1(iii) (I 上の Kn 値関数として) {φi } は 1 次独立nXnci φi 6= 0 (I 上の Kn 値関数として) ということであ
この条件は (∀c ∈ K \ {0})i=1り、結局 (∀c ∈ K \ {0})(∃t0 ∈ I)nnX

ci φi (t0 ) 6= 0 と同値である。

i=1





87

この定理の中の「1 次独立」を「1 次従属」に置き換えた次の命題が成立する (実際、次の命
題の (i), (ii), (iii) の条件は、それぞれ上の定理の (ii), (i), (iii) の否定であるから、上の定理から、互いに同値である。)。



定理 F.5 A : I → Rn×n は連続、φ1 , · · · , φn は微分方程式 x′ (t) = A(t)x(t) の解とする。
このとき次の 3 つの条件は互いに同値である。

(i) (∃t0 ∈ I) {φi (t0 )} は (Kn のベクトルとして) 1 次従属nXn
この条件は (∃t0 ∈ I)(∃c ∈ K \ {0})ci φi (t0 ) = 0 ということである。
i=1(ii) (∀t ∈ I) {φi (t)} は (Kn のベクトルとして) 1 次従属nXn
この条件は (∀t ∈ I)(∃c ∈ K \ {0})ci φi (t) = 0 ということである。
i=1(iii) (I 上の Kn 値関数として) {φi } は 1 次従属nXn
この条件は (∃c ∈ K \ {0})ci φi = 0 (I 上の Kn 値関数として) ということであi=1

り、結局 (∃c ∈ K \ {0})(∀t0 ∈ I)nnX

ci φi (t) = 0 と同値である。

i=1





では定理 F.4 を証明しよう。
(ちなみに Wronskian についての等式が分かっていると、以下の証明の (i) ⇒ (ii) と (iii) ⇒(ii) の証明で、解の一意性定理を使わないですむ。)
証明(ii) ⇒ (i) の証明: I 6= ∅ であるから、∀t ∈ I は ∃t0 ∈ I を導く。
(i) ⇒ (ii) の証明: (i) を仮定すると、任意の c ∈ K \ {0} に対して
n

期値問題の解の一意性から、任意の t ∈ I に対して
n
X

n
X

n
X

ci φi (t0 ) 6= 0. 初i=1ci φi (t) 6= 0. (∵ ある t1 ∈ I で

i=1

ci φi (t1 ) = 0 となるならば、t = t1 を初期時刻とした初期値問題を考えて、任意の t ∈ I に

i=1

対して

n
X

ci φi (t) = 0 が導かれて矛盾が生じる。) ゆえに

i=1(∀t ∈ I)(∀c ∈ K \ {0})nnXci φi (t) 6= 0

i=1が成り立つ (2 つの ∀ の順番を入れ替えただけ)。すなわち (ii) が成り立つ。

(iii) ⇒ (ii) の証明: (iii) を仮定すると、任意の c ∈ Kn \ {0} に対して

n
X
i=1

関数としての式 (定数関数 0 ではない) であるから、ある t0 ∈ I が存在して

ci φi 6= 0. これは
n
X
i=1

88

ci φi (t0 ) 6= 0.

nX

初期値問題の解の一意性から、任意の t1 ∈ I に対して

ci φi (t1 ) 6= 0. ゆえに

i=1(∀t ∈ I)(∀c ∈ K \ {0})nnXci φi (t) 6= 0

i=1が成り立つ。すなわち (ii) が成り立つ。

(ii)⇒(iii) の証明: (ii) が成り立つとする。
nX(∀t ∈ I)(∀c ∈ K \ {0})nci φi (t) 6= 0.

i=1

連続する ∀ の順序は交換可能であるから
nX(∀c ∈ K \ {0})(∀t ∈ I)nci φi (t) 6= 0.

i=1

I 6= ∅ であるから

nX(∀c ∈ K \ {0})(∃t0 ∈ I)nci φi (t0 ) 6= 0.

i=1

すなわち (iii) が成り立つ。
注意 F.6 (しばらく工事中) 1 次独立性は、ふつうnX(∀c ∈ K )n!

ci φi = 0 ⇒ c = 0

i=1と書かれるが、上では(∀c ∈ K \ {0})nnXci φi 6= 0

i=1として議論した。個人的にその方が分かりやすいと思ったのだが、そうする必要はない、という指摘があった。
!
nX(i) (∃t0 ∈ I) (∀c ∈ Kn )ci φi (t0 ) = 0 ⇒ c = 0i=1nX(ii) (∀t ∈ I) (∀c ∈ Kn )ci φi (t) = 0 ⇒ c = 0i=1(iii) (∀c ∈ Kn )(∀t ∈ I)nX!

!
ci φi (t) = 0⇒ c = 0.

i=1(ii)⇒(i)⇒(iii) は明らかであり (と言われた。確かにそうかもしれない。)、(iii)⇒(ii) のみ初期値nXn
問題の解の一意性を使えば良い。実際、任意の t0 ∈ I, c ∈ K が与えられたとき、 ci φi (t0 ) =i=1nX
0 が成り立つとすると、 ci φi (t) は微分方程式の解であるから、(t0 を初期時刻と考えて) 初i=1

期値問題の解の一意性により任意の t ∈ I に対して

n
X
i=1

c = 0.
89ci φi (t) = 0. ここで仮定 (iii) を使うと

F.6

線形同次方程式の解空間の次元

定理 F.7 (線形同次方程式の解空間の次元) n ∈ N, I は R の区間、A : I → Rn×n は連続であるとき、
V = {φ | φ′ (t) = A(t)φ(t) (t ∈ I)}の次元は n である。





証明 1 ≤ i ≤ n に対して、第 i 成分が 1 で、その他の成分が 0 である n 次元ベクトルを ei とする。e1 , · · · , en は 1 次独立である。そこで、任意の t0 ∈ I に対してx′ (t) = A(t)x(t) (t ∈ I),

x(t0 ) = eiの解を φi (i = 1, 2, · · · , n) とすると、φ1 , · · · , φn は 1 次独立である。
x が x′ (t) = A(t)x(t) (t ∈ I) の任意の解とする。このとき

a = (a1 , · · · , an )⊤ := x(t0 ),y(t) := a1 φ1 (t) + · · · + an φn (t)とおくと、

y(t0 ) = a1 φ1 (t0 ) + · · · + an φn (t0 ) = a1 e1 + · · · + an en = a = x(t0 ).
初期値問題の解の一意性によって y = x. すなわち

x(t) := a1 φ1 (t) + · · · + an φn (t) (t ∈ I).
ゆえに φ1 , · · · , φn は解空間の基底である。ゆえに解空間の次元は n である。

F.7

定数係数高階単独線形常微分方程式

p(λ) ∈ R[λ], deg p(λ) = n, D = dtd とするとき、常微分方程式
p(D)x = 0
を考えよう。
p(λ) の相異なる根を α1 , · · · , αr とすると

p(λ) = arY(λ − αj )

j=1

mj

,

r
X

mj = nj=1のように因数分解できる。このとき

x=

j −1
r m
X
X

cjk tk eαj t(cjk は任意定数)

j=1 k=0が一般解である。



補題 F.8 (1) (D − α) (eαt φ(t)) = eαt Dφ(t).

(2) 任意の m ∈ N に対して (D − α)m (eαt φ(t)) = eαt Dm φ(t).



証明

90



(1) 積の微分法により D (eαt φ(t)) = αeαt φ(t) + eαt φ′ (t) である。右辺第 1 項を左辺に移項して、(D − α)(eαt φ(t)) = eαt Dφ(t).
(2) 帰納法による。m = 1 のときは (1) そのものである。m = k のとき成り立つと仮定すると
(D − α)k+1 eαt φ(t) = (D − α) (D − α)k eαt φ(t) = (D − α) eαt Dk φ(t) = eαt DDk φ(t)= eαt Dk+1 φ(t).
すなわち m = k + 1 のときも成立する。ゆえに任意の m ∈ N について成り立つ。

((1) を D = e−αt (D − α) (eαt ·) と書くと、Dm = e−αt (D − α)m (eαt ·) となり (2) が分かりや
すいかも、とか色々考えるけれど、決定版が見つかっていない感がある。)



補題 F.9 (D − α)m x = 0 の一般解はx(t) =m−1
X

Cj tj eαt(C1 , . . . , Cm−1 は任意定数).

j=0



証明



y(t) := e−αtx(t) とおくと、
(D − α)m x(t) = (D − α)m eαt y(t) = eαt Dm y(t).

ゆえに(D − α) x(t) = 0 ⇔ D y(t) = 0 ⇔ (∃C0 , · · · , Cm−1 ) y(t) =
m

m

m−1
X

Cj tj

j=0

⇔ (∃C0 , · · · , Cm−1 ) x(t) =m−1
XCj tj eαt .

j=0

1 ≤ j ≤ r なる j に対して、
φjk (t) := tk eαj t(0 ≤ k ≤ mj − 1)

は (D − αj )mj x(t) = 0 の解である。

p(D) = aY!
(D − αk )mk(D − αj )mjk̸=jであるから、p(D)x = 0 の解でもある。φjk (1 ≤ j ≤ r, 0 ≤ k ≤ mj − 1) は全部で n 個ある
ので、それらの 1 次独立性を示せば、p(D)x = 0 の解空間の基底であることが分かる。
1 次独立性の証明は難しくはないが、一般的に記述しようとすると、(添字が増えて) 式が複雑になり、議論が分かりにくく感じる人が多くなると思われる。以下では、特別な場合に証明する (それが理解出来れば一般の場合に成り立つことが理解できると思われる)。
まず補題を準備する。



d
補題 F.10 (1) f (λ) ∈ C[λ], D = dtとするとき f (D)eαt = f (α)eαt .

(2) k, ℓ を 0 以上の任意の整数とするとき
  k! tk−ℓ eαt (0 ≤ ℓ ≤ k)ℓ k αt(D − α) t e =ℓ!
 0(k > ℓ).




証明

91

(1) Dk eαt = αj eαt を用いる。
(2)



k!
(D − α)ℓ tk eαt = eαt e−αt (D − α)ℓ eαt tk = eαt Dℓ tk = eαt tk−ℓ .
ℓ!

1 次独立性の証明のあらすじ α, β, γ が相異なるとき、eαt , eβt , teβt , eγt , teγt , t2 eγt が 1 次独立であること、すなわち
x(t) = A1 eαt + B1 eβt + B2 teβt + C1 eγt + C2 teγt + C3 t2 eγtについて

x(t) = 0 ⇒ A1 = B1 = B2 = C1 = C2 = C3 = 0を示そう。
C1 = C2 = C3 = 0 を示す (A1 = 0 や B1 = B2 = 0 の証明も同様である)。

pℓ (λ) := (λ − α)(λ − β)2 (λ − γ)ℓ(ℓ = 0, 1, · · · )とおくとき、(D − α)eαt = 0, (D − β)2 eβt = 0, (D − β)2 teβt = 0 であるので、任意の ℓ に対して
pℓ (D)eαt = 0, pℓ (D)eβt = 0, pℓ (D)teβt = 0.
ゆえにpℓ (D)x = pℓ (D) C1 eγt + C2 teγt + C3 t2 eγt .

上で紹介した補題を用いるとp2 (D)t2 eγt = (D−α)(D−β)2 (D−γ)2 t2 eγt = (D−α)(D−β)2

2! 2−2 γtt e = (D−α)(D−β)2 eγt = (γ−α)(γ−β)22!

p2 (D)teγt = (D − α)(D − β)2 (D − γ)2 teγt = (D − α)(D − β)2 0 = 0,
p2 (D)eγt = (D − α)(D − β)2 (D − γ)2 eγt = (D − α)(D − β)2 0 = 0,であるから

p2 (D)x = C3 (γ − α)(γ − β)2 eγt
これが 0 に等しいので、C3 = 0. ゆえに

x(t) = A1 eαt + B1 eβt + B2 teβt + C1 eγt + C2 teγt .
ℓ を 2 から 1 に下げて同様の議論が出来て
0 = p1 (D)x = C2 (γ − α)(γ − β)2 eγt ,∴ C2 = 0.

0 = p0 (D)x = C1 (γ − α)(γ − β)2 eγt ,∴ C1 = 0.

ℓ を 0 に下げて

92

F.8

線形非同次方程式の解空間の構造

非同次方程式の特解 z が一つでも求まれば、z と対応する同次方程式の一般解の和が非同次方程式の一般解となる、というのはよく知られている。きちんと書くと次のようになる。


定理 F.11 (非同次方程式の一般解 = 非同次方程式の任意の特解+同次方程式の一般解)
I は R の区間、n ∈ N, A : I → Rn×n と b : I → Rn は連続とする。

X := {x | x′ (t) = A(t)x(t) + b(t) (t ∈ I)} ,X0 := {z | z ′ (t) = A(t)z(t) (t ∈ I)}
とおくとき(∀y ∈ X) X = y + X0 .
(すなわち X = {y + z | z ∈ X0 }.)



証明



′

y ∈ X, すなわち y (t) = A(t)y(t) + b(t) (t ∈ I) と仮定する。

y + X0 ⊂ X であることの証明 x ∈ y + X0 とすると、ある z ∈ X0 が存在して、x = y + z.
このとき
x′ (t) = y ′ (t) + z ′ (t) = A(t)y(t) + b(t) + A(t)z(t) = A(t) (y(t) + z(t)) + b(t) = A(t)x(t) + b(t).
ゆえに x ∈ X.

X ⊂ y + X0 であることの証明 一方、x ∈ X とするとき、z := x − y とおくとz ′ (t) = x′ (t) − y ′ (t) = A(t)x(t) + b(t) − (A(t)yt + b(t)) = A(t) (x(t) − y(t)) = A(t)z(t).
ゆえに z ∈ X0 . ゆえに x = y + z ∈ y + X0 .
以上より X = y + X0 .

G

力学系についてのメモ

平衡点、安定、不安定などの言葉が良く出て来る。基本的なことは知っておくべきである。

G.1

力学系とは

この文書では、微分方程式として

dx(t) = f (x(t), t)dtという形のもの (1 階正規形微分方程式) を考えた。特別な場合として、関数 f (x, t) が t に依らない場合がある。その場合は単に f (x) と書けば良いので、微分方程式は

(G.1)

dx
dx
(t) = f (x(t)) (手短に書くと= f (x))dtdtという形になる。

93

これを自励系微分方程式 (autonomous diﬀerential equation, autonomous system, 自律系とも訳される) という。力学系 (dynamical system) と呼ばれることも多い。
一つ注意が必要である。力学系という言葉は、微分方程式でない場合にも使われる (広い意
味を持っている、ということである)。ネットで検索したりすると混乱するかもしれない。
上で「特別な場合」と書いたが、応用上出て来る微分方程式で力学系であるものは非常に多い。つまり「完全に一般ではなく、特別ではあるが、応用上は十分一般性が高い」。実際、この文書に現れる微分方程式で力学系でないものは、強制振動の方程式 (4.9) くらいである。従って、力学系に焦点を当てた議論は学ぶ価値がある。

G.2

平衡点とは

a が (G.1) の平衡点 (equilibrium point, equilibrium solution) または不動点 (ﬁxed point) であるとは、
f (a) = 0を満たすことをいう。
このときx(t) := a (t ∈ R)
で定めた x は、(G.1) の解である。
相空間で a を通る解軌道が 1 点だけからなる、ということになる。

G.3

平衡点の安定性、漸近安定性



定義 G.1 (平衡点の安定、不安定) (G.1) の平衡点 a が (リャプノフの意味で) 安定 (stable) であるとは、(∀ε > 0)(∃δ > 0)(∀x : (G.1) の解かつ |x(0) − a| < δ)(∀t ≥ 0) |x(t) − a| < εを満たすことをいう。安定でないことを「不安定である」という。


つまり、任意の正の数 ε に対して、a に十分近いところから出発した任意の解は、a から距
離 ε の範囲に止まる、ということである。

(私が学生で、このあたりのことを勉強したとき、「リャプノフの意味で安定」というのを見て、「そうでない意味で安定というのは、例えばどういうの？」と気になった。でもそういうのを目にすることはなく、最近は「リャプノフの意味で」というのは省略されるのが多くなった。時間が経って、言葉が定着したということなのだろう。)


定義 G.2 (平衡点の漸近安定性) (G.1) の平衡点 a が漸近安定 (asymptotically stable) であるとは、a が (リャプノフの意味で) 安定であり、かつ(∃δ > 0)(∀x : (G.1) の解かつ |x(0) − a| < δ)

lim x(t) = at→∞を満たすことをいう。




94

G.4

線形安定性解析(G.1) の平衡点 a の安定性の判定法として、次の定理が使われることが非常に多い。



定理 G.3 f は C 1 級とする。(G.1) の平衡点 a について、次が成り立つ。

(1) f の a におけるヤコビ行列 f ′ (a) のすべての固有値の実部が負ならば、a は漸近安定である。
(2) f の a におけるヤコビ行列 f ′ (a) の固有値のうち、一つでも正の実部を持つならば、aは不安定である。
「ヤコビ行列って何ですか？」(学生が持っている教科書の索引にヤコビ行列がない…ぶつぶつ ) Ω ⊂ Rn , a ∈ Ω, f : Ω → Rm は微分可能とするとき、f の a におけるヤコビ行列とは、
m × n 型の行列


∂fi
′
(a)
f (a) =
∂xjのことをいう。
力学系 dx= f (x) においては m = n であることに注意しよう (微分方程式の左辺は n 次元,
dt
右辺は m 次元で、それが一致するから)。ゆえに f ′ (a) は n 次正方行列で、(重複度を込めて数えて) n 個の固有値を持つ。行列 f ′ (a) の成分は実数であるが、固有値には虚数が現れることもある。
この定理が、n = 1 の場合にも使えることを注意しておく。n = 1 のとき、f の a におけるf (a + h) − f (a)def.
そのものである。またヤコビ行列は、f の a における微分係数 f ′ (a) = limh→0h
その固有値は、f ′ (a) (これは実数) である (一般に実数 A を、1 × 1 型の実行列とみなすとき、
Ax = Ax が x = 1 に対して成り立つので、A は A の固有値で、1 が固有ベクトルである。)。
ゆえに、定理を n = 1 の場合に限定すると、「f ′ (a) < 0 ならば a は漸近安定、f ′ (a) > 0 な
らば a は不安定」ということになる。
定理が成り立つのを納得したい、という人に向けて: f が a で微分可能であるとは

f (x) − f (a) − f ′ (a)(x − a)=0x→akx − aklimが成り立つことを意味する。ゆえに、x が a に近いとき

f (x) ≒ f ′ (a)(x − a) + f (a)が成り立つということは知っていると思う (m = n = 1 ならば、関数のグラフは、その接線のグラフに近い、ということ)。a が平衡点であれば f (a) = 0 であるから

f (x) ≒ f ′ (a)(x − a)ということになる。すると微分方程式は

dx(t) = f ′ (a)(x(t) − a)dtで近似できる。

e(t) := x(t) − a,x95

A := f ′ (a)

とおくと

de
x
(t) = Ae
x(t).
dt

(G.2)

これは定数係数線形常微分方程式と呼ばれる方程式である。これについての “常識的事項”
を知ると、上の定理が感覚的に納得できると思われるので、次項にまとめておく。

G.5

復習: 定数係数線形常微分方程式

A ∈ Rn×n (つまり、A は実数を成分とする n 次正方行列) とする。このとき
dx(t) = Ax(t)
dt

(G.3)

について調べよう。この形の微分方程式は定数係数線形常微分方程式と呼ばれる。(G.3) の解全体の集合を (G.3) の解空間と呼ぶ。(G.3) の解空間は n 次元線形空間の構造をもつ。その基底を (G.3) の解の基本系とか基本解と呼ぶ。
初期条件

x(0) = x0(G.4)
を満たす解は一意的に存在して、それはx(t) = etA x0である。ここで

e

tA

:=

∞
X
1
n=0

n!

t n An .

(任意の複素数 t に対して、この右辺の級数は収束する。)
0 ∈ Rn は (G.3) の平衡点である。


定理 G.4 (etA → 0 の条件) A ∈ Cn×n に対して、次の (i), (ii) は互いに同値である。

(i) lim etA = 0.
t→+∞(ii) A のすべての固有値 λ は Re λ < 0 を満たす。






系 G.5 (max Re λ < 0 ならば 0 は漸近安定な平衡点) A ∈ Rn×n のすべての固有値がRe λ < 0 を満たすならば、0 は (G.3) の漸近安定な平衡点である。






定理 G.6 A ∈ Cn×n の固有値 λ で、Re λ > 0 を満たすものが存在するならば、あるx0 ∈ Cn が存在して、 lim etA x0 = +∞. A ∈ Rn×n の場合は x0 ∈ Rn と取れる。
t→+∞







系 G.7 (max Re λ > 0 ならば 0 は不安定な平衡点) A ∈ Rn×n の固有値 λ で、Re λ > 0を満たすものが存在するならば、0 は (G.3) の不安定な平衡点である。


証明には、A の Jordan 標準形を用いるのが簡単である。行列の Jordan 標準形は、線形代数の話題であるが、時間や紙数の事情から省略されたり、はしょった説明を強行されたりする96

ことが多い。くわしい (計算例も多くて参考になる) 本として、杉浦・横沼 [42] を紹介しておく。この本は定数係数線形常微分方程式にも詳しい。
証明のメインの部分だけならば、桂田 [43] で読むことができる。
一般に、安定であるが、漸近安定ではない平衡点のことを中立安定な平衡点と呼ぶことがある。(G.3) について、0 が中立安定であるための必要十分条件も知られているが (桂田 [43] に書いておいた)、線形安定性解析には役立たないので (平衡点が、線形化した微分方程式で中立安定であっても、もとの微分方程式で不安定となることがありうる)、ここでは省略する。

H

本文中の図の描き方

コンピューターで使うツールは決定版というものはなく、またどんどん進歩していくので、何か一つのやり方を強制する気はまったくない。色々工夫するとよい。
とりあえず私はどうやったかを説明する。もっと便利 (より簡単とか、結果が分かりやすいとか) な方法を知っていたら、教えて下さい。

H.1

図 1 の描き方

Mathematica でグラフを描いてみる

a = 1; t0 = 0; T = 1;
gMalthus = Plot[
Evaluate[Table[x0 Exp[t - t0], {x0, 0, 1.6, 0.2}]], {t, t0, T},
PlotLegends -> "Expressions", AxesLabel -> Automatic]
Export["malthus.pdf", gMalthus](Evaluate[] は省略できるが、そうすると曲線の色がみな同じになる。)


ゼミで gnuplot で描いた人がいた。gnuplot> に対して



plot [0:2] 0.2*exp(x)replot [0:2] 0.4*exp(x)replot [0:2] 0.6*exp(x)replot [0:2] 0.8*exp(x)replot [0:2] 1*exp(x)replot [0:2] 1.2*exp(x)replot [0:2] 1.4*exp(x)replot [0:2] 1.6*exp(x)replot [0:2] 1.8*exp(x)replot [0:2] 2*exp(x)なるほど。繰り返しの命令 for [] を使って

plot [0:2] 0.2*exp(x)replot for [i=2:10] 0.2*i*exp(x)とできるけれど、これでは凡例が分かりにくい。
(2022/4/16 追記) 解決策を見つけた。title に sprintf() で作った文字列を指定する。


plot [0:2] for [i=1:10] 0.2*i*exp(x) title sprintf("%.1f*exp(x)", i*0.2)




97

16

1

0.2*exp(x)0.2*i*exp(x)0.2*i*exp(x)0.2*i*exp(x)0.2*i*exp(x)0.2*i*exp(x)0.2*i*exp(x)0.2*i*exp(x)0.2*i*exp(x)0.2*i*exp(x)

14

1

12

1

10

1

8

1

6

1

4

1

2

1

0

1

0

0.5

1

1

1

1

1.5

1

2

1

図 20: 図 1 の gnuplot バージョン, 凡例がイマイチ(http://takeno.iee.niit.ac.jp/~foo/gp-jman/data/20170823/gnuplot-ja-div/node187.
html を参考にしました。)

H.2

図 2 の描き方

gnuplot を使って描いた。
1. 次のようなファイルを作る。
euler.data# FIRSTN,LASTN,r=#
10 1.245394e-01
20 6.498412e-02
40 3.321799e-02
80 1.679689e-02
160 8.446252e-03
320 4.235185e-03
640 2.120621e-03
1280 1.061069e-03runge-kutta.data# FIRSTN,LASTN,r=#10 2.084324e-06
20 1.358027e-07
40 8.666190e-09
80 5.473058e-10
160 3.438538e-11
320 2.153833e-12
640 1.358913e-13
1280 1.154632e-14








euler.dat を作るには、https://m-katsurada.sakura.ne.jp/labo/text/num-ode/
node27.html にある reidai5-2.c が使える。


% cc reidai5-2.c% ./a.out > euler.data
10 1280 2%

参考までそれを Python, Julia に書き換えたものを載せておく。

98euler-error.py# euler-error.py ---dx/dt=x (0<t<1), x(0)=1 を Euler 法で解く
import math
# 微分方程式 x’=f(t,x) の右辺の関数 f の定義
def f(t,x):
return xa=0.0; b=1.0
e=math.exp(1.0)
x0=1.0str=input("N0,N1,r=")s=str.split()
N0=int(s[0])N1=int(s[1])r=int(s[2])print(’#’)
n=N0
while n <= N1:
h=(b-a)/n# 開始時刻と初期値のセット
t=ax=x0
for i in range(0,n):x += h * f(t,x)
t += hprint(’%4d %e’ % (n,abs(e-x))) # C 言語の printf("%4d %e\n", n, fabs(e-x)) 相当
n *= r




99



euler-error.jl# euler-error.jl --- dx/dt=x (0<t<1), x(0)=1 を Euler 法で解く
# euler-error.c の真似
using Printffunction f(t,x)
xendfunction compute_euler_error(N0, N1, r)
a=0.0; b=1.0;
e=exp(1.0);
x0=1.0;n=N0;while (n <= N1)
h=(b-a)/n;
t=a;
x=x0;
for i=0:nx += h * f(t, x);t += h;end@printf("%4d %e\n", n, abs(e-x))n *= r;end
end# julia euler-error.jl と実行した場合に実行する。
if abspath(PROGRAM_FILE) == @__FILE__
s=split(Base.prompt("# FIRSTN,LASTN,r="))N0=parse(Int64,s[1]);N1=parse(Int64,s[2]);r=parse(Int64,s[3]);println("#")
compute_euler_error(N0,N1,r)
end








euler-error.py を実行する% python euler-error.py > euler.data
10 1280 2% cat euler.data# FIRSTN,LASTN,r=: #10 1.348349e-01
20 6.768076e-02
40 3.390861e-02
80 1.697167e-02
160 8.490220e-03
320 4.246211e-03
640 2.123381e-03
1280 1.061760e-03
%



100





euler-error.jl を実行する% julia euler-error.jl > euler.datarunge-kutta.dat を作るプログラムは公開していないが、reidai5-2.c を書き換えれば良い。Runge-Kutta 法のコードは、https://m-katsurada.sakura.ne.jp/labo/text/
num-ode/node28.html が参考になる (分かってしまえば簡単, ほぼ該当部分のコピペで済む)。
2. gnuplot を起動して次のようにすれば図が (画面に) 描ける。


gnuplot で描く



set logscale xy
plot "euler.data" with lp
replot "runge-kutta.data" with lp(lp の l は L の小文字。lp は linespoints の略。)


こうして表示されたグラフを保存するには、続いて


gnuplot で描く (続き)

set term pdf
set output "compare_data.pdf"
replotとすればよい。ここでは PDF 形式で保存したが、set term png とか set term postscript
eps color とか、形式は色々なものが選択可能である。
補足 1 ここでは標準出力のリダイレクション (コマンド > ファイル名) を使って、データ・ファイル euler.data を作成したが、C のプログラムの中で fopen(), fprintf(), flose()などの関数を使って直接ファイルを作成することも出来る (その方がずっと融通が効く)。
「常35
微分方程式の初期値問題を解くプログラムの書き方 2.3.1 Euler 法の C プログラム例」 の問2 を参考にすると良い (解答プログラム付き)。
補足 2 C のプログラムから直接 gnuplot を起動して、データを渡してグラフを描くことも出来る。

• 「常微分方程式の初期値問題を解くプログラムの書き方 2.3.1 Euler 法の C プログラム
36の問 3 を参考にすると良い (解答プログラム付き)。
例」
この例では、データは外部ファイルを介して渡している。
• 「gnuplot 入門 C から gnuplot を呼び出す」37
外部ファイルを使わず、すべてパイプを使って処理する例がある。

H.3

matplotlib で日本語タイトルを用いる

日本語フォントを登録して、それを利用するだけである…のはずなのだが、具体的にどうすれば良いのかはケース・バイ・ケースである。OS により用意されているフォントが違うし、どうやってインストールした Python を用いるかでファイルの置き場所が異なる。
35

https://m-katsurada.sakura.ne.jp/labo/text/intro-ode-simulation/node14.html
https://m-katsurada.sakura.ne.jp/labo/text/intro-ode-simulation/node14.html
37
https://m-katsurada.sakura.ne.jp/labo/howto/intro-gnuplot/node21.html
36101以下は、私の場合である。私は Mac で、MacPorts でインストールした python310 を使っている。macOS にはヒラギノ・フォントが用意されている。
次の 2 つのうち、どちらかをすれば良い。一度やれば済むという点で、1 がお薦め。

1. matplotlibrc で

font.family: Hiragino Maru Gothic Proという設定を書き込む。MacPorts で Python をインストールした場合は、matplotlibrc
は ~/.matplotlib/matplotlibrc に置かれる (らしい)。どこに置いてあるかは、一般には次のようにして調べられる。


Python で matplotlibrc のパス名を表示させる

import matplotlib
matplotlib.matplotlib_fname()

2. Python のプログラム中で、例えば次のようにフォントを指定する。


plt.rcParams[’font.family’] = ’sans-serif’
plt.rcParams[’font.sans-serif’] = [’Hiragino Maru Gothic Pro’, ’Yu Gothic’,’Meirio’, ’Takao’, ’IPAexGothic’, ’IPAPGothic’, ’VL PGothic’, ’Noto Sans CJK JP’]

これは某 WWW サイトに載っていたものだが (どこのサイトだったか、メモをなくしてしまった…)、色々なフォントを探すようにしてある。Mac の場合は最初の ’Hiragino Maru
Gothic Pro’ で解決する。Meirio や Yu Gothic は Windows で有効なもの？ IPAexGothic
や IPAGothic は IPA (独立行政法人情報処理推進機構) が開発・公開しているフリーのフォントである。IPA フォントを自分で入手・設定すれば、Windows, Mac, Linux のいずれでも利用可能である (Linux を使っている人がインストールしていそうだ)。
以上、これで良いと思う。
この (とりあえずの) 結論に到達するまでは、IPA フォントを利用するようにしていた。そのときに書いたメモを参考まで以下に引用しておく。
私は MacTeX を使う関係で、IPAexGothic (ipaexg.ttf) が Mac にあるので、それを利用することにした。


ダメ元で Mac にないか、ターミナルで検索

mdfind ipaexg.ttf




私の場合、
/usr/local/texlive/2021/texmf-dist/fonts/truetype/public/ipaex/ipaexg.ttfが見つかった。
これを Python で使っている matplotlib から見えるようにする。
私の場合は、普段 MacPorts の Python を使っているため、/opt/local/Library/Frameworks/Python.framework/Versions/3.10/lib/python3.10/site-packages/matplotlib/mpl-data/fonts/ttf/というフォルダがあり、そこに ipaexg.ttf をおけば動いた。ただしキャッシュもクリアしておく必要がある。

102





sudo ln -s /usr/local/texlive/texmf-local/fonts/truetype/cjk-gs-integrate/ipaexg.ttf \/opt/local/Library/Frameworks/Python.framework/Versions/3.10/lib/python3.10/site-packages/matplotlib/mpl-data/fonts/ttf
rm ~/.matplotlib/*.json

後は、このフォントを使うように指定するだけである。~/.matplotlib/matplotlibrc を
用意して、




font.family: IPAexGothicと書き込んだ。
きちんとやるには、デフォールトの設定ファイルを探し出して、それを~/.matplotlib/matplotlibrcにコピーしてから、上の 1 行を書き足す。デフォールトの設定ファイルは注釈行しかない場合も多く (ユーザーのための参考情報であって、実質的に何も書いていないのと同じ。)、その場
合は上の 1 行だけのファイルを作れば良い。

I

微分多項式についての注意

I.1

個人的な見解・意見

定数係数線形常微分方程式

a0 x(n) (t) + a1 x(n−1) (t) + · · · + an−1 x′ (t) + an x(t) = f (t) (t ∈ I)について論じるときに (I は R のある区間)、微分多項式 (diﬀerential polynomial) が導入されることが多い。
微分作用素

(I.1)

D :=

d
dt

を導入して、f (λ) = a0 λn + a1 λn−1 + · · · + an−1 λ + an ∈ C[λ],(I.2)

x ∈ C n (I)

に対して

(I.3)

f (D)x := a0 x(n) + a1 x(n−1) + · · · + an−1 x′ + an xと定める。f (D)x ∈ C(I) であり、f (D) : C (n) (I) → C(I) とみなせる。
簡単のため、2 階微分方程式の場合、(D − α)(D − β)x = 0という式を考えてみよう。一見単純ではあるが、何かルールを定めないと解釈にあいまいさのある式である。実際、次のような複数の解釈が考えられる。

(i)(D − α) ((D − β)x) = 0と解釈するのか。

103

(ii) f1 (λ) := λ−α, f2 (λ) := λ−β として、f1 (D) = D −α と f2 (D) = D −β の積 f1 (D)f2 (D)を適当に定めることにして、(f1 (D)f2 (D)) x = 0と解釈するのか。f1 (D)f2 (D) としては(a) (f1 (D)f2 (D)) x = f1 (D) (f2 (D)x) (x ∈ C 2 (I)) で定まる f1 (D)f2 (D) : C 2 (I) → C(I).
(b) f (λ) := λ2 − (α + β)λ + αβ (つまり (λ − α)(λ − β) を展開・整理してできる多項式)としたときの f (D) : C 2 (I) → C(I).
微分多項式を扱っている常微分方程式のテキスト・解説はたくさんあるけれど、• どう解釈するのか明言されていないものがある。
• どう解釈するか記されていても、その後で論理が飛躍していると思われるものがある。
この記号がどのように使われるか、典型例を一つあげてみよう。



f (λ) = λ2 + pλ + q ∈ C[λ] が f (λ) = (λ − α)(λ − β) と因数分解されるとき、f (D)x = 0を満たす x に対して x1 := (D − β)x とおくと(D − α)x1 = (D − α)(D − β)x = f (D)x = 0が成り立つので、ある定数 C が存在して x1 (t) = Ceαt が成り立つ。


これを見ると (D − α) ((D − β)x) としているように思われるので、解釈 (i) が有望そうであるが、それを f (D)x に等しいとしているので、(D − α)(D − β) = f (D) とみなしているようでもある。
どのように定義するかは、正しい答えが一つだけある、というのではないことに注意しよう。数学の概念の多くは、複数の仕方で定義することが出来て、ある流儀で定義したとき、別の流儀の定義に現れる式が定理となったりする。全体として同じことができれば構わない。

I.2

高橋 [1] の記述を検討する



以下，

(2.7)

高橋 [1] p. 39

dを D と略記する。このとき，dt(D − α)(D − β)u = (D2 − (α + β)D + αβ)u = (D − β)(D − α)uが成立する．実際， 左辺 = (D − α)Du − βu) = D(Du − βu) − α(Du − βu) = D2 u −βDu − αDu + αβu = 中辺. 対称性より，中辺 = 右辺．したがって，（7）から帰納的に，次の等式が成立することがわかる：(2.8)

P (D)a ≡ Dn u + c1 Dn−1 u + · · · + cn−1 Du + cn u = (D − λ1 )n1 · · · (D − λr )nr u.

さらに，この右辺の，(D − λp )np たちの順序を入替えても P (D) に等しい．


これは (i) の解釈をした上で、(D − α) ((D − β)x) = P (D)x104

であることを証明していると考えられる。
このテキストでは、(この後も) 微分多項式は必ず関数に作用させる (右に関数が書いてある) 形で登場し、(D − α)(D − β) = f (D)のような関数抜きの式は見当たらない。
(D − λr )nr u という式が出て来るが、これはどう解釈すべきだろうか。
(D − λr )j+1 u = (D − λr ) (D − λ)j u(j = 1, · · · , nr − 1)
により定義するのかな？
当たり前であるが、きちんと書かれている。

参考文献[1] 高橋陽一郎：微分方程式入門, 東京大学出版会 (1988), 丸善 eBook にある。https://elib.maruzen.co.jp/elib/html/BookDetail/Id/3000046850.
[2] 桂田祐史：常微分方程式の初期値問題の数値解法, https://m-katsurada.sakura.ne.
jp/labo/text/numerical-ode.pdf (1994〜).
[3] Wikipedia, : https://en.wikipedia.org/wiki/Logistic_function.
[4] 桂田祐史：常微分方程式の初期値問題を解くプログラムの書き方, https://m-katsurada.
sakura.ne.jp/labo/text/intro-ode-simulation/ (2021/4/23).
[5] M. ブラウン：微分方程式 下, シュプリンガー・フェアラーク東京 (2001).
[6] 桂田 祐史：微 分 方程 式入門, https://m-katsurada.sakura.ne.jp/lecture/kiso4/
kiso4ode.pdf, 基礎数学 IV を講義したときの講義ノート (2004).
[7] 今隆助, 竹内康博：常微分方程式とロトカ・ヴォルテラ方程式, 共立出版 (2018), https://elib.maruzen.co.jp/elib/html/Viewer/Id/3000086446.
[8] Murray, J. D.: Mathematical Biology I: An Introduction, Interdisciplinary Applied Mathematics, Springer-Verlag (2002), マレー [44] が邦訳.
[9] Gause, G. F.: Experimental studies on the struggle for existence: 1. Mixed population of
two species of yeast, Journal of Experimental Biology, Vol. 9, pp. 389–402 (1932).
[10] Hutchinson, G. E.: The paradox of the plankton, The American Naturalist, Vol. 95, pp.
137–145 (1961).
[11] May, R. M. and Leonard, W. J.: Nonlinear Aspects of Competition Between Three
Species, SIAM Journal on Applied Mathematics, Vol. 29, pp. 243–253 (1975).
[12] Jaramillo, G., Mrad, L. and Stepien, T. L.: Dynamics of a linearly perturbed MayLeonard competition model, Chaos, Vol. 33, No. 6 (2023), arxiv に置いてある: https://arxiv.org/pdf/2210.04342.pdf.
[13] 河井智康：イワシと逢えなくなる日 — 5 億年の結晶『魚種交代』の謎に迫る, 情報セン
ター出版局 (1988/11/6), 文庫版 (角川ソフィア文庫, 2001/8/25) には「文庫版あとがき」がある。 入手してある。
105

[14] 河井智康：消えたイワシからの暗号, 三五館 (1999/7/30).
[15] 川崎健：イワシと気候変動 — 漁業の未来を考える, 岩波新書, 岩波書店 (2009/6/19), 1983
年レジーム・シフトを提唱した川崎先生の啓蒙本.
[16] Matsuda, H., Wada, T., Takeuchi, Y. and Matsumiya, Y.: Model analysis of the eﬀectof environmental ﬂuctuation on the species replacement pattern of pelagic ﬁshes underinterspeciﬁc competition, Researches on Population Ecology, Vol. 34, No. 1, pp. 309–319(1992).
[17] Takeuchi, Y., Oshime, Y. and Matsuda, H.: Persistence and Periodic Orbits of a ThreeCompetitor Model with Refuges, Mathematical Biosciences, Vol. 108, pp. 105–125 (1992).
[18] 松田裕之：環境生態学序説, 共立出版 (2000/12/1).
[19] Kermack, W. O. and McKendrick, A. G.: A Contribution to the Mathematical Theory
of Epidemics, Proceedings of the Royal Society of London. Series A, Vol. 115, No. 772,
pp. 700–721 (1927).
ふさお[20] 佐藤總夫：自然の数理と社会の数理 II, 日本評論社 (1987), 版元在庫切れだったりします
が、図書館で借りるとか、 アマゾンで中古をぽちするとか、 読む方法はあると思います。
ゼミ生は見つからなければ相談して下さい。
ひさし[21] 稲葉 寿 ：感染症の数理モデル, 培風館 (初版 2008, 増補版 2020/12/15).
[22] 神永正博：Python と実例で学ぶ微分方程式 〜はりの方程式から感染症の数理モデル
まで〜, コロナ社 (2021/10/6), https://elib.maruzen.co.jp/elib/html/Viewer/Id/
3000112217.
[23] Harko, T., Lobo, F. S. N. and Mak, M. K.: Exact analytical solutions of the SusceptibleInfected-Recovered (SIR) epidemic model and of the SIR model with equal death andbirth rates, Applied Mathematics and Computation, Vol. 236, pp. 184–194 (2014).
こうじ[24] 笠原晧司：微分方程式の基礎, 数理科学ライブラリー, 朝倉書店 (1982).
[25] 笠原晧司：新微分方程式対話, 日本評論社 (1981/1/1, 改訂版 1992/3/24).
[26] 竹之内脩：常微分方程式, ちくま学芸文庫, 筑摩書房 (2020/12/10), 1981 年に秀潤社から出版された書籍の文庫化。Kindle もある。
[27] M. ブラウン：微分方程式 上, シュプリンガー・フェアラーク東京 (2001).
[28] 俣野博：常微分方程式入門 — 基礎から応用へ —, 岩波書店 (2003), 俣野 博, 微分方程式
II, 岩波講座 応用数学, 岩波書店 (1993) の単行本化． 2015/6/10 にオンデマンドブック
スとして復刊された．
[29] 高橋陽一郎：力学と微分方程式, 岩波書店 (2020/1/10), 岩波講座 現代数学への入門「力学と微分方程式」(1996/05/29) の単行本化。 丸善 eBook にある。https://elib.maruzen.
co.jp/elib/html/BookDetail/Id/3000007597.
[30] コディントン・レヴィンソン：常微分方程式論 上, 下, 吉岡書店 (1968, 1969).
106

たけとも

[31] 三井斌友：数値解析入門, 朝倉書店 (1985).
[32] 三井斌友：常微分方程式の数値解法, 岩波書店 (2003), 「微分方程式の数値解法 I」岩波講座応用数学 (1993) の単行本化.
たけとも

こ と う としゆき

よしひろ

[33] 三井斌友, 小藤俊幸, 齊藤善弘：微分方程式による計算科学入門, 共立出版 (2004/2/25).
[34] E. ハイラー, S. P. ネルセット, G. ヴァンナー：常微分方程式の数値解法 I 基礎編, シュプ
リンガー・ジャパン (2007).
[35] E. ハイラー, G. ヴァンナー：常微分方程式の数値解法 II 発展編, シュプリンガー・ジャ
パン (2008).
[36] E. Hairer, C. L. and Wanner, G.: Geometric Numerical Integration: Structure-Preserving
Algorithms for Ordinary Diﬀerential equations, Springer Series in Computational Mathematics, Springer (2010/3/11), ISBN: 978-3642051579, http://www.unige.ch/~hairer/
software.html.
[37] 一松信：微分積分学入門第二課, 近代科学社 (1990).
[38] 小川知之, 宮路智行：数理モデルとシミュレーション, サイエンス社 (2020/12/19).
[39] 坂井秀隆：常微分方程式, 大学数学の入門, 東京大学出版会 (2015), 内容が非常に豊富.
やや粗いけれど優れたテキスト. 丸善 eBook では https://elib.maruzen.co.jp/elib/
html/BookDetail/Id/3000028058 でアクセス可能である.
[40] 高野恭一：常微分方程式, 朝倉書店 (1994), 線形微分方程式のモノドロミー表現やフックス型微分方程式等、複素領域における微分方程式の話が載っているのが特徴.
[41] 桂田祐史：常微分方程式の初期値問題の解の延長, https://m-katsurada.sakura.ne.
jp/labo/text/continuation_of_solution.pdf (2022/3/19).
[42] 杉浦光夫, 横沼健雄：Jordan 標準形・テンソル代数, 岩波書店 (1990), この前半は、杉浦光夫, Jordan 標準形と単因子論 I, II, 岩波講座 基礎数学 (1976,1977) が元になっている。
[43] 桂田祐史：定数係数線型常微分方程式の解の漸近挙動, https://m-katsurada.sakura.
ne.jp/labo/text/exptA.pdf (2022/3/2).
[44] ジェームズ D. マレー：マレー数理生物学入門, 丸善出版 (2014/1/18), 三村 昌泰 (監修,翻訳), 瀬野 裕美 (監修), 河内 一樹 (監修), 中口 悦史 (監修).

107