第５・６回（第３週）

グローバル変数とローカル変数

これまで関数と変数について習いました。ここでは、プログラム中での変数の定義の位置によって、変数の意味合いが変わることを見ていきましょう。まずは基本となるプログラム（sample_global1.py）を考えましょう。

sample_global1.py

a = 10

print(f'a={a} in initial')

もはやこのプログラムの説明は不要でしょう。実行結果はa=10 in initialです。では次のような関数が含まれる場合（sample_global2.py）はどうでしょうか。

sample_global2.py

a = 10

def func1():
    print(f'a={a} in func1')

print(f'a={a} in initial')

func1()

実行してみると、

a=10 in initial
a=10 in func1

という結果が返ってきます。計算数学演習でC言語を学習した人は「何故関数func1の中で、引数でもなく、定義もされていない変数aの値が得られるのだろうか？」という疑問を持つ人がいるかもしれません（疑問に思ってくれると助かります）。

※ C言語では、main関数で定義した変数の値を別の関数で使うには、引数として使うことを定義しなければならない、と教えました。一方、Pythonでは、このmain関数というものがありません（あるいは、1つのpythonプログラムそのものがmain関数という考え方もできますが）。

今回の変数aのように、関数の定義の外に定義している変数をグローバル変数といい、定義されている関数内でも参照することができます。では次のプログラム（sample_global3.py）はどうでしょうか？

sample_global3.py

a = 10

def func1():
    print(f'a={a} in func1')
    a = 20

print(f'a={a} in initial')
func1()
print(f'a={a} in final')

もしグローバル変数である変数aを関数内で読み書きすることできるのであれば、

a=10 in initial
a=10 in func1
a=20 in final

という結果が返ってくることを期待するはずですが、残念ながらエラーが返ってくるはずです。これは、「関数内で変数aの値の代入（書き換え）を伴っているため、変数aはグローバル変数ではない」と解釈されてしまうからです。このときの変数aは、ローカル変数といい、関数内でのみ利用可能な変数になってしまっています。そのため、sample_global3.pyのprint(f"a={a} in func1")の段階で、「変数aが定義されていないのに使われているよ！」とエラーがでるわけです。実際、次のプログラム（sample_global4.py）を実行すると、関数内での処理が終わった後は、さも代入などなかったかのように元の値に戻っています（正確にはそもそも変わっていないのですが…）。

sample_global4.py

a = 10

def func1():
    a = 20
    print(f'a={a} in func1')

print(f'a={a} in initial')
func1()
print(f'a={a} in final')

a=10 in initial
a=20 in func1
a=10 in final

では、グローバル変数として関数内で変数aの値を変更するにはどうすれば良いのでしょうか？例えば次のプログラムのように関数内で変数aがグローバル変数であることを明示すれば良いでしょう。

sample_global5.py

a = 10

def func1():
    global a
    a = 20
    print(f'a={a} in func1')

print(f'a={a} in initial')
func1()
print(f'a={a} in final')

global aによって、変数aがグローバル変数であることが宣言されています。

これらのプログラムを図的に表すと次のようなイメージです。

このように、グローバル変数とローカル変数という考え方を上手く使うと、効率的なプログラムを書くことができます。例えば、今後行うシミュレーションにおいて、定数パラメータをグローバル変数で定義すれば、いちいち引数などを使わなくても、プログラム全体で使用することができます。一方、広い範囲で使える、ということは、どこかで変更すると、その変更が他の部分に渡って反映されてしまう、ということでもあります。この場合、プログラムがこちらの想定外のところで予想と反した挙動をとるかもしれませんので、使用には注意が必要です。また、sample_global4.pyのように、同じ名前のグローバル変数とローカル変数が混在するようなプログラムも極力避けたほうが良いでしょう。

ちなみにC言語においてもグローバル変数を定義することができます。ヘッダーファイルなどのincludeが終わったあと（関数の定義の前）に変数を定義すれば、そのプログラム全体に渡って使うことができます。

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

今回は次のような常微分方程式の初期値問題を、オイラー法と呼ばれる方法を用いて数値的に解いて結果をグラフにしてもらいます。オイラー法は計算数学演習で扱っていますが、おさらいをしておきます。次の問題を考えます。

(1) $\begin{equation*} \left\{ \begin{array}{ccc} \displaystyle\frac{dy}{dt} &=& f(t, y) \\ y(0) &=& a \end{array} \right. \end{equation*}$

求めたいのは、 $t>0$ の範囲の $y(t)$ であり、関数 $f(t, y)$ および定数 $a$ は与えられているとします。また、 $f(t, y)$ の値は $t, y$ が与えられればいつでも計算できるものとします。（関数として定義するとよい。）

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

差分商と差分方程式

$f(x)$ を $x$ の関数としたとき、 $f(x)$ の $x=a$ における微分係数は

(2) $\begin{equation*}\lim_{h\rightarrow 0} \frac{1}{h}\left\{ f(a+h) - f(a) \right\}\end{equation*}$

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

(3) $\begin{equation*}D = \frac{1}{h}\left\{ f(a+h) - f(a) \right\}\end{equation*}$

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

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

(4) $\begin{equation*}\frac{1}{h}\left\{ y(t+h) - y(t) \right\} = f(t, y(t))\end{equation*}$

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

(5) $\begin{equation*}\frac{1}{h}\left\{ Y(t+h) - Y(t) \right\} = f(t, Y(t))\end{equation*}$

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

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

(6) $\begin{equation*}Y(0) = a\end{equation*}$

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

(7) $\begin{equation*}Y(t+h) = Y(t) + hf(t, Y(t))\end{equation*}$

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

(8) $\begin{equation*}Y(h) = Y(0) + hf(0, Y(0))=a+hf(0, a)\end{equation*}$

となり、 $Y(0)$ から直ちに $Y(h)$ が求まります。同様にして（7）式を繰り返し用いると、

(9) $\begin{equation*} \left\{ \begin{array}{rll} Y(2h) &=& Y(h)+hf(h, Y(h)), \\ Y(3h) &=& Y(2h)+hf(2h, Y(2h)), \\ Y(4h) &=& Y(3h)+hf(3h, Y(3h)), \\ &\vdots& \\ Y((j+1)h) &=& Y(jh)+hf(jh, Y(jh)), \\ &\vdots& \\ \end{array} \right. \end{equation*}$

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

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

(10) $\begin{equation*} \left\{ \begin{array}{rll} \frac{1}{h}\left(Y_{j+1}-Y_{j}\right) &=& f(t_j, Y_j), \\ Y_0 &=& a \end{array} \right. \end{equation*}$

となり、結局（1）式の常微分方程式の初期値問題が、 $Y_j$ に関する漸化式の問題（10）式に置き換えられました。この方法をオイラー法と呼びます。見やすいように（10）を書き直すと

(11) $\begin{equation*} \left\{ \begin{array}{rll} Y_0 &=& a, \\ Y_{j+1} &=&Y_{j} + hf(t_j, Y_j) \end{array} \right. \end{equation*}$

となります。

1階常微分方程式に対するオイラー法のアルゴリズム

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

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

N, j は整数型(int)

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

T, N, aを設定
配列Y[N+1], t[N+1]を用意する
h = T/N　（TおよびN双方が整数型でないか注意）
Y[0] = a　（初期条件の設定）
t[0] = 0
j = 0, 1, 2, ..., N-1の順に（for文）
    t[j+1] = (j+1) * h
    Y[j+1] = Y[j] + h * f(t[j], Y[j])
以上をj=N-1まで繰り返す。
matplotlibを使ってグラフを描く

f(t,Y) は関数として定義するとよいでしょう。ここではfor文の中でtの具体的な値を代入していますが、for文の外でt = np.linspace(0,T,N+1)としてももちろんOKです。

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

さて、実際にオイラー法をつかって常微分方程式を解いてみましょう。今回も、計算数理Aの講義に出てきた次の常微分方程式を用いて演習を進めます。(ただし講義では $N$ としていた変数名を $y$ と書き換えました。上記アルゴリズム等の分割数 $N$ とまぎらわしいので。)

(P) $\begin{equation*} \left\{\begin{array}{lll} \frac{d}{dt} y(t) &=& - \lambda y, \\ y(0) &=& y_0 \end{array}\right. \end{equation*}$

前回（第3・4回）は、この微分方程式を解析的に解いて得られた解をグラフ化してもらいました。今回は(P)式をオイラー法を用いて数値計算で解いて、得られた数値解をグラフ化してもらいます。

課題1

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

$Y_0 = 1, \lambda = 1, 0 \le t \le 10$

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

先週、解析的に求めた解をグラフとして表示したものは以下であるが、殆ど違いがない事に注意。（課題2でみるが、微妙に異なっている）

課題2

前回解析的に求めた解(S)とオイラー法で求めた数値解を重ねて描くプログラムを作成せよ。例えば、次のようになる。この例では、解(S)を青色で描き、数値解を赤色で描いている。（この例では、0≦t≦10 を 20等分している）

刻み幅を小さくすると両者の差は縮まる。（次の例では、0≦t≦10 を 300等分している）

課題3

次の常微分方程式の初期値問題をオイラー法で解き、t と y の関係をグラフ化せよ。（次の小ネタも参照）

(1) dy/dt = -y + sin(t),　y(0) = 1 　(時刻 0≦ t ≦20 の範囲を20000等分割、100ステップ毎に値を表示)

(2) dy/dt = y(5-y),　y(0) = 1 　(時刻 0≦ t ≦20 の範囲を20000等分割、50ステップ毎に値を表示)

(3) dy/dt = y – 2y³,　y(0) = 0.1 もしくは　y(0) = -0.1　(時刻 0≦ t ≦20 の範囲を20000等分割、100ステップ毎に値を表示)

小ネタ：配列の要素を部分的に抜き出す

pythonの配列の便利さとして、次のような使い方ができます。

x = np.linspace(0,1,101)
x[[0,4,6]] # 配列の0・4・6番目を抜き出す
x[3:10:2] # 配列のうち3以上10未満の番号のものを2つ飛ばしで抜き出す
x[:4] # 配列のうち4未満の番号のものを抜き出す
x[4:] # 配列のうち4以上の番号のものを抜き出す
x[4:-1] # 配列のうち4以上の番号でかつ最後でないものを抜き出す

これを使うと、必要以上にデータを描画したくないときに役に立ちます。

レポート課題

課題2で作成したプログラムを少し変更して、「λの値の入力を促し、そのλの値における解析解とオイラー法で求めた数値解を重ねて表示する」プログラムを〇〇.py形式（名前は自由）で作成し、Bb9から提出して下さい。初期値、x、yの表示範囲は課題2と同じでかまいません。 0≦t≦10 を 100等分とすること。数字や軸の名前が描いてあるかどうかも評価のポイントです。関数を作ってmain関数をシンプルに（読み易く）記述しているものは高く評価します。コメントは必ず行うこと（#以降はコメント文になります）。また、プログラムの最初の方にコメントとして学生番号と氏名は（最低限）必ず記述すること。（# bxxxxxx 氏名） 提出期限2020年5月29日（金）18時00分

自由課題(上記の課題が終わった人へ）

セルオートマトン

一次元セルオートマトンについて調べて描いてみてください。ただし、ある時刻jのときの各セルiがもつ値x[i][j]（0か1）は自分の値と両隣の値をもとに、次の時刻j+1の自分の値を決定する物とします。（セルの数は自由に設定してください）

例としてRule90と呼ばれるもの以下に示します。（iがセルの番号で、jが時刻を表すと思ってください）

多次元配列を使っても良いですし、1次元配列を使っても良いです。1のところだけ●をつけると、次のような図が得られるはずです。

今回のように、（自分自身と）両隣のセルの値に応じて次のステップの値が決まるようなものを1次元3近傍系のセルオートマトンと言います。取りうる値の組合せは、3つのセルが2状態（0 or 1）を取るので、通りになります。それぞれの値の取り方によって、以下のような決まりによってRuleの番号が定められます。ちなみに藤井が個人的に好きなRuleは184です。（何故でしょう？）

講義で出てきた１階常微分方程式

いくつか講義で出てきた問題をあげておきます。オイラー法で解いてみて、解をグラフとして表示してみてください。初期値やパラメータを適当にきめて、いくつかシミュレーションしてみてください。

(12) $\begin{equation*} \left\{ \begin{array}{rcl} \frac{dx}{dt} &=& \gamma x \\ x(0) &=& x_0 \end{array} \right. \end{equation*}$

(13) $\begin{equation*} \left\{ \begin{array}{rcl} \frac{du}{dt} &=& u(1-u) \\ u(0) &=& u_0 > 0 \end{array} \right. \end{equation*}$

２階常微分方程式の数値解法

２階常微分方程式もオイラー法を用いて解くことができます。次の問題を考えて行きましょう。

(14) $\begin{equation*} \left\{ \begin{array}{rcl} y'' &=& -\omega^2 y \\ y(0) &=& 1 \\ y'(0) &=& 0 \end{array} \right. \end{equation*}$

講義の説明の通り(下の枠中)、これは単振動を記述する２階の常微分方程式です。ただし簡単の為、時間 $t$ での微分を ‘ を用いてあらわしています。バネ定数 $k$ のバネの一端を固定し、他端に質量 $m$ のおもりをつけた状況を考えます。

このおもりを少しだけ右に引っ張ってから手を放します。おもりの中心の時刻 $t$ におけるつり合いの位置からのずれを $y(t)$ とすると、 $y(t)$ が満たす微分方程式は

(15) $\begin{equation*} \frac{d^2 y}{dt^2} = - \omega^2 y \end{equation*}$

ただし、

(16) $\begin{equation*} \omega = \sqrt{k/m} \end{equation*}$

としました。

２階の常微分方程式は次のように連立の１階の常微分方程式に帰着できます。 $y_1(t)=y(t), y_2(t)=y'(t)$ とおくと、式(14)は、

(17) $\begin{equation*} \left\{ \begin{array}{rcl} y_1' &=& y_2 \\ y_2' &=& -\omega^2 y_1 \\ y_1(0) &=& 1 \\ y_2(0) &=& 0 \end{array} \right. \end{equation*}$

となります。この形になれば、前回紹介したオイラー法を用いて解を求めることができます。アルゴリズムは次のようになります。但し、 $T$ を解を求める時間の上限、 $h$ を時間刻み幅、 $N$ を時間刻み数とします。また、変数は $y_1, y_2$ として、 $a_1, a_2$ を初期値とします。また、 $f_1(t, y_1, y_2)$ および $f_2(t, y_1, y_2)$ はそれぞれ関数で式(17)の場合には $f_1$ が第一式の右辺、 $f_2$ が第二式の右辺となります。

2階常微分方程式に対するオイラー法のアルゴリズム

解を求める時刻の上限を T とし、適当な自然数 N を定めて h=T/N とする。 T, h, a1, a2, y1, y2, y1_new, y2_new, t 及び関数 f1, f2 は実数型 N, j は整数型

T,N,h を定める
配列y1[N+1], y2[N+1], t[N+1]を用意する
f1(t,y1,y2), f2(t,y1,y2) を関数として定義
y1[0] = a1, y2[0] = a2, t[0] = 0
j = 0,1,2,....,N-1の順に
    t[j+1] = (j+1)*h
    y1[j+1] = y1[j] + h*f1(t[j], y1[j], y2[j])
    y2[j+1] = y2[j] + h*f2(t[j], y1[j], y2[j])
を繰り返す
matplotlibでグラフを描画

課題4

式(14)を解析的に解け。

課題5

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

色々な h に対する計算結果は次のようになる（青はh=0.1、オレンジはh=0.01、緑はh=0.001、赤はh=0.0001のときの結果。紫は解析解のグラフ）。この結果から、h は十分小さい値でないとまずいことがよく分かる。

オイラー法では h を十分小さくとらないと本当の解をうまく近似できません。h を小さく取るということは、計算時間が多くかかることを意味します。そこでより高精度で大きな h でもうまく解を近似できるような解法が求められます。そのような要求に応える解法としては、ルンゲ・クッタ法があります。（計算数学演習でもやっていますが、次回改めて簡単に紹介します。）

課題6

課題5のプログラムを用いて、

(18) $\begin{equation*} E = \frac{1}{2}\left( y' \right)^2 + \frac{1}{2} y^2 \end{equation*}$

を計算せよ。結果は、横軸を $t$ 、縦軸を $E$ とするグラフで確認すること。この $E$ は $m=k=1$ のときの、バネの力学エネルギーに対応する。

課題7

実際のバネではバネの復元力以外におもりの速度 $y'$ に比例する抵抗が働く。その抵抗を考慮すると、式(14)ではなく次の式になる。( $2\gamma y'$ の項が新たに加わった。)

(19) $\begin{equation*} \left\{ \begin{array}{rcl} y'' + 2\gamma y' + \omega^2 y &=& 0 \\ y(0) &=& 1 \\ y'(0) &=& 0 \end{array} \right. \end{equation*}$

ここで $\gamma$ は抵抗力の強さを表す正の定数です。式(19)の問題を連立１階常微分方程式に変換し、オイラー法を用いて解き、下のように「左端から線の繋がった赤い質点が振動する」アニメーションを作成せよ。

ただし、 $\omega = 2, \gamma = 0.2$ とし、時間の刻み幅 $dt$ は、0.01とする。また、左上の”t= ..”は現在の $t$ を表示している。デモプログラムでは、時間の刻み幅は $dt=0.01$ として、 $t=50$ まで計算している。

今回は初めて微分方程式を解いてアニメーションを作成する、ということをやります。前回の資料も参考にしつつ、可能な限り自力で頑張ってみてください。どうしても無理！という人は、以下のリンク先に解答例を置いておきますので、参考にしてみてください（文字化けする場合は、テキストエンコーディングをUnicodeにするか、ダウンロードしてエディタで見てください）。

解答例1：逐次描画する場合（推奨、簡単）
解答例2：animation.ArtistAnimation関数を使う方法（ややコツが必要）
解答例3：animation.FuncAnimation関数を使う方法（トリッキーさが必要）

今後もこのようなアニメーションを作成するプログラムを作る機会が増え、レポートでも提出してもらうことが増えます。提出するプログラムでは、アニメーションの描画をが動く毎におこなうようにしてください。（ヒント１： "%"をつかったif文） これは、無駄にプログラムの実行時間が長くなるのを防ぐためです。

"t="の部分のアニメーションのさせ方：
他の点や棒の表記と同様にオブジェクトとしてテキストを作り、それをset_text関数を使って更新していきます。簡単な例は以下のようになります。

注意点はax.plotではline, のようにカンマをつけていましたが、textにはいりませんのでよくよく注意してください。

%matplotlib tk
import numpy
import matplotlib.pyplot as plt

T = 100
h = 0.01

fig = plt.figure()
ax = fig.add_subplot(111)

t = 0.0
txt = ax.text(0.3, 0.5, f"t={t:.2f}", transform=ax.transAxes)

for i in range(T):
    t = i*h
    txt.set_text(f"t={t:.2f}")
    plt.pause(0.01)

課題8

下図のような重力下で棒の長さ $L$ が一定の振り子運動について考える。振り子の角度(ラジアン）を $X$ とすると $X$ がある程度小さければ、 $X$ は

(20) $\begin{equation*} LX'' = -g X \end{equation*}$

という運動方程式を満たす（ $g$ は重力加速度、摩擦は無視しています）。 $L=1, g=9.8$ として、この方程式を適当な初期条件、 $X(0) = 0, X'(0) = a$ でシミュレーションし振り子の動画を作成せよ。

今回の課題は、オイラー法では十分に小さい $h$ を用いないと計算結果がすぐにずれてしまいますが、練習としてやってみてください。

次回、オイラー法より計算精度の高いルンゲ・クッタ法を紹介しますが、そのときに比較してみてください。

課題9

課題8の振り子が２つ相互作用する場合を考える. 例えば, ２つの振り子の角度を $2\pi X_1, 2\pi X_2$ とすると、 $X_1, X_2$ が

(21) $\begin{equation*} \left\{ \begin{array}{rcl} L_1X_1'' &=& - g \sin(2\pi X_1) - C \sin(2\pi(X_1 - X_2)) \\ L_2X_2'' &=& - g \sin(2\pi X_2) - C \sin(2\pi(X_2 - X_1)) \end{array} \right. \end{equation*}$

という運動方程式を満たすとする。この方程式を適当な $L_1, L_2, C$ と初期条件、 $X_1(0), X_1'(0), X_2(0), X_2'(0)$ でシミュレーションし、振り子のアニメーションを作成せよ. また振り子の数をもっと多くした場合も考察せよ。

Masashi FUJII

計算数理A演習第5・6回