Python を使って微分方程式を解いてみる
みなさんこんにちは。明けましておめでとうございます。
年末から結構な頻度でお酒を飲むことができて幸せな佐野です。
新年一発目の記事ですが、前回に引き続き、Python を用いてやっていきます。
今回は前回とはタイトルの趣向を変えて、「微分方程式 を解いてみる」としました。
さてみなさん、 微分方程式 と言うものをご存じでしょうか？方程式と名が付いている以上、方程式であることは当然なのですが、「なんか特殊な方程式なんだろうなあ」と言う方がほとんどだと思います。
以下で本記事に関する数学的および物理学的知識を説明しますが、本質的ではないため、極力単純に書いていこうと思います。そうなっていなかったらすみません。
方程式の種類
方程式と一言で言っても、実はいくつかの種類があります。大まかに分類分けを行い、以下にまとめます。
代数方程式
ある未知数に対する関係式のこと
一般的に「方程式」と呼ばれるものの正式名称-
関数方程式
ある関数に対する関係式のこと
微分方程式は 関数方程式 の方に分類分けされます。
連続的な値と離散的な値
ここでは 連続的である と言う概念と、離散的である と言う概念に関して説明します。
連続とは、おそらく簡単にイメージできると思うのですが、切れ目がなく連なって続くことを言います。離散的もその言葉からイメージできると思いますが、離れて散っていることを言います。
例えば、0 と 1 という値の間には 0.1 とか 0.2 とか 0.8 とかいう値が存在します。では、0.1 と 0.2 という値の間にはどの様な値があるでしょうか？そうです、みなさんの想像通り、0.11 とか 0.12 とか 0.15 とかいう値が存在しますよね。
この様にどんどん突き詰めていくと、 0.1 と 0.2 と言う値の間にでさえ無限に値が存在します。この様な状態にある値を 連続的な値 と言います。逆に、0.1 と 0.2 の間に値が存在しない、0.1 の次の値は 0.2 であるといったような、飛び飛びの状態にある値を 離散的な値 と言います。
一般的に、数値計算等を行う際は基本的に離散的な値を用います。
特定の例外はあるものの、通常数学で考慮すべきは連続的な値なので、連続的な値で関係付けられている方程式を離散的な値でも扱える様に近似しなければいけません。代数方程式であれば元の式とほとんど変化がないのですが、関数方程式や微分方程式となると話は変わってきます。その際の手法として、オイラー法とルンゲ=クッタ法と言う近似手法をこの記事では使用しますが、この記事ではそれらの説明までは行いません。
では早速解いて行きましょう。
微分方程式を解く
例題.1 一階微分方程式を解く
一階微分方程式と言うのは、 1回のみ微分された関数を含む方程式である と言う認識でここでは大丈夫です。決して誤字ではありません。
以下の様な一階微分方程式を解いてみましょう。
x は t の関数です。右側にある x(0) = 1 というのは初期条件と呼ばれるものです。ここではオイラー法という近似手法を使用します。計算プログラムは以下になります。
import math
import matplotlib.pyplot as plt
# 陽的オイラー法
x = 1
dt = 0.1
xAxis1 = []
tAxis1 = []
# 時間発展
for i in range(20):
xAxis1.append(x)
t = i * dt
x = (1 + dt) * x
tAxis1.append(t)
# 陰的オイラー法
x = 1
dt = 0.1
xAxis2 = []
tAxis2 = []
# 時間発展
for i in range(20):
xAxis2.append(x)
t = i * dt
x = x/(1-dt)
tAxis2.append(t)
# 厳密解
x = 1
dt = 0.1
xAxis3 = []
tAxis3 = []
# 時間発展
for i in range(20):
t = i * dt
x = math.e**t # 厳密解 : x = e^t
xAxis3.append(x)
tAxis3.append(t)
# プロット
plt.plot(tAxis1, xAxis1, '-o', label='exp: Euler(+)')
plt.plot(tAxis2, xAxis2, '-o', label='exp: Euler(-)')
plt.plot(tAxis3, xAxis3, label='exp: Exact')
# 判例の表示
plt.legend()
# 罫線の表示
plt.grid()
# タイトル、軸名の設定
plt.title('practice')
plt.xlabel('t')
plt.ylabel('x(t)')
# プロットの表示
plt.show()
基本的にはまず初期値を設定し、x の時間発展を計算している、という感じになります。
厳密解とは、オイラー法を用いずに手で微分方程式を解いた時の解になります。近似を用いないので厳密な解ということ厳密解ですね。
計算結果をプロットしたのが以下のグラフになります。
オイラー法を用いた二つのグラフは厳密解に対してずれているのが分かると思います。
近似とは言え、結構ずれるものだなあと感じました。
では、もう一つ微分方程式を解いてみましょう。
例題2. 二階微分方程式を解く
以下の様な二階微分方程式を解きます。
この方程式が言っていることを図示すると、以下の様な図になります。
おそらくみたことがある人もいるのではないでしょうか？これは物理の中でも結構有名な問題で、いわゆる 単振動の運動方程式 と呼ばれる問題です。壁に取り付けられたバネの先端に質量 m の重りが取り付けられているというものです。重りを引っ張ると、バネが元に戻ろうとする力で重りが壁側に引き寄せられます。さらに、ある程度壁側に引き寄せられると、縮まったバネが元の長さに戻ろうとして、今度は重りが壁から引き離される方向に力が働きます。この時の重りの軌跡を計算してみましょう。ただし、バネ定数を k とし、抵抗は考えないものとします。
ここではルンゲ=クッタ法を用いて計算します。
import math
import numpy as np
import matplotlib.pyplot as plt
# 定数の設定
k = 1.0
m = 1.0
# ベクトル解を求める関数
def vectorSol(t, x):
return np.array([x[1], -k/m * x[0]])
# ルンゲ=クッタ法で求める単振動の解
xAxis = []
tAxis = []
# 初期条件
x = np.array([1.0, 0.0])
dt = 0.1
for i in range(126):
xAxis.append(x[0])
t = i * dt
k1 = vectorSol(t, x)
k2 = vectorSol(t + dt/2, x + k1 * dt/2)
k3 = vectorSol(t + dt/2, x + k2 * dt/2)
k4 = vectorSol(t + dt, x + k3 * dt)
x += (k1 + 2 * k2 + 2 * k3 + k4) * dt/6
tAxis.append(t)
# 厳密解
tAxis2 = []
xAxis2 = []
for i in range(0, 126):
xAxis2.append(math.cos(i * 4 * math.pi/(126)))
tAxis2.append(i * dt)
# プロット
plt.plot(tAxis, xAxis, '-o', label='cos(ωt)')
plt.plot(tAxis2, xAxis2, label='exp: Exact')
# 判例の表示
plt.legend()
# 罫線の表示
plt.grid()
# タイトル、軸名の設定
plt.title('Harmonic Oscillation')
plt.xlabel('t')
plt.ylabel('x(t)')
# プロットの表示
plt.show()
ここではベクトルを用いて計算を行いました。これは、 v = dx/dt と置くと、以下の様に表すことが出来るためです。
近似の手法は異なるものの、計算の流れは例題 1 と同様ですね。
計算結果をプロットしたのが以下になります。
綺麗な周期的な動きをするということがわかりますね。
ルンゲ=クッタ法で求めた近似解が厳密解とほぼ同じであることもわかりますね。オイラー法よりも精度が良いことが見て取れると思います。
まとめ
いかがでしたでしょうか。
個人的には、学生時代ぶりの微分方程式に触れることができてとても楽しかったです。反面、説明が多くなってしまった印象です、そこはすみません。これを機にみなさんも数学と物理を再度学んでいただけたら嬉しい限りです。
また機会があれば別の方程式も解いてみたいですね。
それではみなさん、本年も何卒よろしくお願いします。
参考
