第2回
matplotlib で折れ線を描く、線形補間
（2021/03/02 最終修正）
簡単な常微分方程式とその解のグラフ化、線形補間も
次の簡単な常微分方程式を用いて演習を進めます。この方程式は、例えば放射性同位体の崩壊を記述するモデル方程式として導出されます。例えば放射性炭素による年代測定（第3回で解説）にはこういった微分方程式が役立ちます。N(t) を時刻 t における放射性同位体の原子数とすると下記のような微分方程式が得られます。
(P)
さて、この微分方程式は簡単に解けて、一般解は次のようになります。
N(t) = N0 e-λt (S)
最初の目標は、解（S）をグラフとして表示するプログラムを作成することです。パラメータλおよびN0を与えれば、（S）のグラフをかけるはずです。コンピュータを用いてこのようなグラフを描く場合には折れ線近似で描く必要があります。
課題1
１階常微分方程式の初期値問題(P)を解け。（微分方程式の教科書等を参照）
準備1
Python でグラフを描くために、matplotlibを用います（マニュアル）。また、数学関数を用いるので、mathモジュールも利用します。その為、Pythonコードのはじめに、以下のコードが必要となります。
import math
import matplotlib.pyplot as plt
課題2（演習）
matplotlibを用いて、
N(t) = N0 e-λt (S)
のグラフを描け。ただし、
N0 = 1, λ = 1, 0 ≦ t ≦ 10
とする。（この例では、0≦t≦10 を 20等分している。つまり時間刻み幅を dt と書くことにすると、 dt = 10.0/20.0 = 0.5 である。曲線を20本の直線のつながりである折れ線で近似するのである。）
Pythonのmathモジュールでは，ea は math.exp(a) で求まる。また、matplotlibのlabelにTeX記法が使える。例えば、以下の様。
label="$N(t)=\mathrm{e}^{- \lambda t}$"
Code2-2.py : グラフの描画
# Code2-2.py
import math
import matplotlib.pyplot as plt
N0 = 1.0
Lambda = 1.0
T = 10.0
N = 20
dt = T/float(N)
t, y = [0.0]*(N + 1), [0.0]*(N + 1)
for i in range(0, N + 1):
t[i] = dt*i
y[i] = N0*math.exp(-Lambda*t[i])
plt.title("2-2")
plt.plot(t, y, color="red", linewidth=2, marker="o", label="$N(t)=\mathrm{e}^{- \lambda t}$")
plt.xlabel("t", fontsize=20)
plt.ylabel("N(t)", fontsize=20)
plt.legend()
plt.show()
#plt.savefig("fig2-2.png")
課題3
N0=1, λ=1のとき，N(τ) = 0.5 となる τ を求めよ。（(S)を用いて具体的に求めよ）τは半減期とも呼ばれる。
答え： log2 (約 0.693)
Pythonのmathモジュールでは log2 は math. log(2.0) で求まる。
課題4（演習）
課題2で作成したプログラムを用いて課題3のτを近似的に次のような手法で求めたい。そのようなプログラムを作成せよ。（課題2でわかっていることを何故近似的な方法を使って求めるのか？この場合はそうであるが、例えば飛び飛びの時間に得られた実データについてτを求めたい場合がある。そのように想像すべきである。）
方法
N(t) の値を（S)を用いて飛び飛びの t の値 ti = dt*i （ただし dt は時間刻み幅）について求めるのであるが、(S)と課題2のグラフから明らかなように N(t) は単調減少関数である。であるから、
N(ti) ≧ 0.5 > N(ti+1)
となる i がただ一つあるであろう。このとき、 ti を τ の近似値 τ* とする。
小問１：τ* を求めよ．
小問２：
0 ≦ t ≦ 10 を 100等分して上記方法にてτ* を求め、課題3で求めた値との差（誤差と呼ぶ）
err* = |τ- τ*|
を求めよ。
注意：当然ながら刻み数を 100 から 1000 に増やせばより正確な値が求まる。
小問３：
err* は刻み数を変更するとどのようにかわるか？刻み数100, 1000, 10000に関してそれぞれ調べ，横軸を刻み数，縦軸をその刻み数に対するerr* として，グラフにしてみよ(例えばExcelを用いる)．
Code2-4.py : 条件を満たす区間を見つける
# Code2-4.py
import math
import matplotlib.pyplot as plt
N0 = 1.0
Lambda = 1.0
T = 10.0
errs = []
NN = [100, 1000, 10000]
for N in NN:
dt = T/float(N)
t, y = [0.0]*(N + 1), [0.0]*(N + 1)
for i in range(0, N + 1):
t[i] = dt*i
y[i] = N0*math.exp(-Lambda*t[i])
itau = 0
for i in range(0, N):
if y[i] >= 0.5 and y[i + 1] < 0.5:
itau = i
taus = itau*dt
errs.append(abs(math.log(2.0) - taus))
print(taus, N, errs)
plt.title("2-4")
plt.plot(NN, errs, color="red", linewidth=2, marker="o", label="err*")
plt.yscale("log")
plt.xscale("log")
plt.xlabel("N", fontsize=20)
plt.ylabel("err", fontsize=20)
plt.legend()
plt.show()
#plt.savefig("fig2-4.png")
課題5（演習）
課題4の方法ではあまり良い精度で近似値を求めることができない。次の方法でより精度良く近似値をもとめることを考える。課題4と同様に，
N(ti) ≧ 0.5 > N(ti+1)
なる i を求め、２点 (ti, N(ti)), (ti+1, N(ti+1)) を結ぶ直線を考える（下図の青丸および青線）．その直線と N = 0.5 の直線の交点の t の値 τ** をもって、τ の近似値をする（下図赤丸）．このような手法を線形補間と呼ぶ．下図において，黒い曲線が例えばN(t)であって，黄緑色の丸がτである（ただし，このグラフは適当に描いたものである）．
小問１：τ** を求めよ．
小問２：
0 ≦ t ≦ 10 を 100等分して上記方法にて τ** を求め、課題3で求めた値との差
err** = |τ- τ**|
を求めよ。課題4の結果よりも良いだろうか？
注意：
課題4の方法でも刻み数を増やせばより正確な値が求まるのであるが、課題5の補間を使えば、あまり多くない刻み数でも正確な近似値が得られることを実感してほしい。
小問３：
err**は刻み数を変更するとどのようにかわるか？刻み数100, 1000, 10000に関してそれぞれ調べ，横軸を刻み数，縦軸をその刻み数に対するerr**として，グラフにしてみよ(例えばExcelを用いる)．
小問４：
課題4の小問３および，課題5の小問３にて得られたグラフを重ねて表示してみよ(例えばExcel を用いる）．
Code2-5.py : 線形補間を用いる
# Code2-5.py
import math
import matplotlib.pyplot as plt
N0 = 1.0
Lambda = 1.0
T = 10.0
errs, errss = [], []
NN = [100, 1000, 10000]
for N in NN:
dt = T/float(N)
t, y = [0.0]*(N + 1), [0.0]*(N + 1)
for i in range(0, N + 1):
t[i] = dt*i
y[i] = N0*math.exp(-Lambda*t[i])
itau = 0
for i in range(0, N):
if y[i] >= 0.5 and y[i + 1] < 0.5:
itau = i
taus = itau*dt
errs.append(abs(math.log(2.0) - taus))
dy = y[itau + 1] - y[itau]
tauss = dt*(0.5 - y[itau])/dy + dt*itau
errss.append(abs(math.log(2.0) - tauss))
#print(taus, N, errs)
plt.title("2-5")
plt.plot(NN, errs, color="red", linewidth=2, marker="o", label="err*")
plt.plot(NN, errss, color="blue", linewidth=2, marker="o", label="err**")
plt.yscale("log")
plt.xscale("log")
plt.xlabel("N", fontsize=20)
plt.ylabel("err", fontsize=20)
plt.legend()
plt.show()
#plt.savefig("fig2-5.png")
課題6
課題5において刻み数100で近似値を求めた場合の誤差と同程度に小さい誤差とするためには課題4において刻み数をどの程度にしなければならないか。実験してみよ。
注意：
課題4、課題5の方法は、 N(t) が具体的に関数として与えられていない場合（例えば実測データなど）にも有効であることに注意。
課題4、課題5は２分法の離散版とみることもできます。N(t)の具体的な関数形が分かっている場合（例えば今回の場合）には、２分法やニュートン法等の方法で根を求める事ができます。N(t)の具体的な関数系がわからず、単なる観察データである場合には、今回の手法が良い手段となります。
今回は、微分方程式を数値的に解いたわけではないことに注意してください。解析的にもとまった解 (S) をグラフにしたにすぎません。次回は、常微分方程式を数値的に解いてもらい、それをグラフ化してもらいます。
課題4，課題5の内容の発展的内容として次の問題も考えてみて下さい．
発展問題
簡単な発展問題であるが，実際に取り組むことが重要である．いくつか重要な事を学ぶことができる．（例えば，課題4, 5では，関数は単調減少関数であったが，これは違う為，交点部分の判定に工夫がいるであろう．）
課題5の方法を用いて，
y = sin(x*x)
と直線
y = 0.5
の交点の x 座標を 0 ≦ x ≦ πの範囲で全て求めよ．下のグラフに示すように（縦軸をy，横軸をx とした），交点は4つ存在する．
