非線形の微分方程式だって解けるケロ（差分の基礎編）

非線形の微分方程式だって解けるケロ

 

 

差分法を用いた、次の非線形微分方程式の初期値問題について考える。

　　

 

前進差分を用いて（１）を差分方程式に近似すると次のようになる。

　　

これから、

　　

という漸化式を得ることができ、x₀=0におけるyの値y₀を計算の出発点にし、前進的にの値を次々に求めることができる。

 

これに対して、後退差分を用いると、（１）の微分方程式に関する差分方程式は

　　

となり、次の２次方程式を得ることができる。

　　

この２次方程式を解くと、

　　

という漸化式が得られ、（３）式を用いて前進的にの値を次々に求めることができる。

 

（１）の場合、後退差分は２次方程式を解く必要があるんだケロよ。２次方程式を解くというひと手間が加わるだけではなく、（２）と（３）を比較すればわかるように、漸化式が複雑になる。

 

　　

だったら、３次方程式を解かないといけない。３次方程式ならば解の公式があるので代数的に解けないことはないけれど、n≧5のとき

　　

この式から得られる

　　

というのn（≧５）次方程式を代数的に解くことは一般にできないので、ニュートン法などを用いてを解く必要がある。

ということで、後退差分（陰解法？）は前進差分（陽解法？）とは違って難しいんだケロよ。だから、常微分方程式の初期値問題で後退差分を用いた解法（陰解法？）は一般に使われない。

陰解法が威力を発揮するのは、

　　

といった偏微分方程式の数値解法なんだケロよ。

 

ではありますが、Δx=0.1とし、（２）式と（３）式を使って、（１）の初期値問題を解いた結果は次のとおり。

 

微分方程式（１）の（厳密）解は簡単に求まりそうに見えるだろうが、これはリッカチ形の微分方程式で、この解は三角関数、指数関数、対数関数を用いて表すことができない、つまり、解けないんだケロ。

なので、比較参考のために、４次のルンゲ・クッタ法を用いてΔx=0.1で計算した値を厳密解のかわりにあわせて図示している。

この計算結果を見ると、陽解法による計算結果は厳密解よりも小さく、陰解法による数値計算結果は厳密解よりも大きいこと、そして、厳密解がこの間に位置していることがわかる。

ということで、陽解法による計算結果を陰解法による数値計算結果を足して２で割る、つまり、両者の平均をとると、厳密解とかなりよく一致することがわかる。

このことは、数値計算結果の絶対誤差を図示したグラフを見るとよくわかると思う。

 

 

前回に引き続き、陽解法と陰解法の計算結果を足して２で割ると、誤差のプラスマイナスが互いに消し合って、誤差が小さくなって、厳密解に近い値が出るようになるのであった。

Δx=0.1のときに、誤差がおよそ1/10になっており、また、陽解法、陰解法の誤差の程度がO((Δx)²)であるので、この解法の誤差はO((Δx)³)と推測される。だから、修正オイラー法、２次のルンゲ・クッタ法と誤差のオーダーと同じであるに違いない。

 

より進んだ解法は、

　　

①と②の両辺を足しあわせて、両辺を２で割ると、

　　

となるにゃ。

についての２次方程式（４）を解き、こうして新たに得られた漸化式を用いて計算してみるといい。

 

（４）を使った計算結果は、（２）と（３）を使って求めた値の平均値と同じ値になると思うかもしれないけれど、これが一致しないんだね〜。

（４）を使ったもののほうが精度はさらに向上するはずなんだよね〜。

嘘だと思うならば、の２次方程式（４）を解き、の漸化式を求め、その漸化式で計算してみな。

というか、お前ら、これくらいはやれよな。

 

何故、（４）式は、（２）や（３）より精度よく計算できるのかについては、尤もらしい説明をすることはできなくはない。

というのは、との中点をと表すことにし、それに対応するyの値をで表し、

　　

③と④は、それぞれ、

　　

となる。

⑤と⑥を足すと、が消えて、（４）式になる。

③は、Δxより小さいΔx/2という分割幅で仮想のの値を求めるいることを意味し、④式はこうして求められたを元にΔxよりも小さい分割幅でを求めることを意味している。

分割幅を小さくすると一般に計算精度が向上する。

また、仮想の点とそれに対応するyの値を考慮に入れると、（４）の左辺のはにおける中心差分とみなすこともできる。

したがって、（４）は、（２）や（３）よりも高精度に計算できる・・・。

 

 

感覚的にこう理解すればいいのではないか。

これにテーラー展開の話を交え、「（２）や（３）の誤差の程度がO(h²)であるのに対し（４）式はO(h³)で修正オイラー法、２次のルンゲ・クッタと誤差の程度が同オーダーである」なんて尤もらしい説明もできるけれど、かなり胡散臭い話だから、ここではしないにゃ。

 

（２）、（３）、（４）式は、０≦φ≦１である重みφを使うと、

　　

と、１つの式で表すことができる。

φ=0のときは陽解法、φ=1/2のときは半陰解法、φ=1のときは純陰解法と分類分けすることも可能。

 

そして、プログラムを作ってこの問題を解く場合、0≦φ≦1のすべてのφについて解けるものを作った方がいいにゃ。それが賢い方法だにゃ。