第3回

内容

• 1階常微分方程式の数値解法（オイラー法）

• 2階常微分方程式の数値解法（オイラー法）

目標

• 1階常微分方程式を解く方法を知る

• N階常微分方程式についてもN連立1階常微分方程式に帰着できること知る

放射性炭素・年代測定法

前回の常微分方程式について少し解説します。次の常微分方程式でした。

また，この常微分方程式の一般解が次のようになることから，

N(t) = N₀ e^-λt (S)

このグラフを 0 ≦ t ≦ 10 の範囲で描画し，さらに，

N(τ) = N₀/2

となる　τ　は log2/λ ですが（τ に N₀ が含まれない事に注目），それを近似的に求めることをしました。さて，この方程式 (P) は様々な現象のモデル方程式と見ることができます。先の τ を求める問題との関係から言えば，(P) は放射性炭素・年代測定法との関係を説明するのが適切でしょう。

放射性炭素・年代測定法は，W.F. Libby が 1940年代後半に開発した方法で（1960年にノーベル化学賞を受賞），今世紀初めに英国の実験物理学者ラザフォード等によって発見・開発された放射能の原理に基づくものです。ある種の原子はもともと不安定であるため，外部からの影響が無くとも時間が経つと放射性物質を放射し，他の新しい元素に遷移します。実験結果から，ラザフォードは，放射性物質のいくつかについて，単位時間あたりの元素の崩壊個数は現在の原子数に比例する事を示しました。すなわち，時刻 t におけるある放射性物質中の原子の個数を N(t) とすると，

dN/dt = -λN

が成り立ちます。ここで，λは正の定数で崩壊定数と呼ばれ，物質毎に決まるものです。ある特定の物質のλは実験によって定める必要があります。実際には，実験的に崩壊定数を求めるのではなく，半減期τを測定し，それから崩壊定数λを決めるようです。

時刻 0 における，原子の個数が N₀ であるとき，つまり (P) を解けば，

N(t) = N₀ e^-λt

が解として求まります。したがって， t = τ において，N = N₀/2 ならば，

N₀/2 = N₀ e^-λτ

よって，

τ = log2/λ

が分かります（τ に N₀ が含まれない事に注意。つまり初期の放射性物質の量によらないということ）。実験的に知られているいくつかのサンプルの半減期を以下に示します。

放射性炭素・年代測定法においては，炭素14が重要です。これについて，崩壊定数を求めると，

λ = log2/τ = 1.245×10^-4 /年

となります。さて，放射性炭素・年代測定法がなぜに成り立つかを少しだけ説明しましょう。

自然界には3種類の炭素があります。それぞれ，炭素12，炭素13，炭素14と呼ばれ，それぞれの存在比率は

炭素12：98.9%

炭素13：1.1%

炭素14：0.00000000012%　（おおよそ一兆分の一）

です。また，炭素12，炭素13は安定な物質だが，炭素14は不安定な物質で，時間と共に放射線を放出して窒素14になり，その半減期は，前述の5568年です。

炭素14は，宇宙から降り注ぐ宇宙線によって，窒素14から生じます。生成された炭素14は，生成後直ちに酸素と結合し二酸化炭素として通常の二酸化炭素と共に大気中に拡散します。炭素14の生成速度と，先の崩壊速度がおおむね釣り合っており，その結果大気中の炭素14はほぼ一定と見なすことができることが知られています。（ほぼ一定ですが、様々な要因で変化しており、その変化を考慮して補正することも行われるようです。）

さて，その炭素14を含む二酸化炭素は，光合成によって通常の二酸化炭素と共に植物に取り込まれ，食物連鎖によって動物にも広がります。例えば，木の場合，木が生きている間は二酸化炭素の形で炭素14と他の炭素12，炭素13を吸収し続けるため，木に含まれる炭素14の比率は大気中と同じであると考えられます。つまり，おおよそ一兆分の一の比率で炭素14が含まれているはずです。

その後木が死ぬとどうなるか？例えば，木の加工品を作ろうと人が木を伐採するとそのとき木は死に，光合成の停止と共にあらたな炭素14の取り込みは停止します。その後は，炭素14は，その崩壊速度に従って崩壊するのみとなります。（他の炭素12，炭素13は安定な物質なので崩壊しない。）

この事実から，例えば測定したい木でできた加工品について，精密に炭素14と他の炭素12，炭素13との比率を測定したとして，その比率が例えば二兆分の一であったなら，その加工品はおおよび5500年前のものであることが分かるという理屈です。

1階常微分方程式の数値解法（オイラー法）

次のような常微分方程式の初期値問題を、オイラー法を用いて数値的に解き、matplotlib を用いて結果をグラフにしてもらいます．

次の初期値問題を考えます．

求めたいのは、ｔ＞０の範囲の関数 ｙ（ｔ）であり、関数 f(t, y) および定数 a は与えられているとします．また、 f(t ,y) の値は t, y が与えられればいつでも計算できるものとします．（プログラム中では関数として定義するとよい．）

コンピューターでは数値を離散的に（とびとびに）しか扱うことができません．そこで、（１）式をコンピューターで扱えるようにする必要があります．（１）をコンピュータで解く方法には色々とありますが、この演習では代表的なオイラー法（とルンゲ・クッタ法）を用います．その前に、準備として差分商について説明します．

差分商

f(x) を x の関数としたとき（ここでの f(x) は上記(1)式中の f とは関係ありません）、f(x) の x = a における微分係数は

で定義されます．ここで、定義中の極限操作を取り払い、h を有限にとどめた

を考えると、h を十分小さな正の実数にとれば、（３）式は（２）式の近似となっていると考えられます．この D のように、関数のいくつかの点における値の差を用いてその関数の微分係数を近似することを差分近似といい、（３）式の右辺のような量を差分商といいます．いまの場合は１階微分係数を近似する１階差分商です．

では、このような差分商を用いて（１）式を離散化してみましょう．まず、（１）式の第一式の微分係数を上記の差分商で置き換えてみます．

ただし、h は正の数とします．（４）式と（１）式の第一式は異なる方程式なので、（１）式の第一式の解 y(t) は一般に（４）式を満たしません．そこで、混乱を防ぐために（４）式の y(t) を Y(t) と書き換えましょう．

（５）式のように差分を含む方程式を差分方程式といいます．

次に、（１）式の y を Y に置き換えた初期条件

の下で（５）式を解くことを考えます．（５）式は、

と変形できるので、t = 0 を代入すると、

となり、Y(0) から直ちに Y(h) が求まります（注意：Y(0)は初期条件として与えられているので既知）．同様にして（７）式を繰り返し用いると、

というように、j=1,2,3,.... とすると Y((j+1)h) を Y(jh) から計算できることがわかります．

h をいったん決めると、 t = jh 以外の時刻の Y の値は（７）式からは求めることができません．このように、差分方程式を解くと従属変数はとびとびの時刻で値が定まります．そのようなとびとびの時刻を格子点と呼びます．Y は格子点でのみ意味があるので、そのことを明示するために Y(jh) を Yj と書き換え、 Tj = jh とすると、（５）式と初期条件は、

となり、結局（１）式の常微分方程式の初期値問題が、Yjに関する漸化式の問題（１０）式に置き換えられました．この方法をオイラー法と呼びます．

オイラー法のアルゴリズム１

解を求める時刻の上限を T とし、適当な自然数 N を定めて h=T/Nとする．

T, h, a 及び 関数 f は実数型

N, j は整数型

Y[N+1], t[N+1] は実数型配列

T, Nを設定

配列 T[N+1], t[N+1] を用意する

h = T/N

a の入力を受ける (input 関数）

Y[0] = a （初期値の設定）

j = 0,1,2,.......,N-1 の順に (for 文）

t = j*h

Y[j + 1] = Y[j] + h*f(t[j], Y[j])

以上を繰り返す

点(t[0], Y[0]), (t[1], Y[1]), ... , (t[N], Y[N]) を直線で結ぶことにより数値解をグラフとして描画 (matplotlib)

オイラー法のアルゴリズム１は，計算結果を全て保存しています。しかし、グラフ描画のために全ての計算結果を保管する必要は無く、次のオイラー法のアルゴリズム２で十分です．

オイラー法のアルゴリズム２

解を求める時刻の上限を T とし、適当な自然数 N を定めて h=T/N とする。また、一定間隔 INTV 毎に計算結果を記憶することとする。

T, h, a, Y, t 及び 関数 f は実数型

N, j は整数型

T, N, INTVを設定

空のリスト Yg, tg を用意する

h = T/N

a の入力を受ける (input 関数）

Y = a

t = 0 での値を記憶（リストに追加）

j = 0,1,2,.......,N-1 の順に (for 文）

t = j*h

Y = Y + h*f(t, Y)

一定間隔（INTV）毎に値を記憶（リスト tg, Yg に追加）

(求まった Y は時刻 t+h での値であることに注意）

以上を繰り返す

点(tg[0], Yg[0]), (tg[1], Yg[1]), ... , (tg[N/INTV], Yg[N/INTV]) を直線で結ぶことにより数値解をグラフとして描画 (matplotlib)

f(t, Y) は関数として定義するとよい．

例えば、INTV=100とすると、作成されるデータファイルの容量は1/100となる．これは，先に見たようにある程度の折れ線数で得あれば十分なめらかに見えるので，グラフを描く為のデータとして巨大になりすぎることを防ぐ為である．ここで100という数値は適当に与えた物であって，グラフが十分なめらかに見えかつデータが小さくなるような適当な値を探す必要がある．

簡単な常微分方程式とその数値解のグラフ化

さて，実際にオイラー法をつかって常微分方程式を解いてみましょう．

前回からの常微分方程式を用いて演習を進めます。(ただし N を y と書き換えた．上記アルゴリズム等の分割数 N とまぎらわしいので．)　

前回は，この微分方程式を解析的に解き，得られた解をグラフ化してもらいましたが，今回は(P)をオイラー法を用いて数値計算で解いてもらい，得られた数値解をグラフ化してもらいます（違いを確認）．

課題1（演習）

(P)をオイラー法で解き，得られた数値解をグラフとして表示するプログラムを作成せよ．ただし，

Y₀=1, λ=1, 0≦t≦10

とする。

また、第2回にて扱った(P)の解

y(t) = y₀ e^-λt (S) (Nをyで書き換えた)

を重ねて描く事により、それらの違いを確認せよ。

例えば次のようになる。（この例では、0≦t≦10 を 100等分している。つまり時間刻み幅を h と書くことにすると、 h = 10.0/100.0 = 0.1 である。）

次の動画で、input関数で初期値を得るように変更し、半減期も計算結果から求め、初期値によらないことを確かめる。最後に、オイラー法のアルゴリズム2に書き換える。最終的なCodeには半減期の部分は消去した。細々やってると長くなってしまった。

課題2（演習）

ある種の生物の増殖に関する大変有効なモデルとして次のロジスティック方程式がある．

dx/dt = rx(1 - x/K) (L)

ここに、x(t) は時刻 t における、生物集団のサイズであり，K > 0 を環境収容力と呼び，r > 0 を内的自然増加率と呼ぶ．

今、 r = 1, K = 100 とし、初期条件として、次の３種類のものを考える．

1. x(0) = 0

2. x(0) = 10

3. x(0) = 150

このとき、それぞれの初期値に対する (L) の解を 0≦ t ≦10 についてオイラー法で求め，結果をmatplotlibでグラフ化せよ．ただし、N = 10000 （時間刻み数）とする．データは100毎に出せば十分である．グラフは、横軸を t 、縦軸を x(t) とする．

なお、ロジスティック方程式については、別途調べる事。扱う方程式に興味を持つことは大事である。

動画には途中間違いもあったが、そのまま載せる。間違いも参考にして欲しい。

２階常微分方程式の数値解法（オイラー法）

つづいて、２階常微分方程式をオイラー法を用いて解いてもらいます．次の問題を考えます．

これは単振動を記述する２階の常微分方程式です．ただし簡単の為、時間 t での微分を ' を用いてあらわしています．２階の常微分方程式は次のように連立の１階の常微分方程式に帰着できます．y1(t)=y(t), y2(t)=y'(t) とおくと、（１１）式は

となります．こうなれば、前回同様オイラー法を用いて解を求めることができます．アルゴリズムは次のようになります．但し、 T を解を求める時間の上限とし、h を時間刻み幅、N を時間刻み数とします．また、変数は y1, y2 として、a1, a2 を初期値とします．また、f1(t, y1, y2) および f2(t, y1, y2) はそれぞれ関数で（１２）式の場合には f1 が第一式の右辺、 f2 が第二式の右辺となります．

このように、一般にN階の常微分方程式は、N連立の1階常微分方程式系に変換されますので、次のオイラー法のアルゴリズムを発展させることで、解く事ができます。よって、ここでは2階までしか解説しません。

オイラー法のアルゴリズム

解を求める時刻の上限を T とし、適当な自然数 N を定めて h=T/(N-1) とする．

T, h, a1, a2, y1, y2, y1_new, y2_new, t 及び 関数 f1, f2 は実数型

N, j は整数型

T,N,h を定める

空のリスト y1g, y2g, tg を用意する

y1 = a1, y2 = a2

t = 0 での値の記録 (y1, y2, t をリストに追加)

j = 0,1,2,....,N-1の順に

t = j*h

y1_new = y1 + h*f1(t, y1, y2)

y2_new = y2 + h*f2(t, y1, y2)

y1 = y1_new

y2 = y2_new

一定間隔毎にリストに追加(y1, y2, t をリストに追加)

（ここで求まったy1,y2 は時刻t+hのものであることに注意）

を繰り返す

y1g, t2g, t を matplotlib で描画

課題3（演習）

(12) をオイラー法で解き、グラフを描け．グラフは横軸が t であり、縦軸を y(t) とせよ．また、時間刻み幅 h を変えたときに結果がどのようにかわるか観察せよ．0≦t≦20の範囲において、h=0.0001,0.001, 0.01, 0.1それぞれの場合につい微分解 y(t)=cos t とどの程度あっているか各自調べよ．

色々な h に対する計算結果は次のようになる．この結果から、h は十分小さい値でないとまずいことが分かるだろう．

上の課題の結果から分かるように、オイラー法では h を十分小さくとらないと微分方程式の解をうまく近似できません．h を小さく取るということは、計算時間が多くかかることを意味します．そこでより高精度で大きな h でもうまく解を近似できるような解法が求められます．そのような要求に応えるルンゲ・クッタ法があります．（ただし、ここで扱った例が特にオイラー法に不利な例であることは注意しておく。ズレの方向が一方向で、それが蓄積してどんどん拡大する。一般的には、オイラー法で十分な場合も多い。）

課題4

下の図のように、壁に一端を固定されたバネの他端におもりがつけられている状況を考えます．

このおもりを少しだけ右に引っ張ってから手をはなす．おもりの中心の時刻 t におけるつり合いの位置からのずれを y(t) とします．また、時刻 0 でのおもりの位置を y(0)=1、おもりの速度を y'(0)=0 とすると、課題3と全く同じ、つぎの単振動の問題となります．

これではおもりに働く力はバネの復元力のみで、ある意味非現実的です．今回は復元力以外におもりの速度 y' に比例する抵抗が働く場合を考えましょう（おもりと床との間の摩擦力のようなものをイメージ）．改良されたモデルは次のようになります(2ky' の項が新たに加わったもので、摩擦のようなものをあらわしている)．

ここで k,ω はそれぞれ抵抗力、復元力の強さを表す正の定数です．課題は、（１４）式の問題を連立１階常微分方程式に変換し、オイラー法を用いて解き、matplotlib をつかってグラフにしてください．k < ω でのグラフと k > ω でのグラフを見比べて違いに関して考察してみて下さい．h を色々と変化させると数値解はどう変化するか？

課題5

次の、ロトカ・ボルテラ競争系について、図書館の書籍、WEB上の情報などをもちいて調べまとめよ．どういった現象を表すモデルであるか？

dx/dt = r1 x(1 - (x + ay)/K1)

dy/dt = r2 y(1 - (bx + y)/K2)

課題6

課題5中のロトカ・ボルテラ競争系モデルを数値計算を用いて解け．

課題7

（こちら、特に数値計算を目的とする人は、問題3をやって下さい。問題1は既に課題2でやりました。

ある種の生物の増殖に関する大変有効なモデルとして次のロジスティック方程式がある．

dx/dt = rx(1 - x/K) (L)

ここに、x(t) は時刻 t における、生物集団のサイズであり，K > 0 を環境収容力と呼び，r > 0 を内的自然増加率と呼ぶ．

今、 r = 1, K = 100 とし、初期条件として、次の３種類のものを考える．

1. x(0) = 0

2. x(0) = 10

3. x(0) = 150

問題１：（課題2と同じ）

このとき、それぞれの初期値に対する (L) の解を 0≦ t ≦10 についてオイラー法で求め，結果をmatplotlibでグラフ化せよ．ただし、N = 10000 （時間刻み数）とする．データは100毎に出せば十分である．グラフは、横軸を t 、縦軸を x(t) とする．次に示す画像を参考にせよ．（こちらは、 gnuplot というソフトで、数値データをグラフ化したもの。興味がある人は調べると良い。）

問題２：

ロジスティック方程式(L)の解を解析的に求めよ．また得られた解が，確かに(L)の解であることも確かめよ．

ちなみに，解は以下のようになる．

x(t) = K/(1 + C *exp(-r*t))

ここに C は，初期値 x(0) から決まる定数で，

C = K/x(0) - 1

で与えられる．

問題３：

オイラー法で求めた近似解と問題2で求めた真の解を比べたい．h を小さくすると近似解は真の解に近づくであろうがどのように近づくであろうか．問題２で求めた真の解を x(t)，刻み幅 h によってオイラー法で求めた近似解を X(t;h) と書くことにする．今，初期値は 10 すなわち x(0) = 10, X(0;h) = 10 とする．次のような尺度をもって近似解の正確さはかることにする．

E(t;h) = |x(t) - X(t;h)|

E(t;h) は刻み幅が h である時の時刻 t における真の解と近似解の差を与える．ここでは，t = 2 における E を様々な h について調べることにする．すなわち，h = 0.1, 0.01, 0.001, 0.0001 それぞれについて，E(2; 0.1), E(2; 0.01), E(2; 0.001), E(2; 0.0001) を求めグラフで検討せよ．近似解と真の解の差は h に関してどのように減少するか考えよ．

目次に戻る