前回と今回で「数値計算」シリーズ第3弾とて「常微分方程式の数値解法」を。 最後にごく簡単ながら、テキスト第7.5章に沿って、偏微分方程式についても触れておく。
常微分方程式の初期値問題その2
常微分方程式 \(\dfrac{d{x}}{d{t}}=f(t,x)\) を \({a}\le{t}\le{b}\) の範囲で、初期条件 \(x(t=a)=x_0\) のもとで解くことを考える。 解を求めたい範囲 \({a}\le{t}\le{b}\) を \(N\) 等分することにより、\(t=t_0,t_1,\cdots,t_N\) をとってあることにする。 点の間隔を \(h=\dfrac{b-a}{N}\)とかくと、 \begin{equation} t_i=a+{h}{i}=a+\dfrac{b-a}{N}i\quad(i=0,1,\cdots,N) \tag{10.10} \end{equation} 当然ながら、「小区間の幅」が小さいほうが計算の正確さは増すものの、実際の問題に適用するには正確さの増す「スピード」が重要になる。
第12回では手始めとして、最も簡単かつ基本的なオイラー法 \begin{equation} x_{i+1}=x_{i}+{f(t_i,x_i)}\times{h} \tag{10.11} \end{equation} をとりあげたのだが、
- オイラー法では、真の答からのずれが「小区間」の数の1乗に反比例して減少 する。
- この「スピード」では実用的な問題には向かない。 正確さを10倍増すには、「小区間」の数を10倍にする必要があり、その分だけ手間と時間が増える。 かといってその分だけ正確さが本当に増すとも限らない。
ホイン法
オイラー法では点 \((t_i,x_i)\) での導関数の値 \(f(t_i,x_i)\) だけを使って \(x_{i+1}\) を求めたが、ホイン法 (または修正オイラー法) では2つの点における導関数の値を用いて計算する。 具体的にはホイン法では、次の手順により \(x_{i+1}\) を求める。 \begin{equation} \begin{split} k_1&=f(t_i,x_i)\times{h} \\ k_2&=f(t_i+{h},x_i+k_1)\times{h}\\&=f(t_{i+1},x_i+k_1)\times{h} \\ x_{i+1}&=x_{i}+\dfrac{1}{2}(k_1+k_2) \end{split} \tag{10.15} \end{equation} ホイン法の意味は以下の通り。 図10.3の模式図も参照のこと。
- 点 \((t_i,x_i)\) における接線の傾きだけを使うのではなく、点 \((t_{i+1},x_{i+1})\) における傾きとの平均 をとってやれば、その区間での傾きがうまく表わされるはず。
- でも「点 \((t_{i+1},x_{i+1})\) における傾き」は分からない。
- \(t_{i+1}\) は分かるけれど、\(x_{i+1}\) の値は分からない。 (それも当然で、そもそも \(x_{i+1}\) を求めるための計算をしているのだから)
- そこで、「\((t_{i+1},x_{i}+k_1)\) における傾き」を代わりに使う ことにして、その区間での傾きを計算する。
中点法
中点法 (改良オイラー法) と呼ばれる方法では、ホイン法 (修正オイラー法) とはまた違った発想により、\((t_i,x_i)\) と \((t_{i+1},x_{i+1})\) の区間における傾きを表わそうとしたものである。 具体的には中点法では、次の手順により \(x_{i+1}\) を求める。 \begin{equation} \begin{split} k_1&=f(t_i,x_i)\times{h} \\ k_2&=f(t_i+{h}/2,x_i+k_1/2)\times{h}\\ x_{i+1}&=x_{i}+k_2 \end{split} \tag{10.16} \end{equation} 中点法の意味は以下の通り。 図10.4の模式図も参照のこと。
- 点 \((t_i,x_i)\) における接線の傾きを使うのではなく、その点と点 \((t_{i+1},x_{i+1})\) との中間の点における傾き をとってやれば、その区間での傾きがうまく表わされるはず。
- でも「その中間の点の場所」は分からない。
- \(t_{i+1}\) は分かるけれど、やはり \(x_{i+1}\) の値は分からない。
- そこで、「点 \((t_i,x_i)\) と \((t_{i+1},x_{i}+k_1)\) との中点」を代わりに使う ことにして、その区間での傾きを計算する。
ルンゲ-クッタ法
普通「ルンゲ-クッタ法」といえば以下の「古典的」ルンゲ-クッタ法のことを指す。 \begin{equation} \begin{split} k_1&=f(t_i,x_i)\times{h} \\ k_2&=f(t_i+\dfrac{1}{2}{h},x_i+\dfrac{1}{2}{k_1})\times{h} \\ k_3&=f(t_i+\dfrac{1}{2}{h},x_i+\dfrac{1}{2}{k_2})\times{h} \\ k_4&=f(t_i+{h},x_i+{k_3})\times{h} \\ x_{i+1}&=x_{i}+\dfrac{1}{6}({k_1}+2{k_2}+2{k_3}+{k_4}) \end{split} \tag{10.17} \end{equation} このルンゲ-クッタ法は、テイラー展開の最初の5つの項 (\(h^4\) まで) を正しく取り込んだ表式になっている。 多くの場合は、このルンゲ-クッタ法を使えばOK。
なお厳密にいえば「ルンゲ-クッタ法」とは、複数の点における導関数の値を計算しておき、その重みつき平均として \(x_{i+1}\) を求める方法の総称。 先に挙げた「ホイン法」も「中点法」もルンゲ-クッタ法の一種である。
偏微分方程式
独立変数が2つ以上ある微分方程式は「偏微分方程式」という。 その例として、熱伝導方程式をとりあげてみる。
長さ \(L\)、温度が \(T_1\) の棒がある。 時刻 \(t=0\) で棒の一方の端の温度を瞬間的に \(T=T_0\) に変化させ、\({t}>0\) の間ずっと維持するとする。 このとき、棒の中の温度の分布は時間とともにどう変化するか?棒の中の温度 \(T\) は、
- 経過した時間 \(t\)
- 棒の中の位置 \(z\) (例えば「一方の端からの距離」)
その他、物理現象のほとんどは偏微分方程式によって書かれる。
- 流体の運動を記述する方程式 (ナビエ-ストークス方程式)
- 波が伝わることを記述する方程式 (波動方程式)
- 地下の質量分布から重力ポテンシャルを求める方程式 (ポアソン方程式の一種)