無意味なブログ

勉強のこと、趣味のこと、日記など適当に

MENU

微分方程式をプログラムで解こう!(2):差分法[1]「sin(x)を微分する」

前回(第1回)
微分方程式をプログラムで解こう!(1):導入 - 無意味なブログ

今回は"差分法"について説明します。
基礎的な3つの差分法を用いてsin(x)を微分していきます。

差分法とは

差分法は、微分方程式を解くために微分を有限差分近似(差分商)で置き換えて得られる差分方程式で近似するという離散化手法を用いる数値解法である。
差分法 - Wikipedia

微分方程式をプログラムで解くには微分値を計算する必要がでてきます。そりゃそうですよね。"微分"方程式なんだから。

例えば関数f(x)=sin(x)について、関数f(x)をxで微分するとf'(x)=cos(x)となります。
sin(x)みたいな簡単な式で表される関数であればこうやって解析的に解くこともできますね。
ただ、コンピュータの世界では以下のグラフに示されるような、とびとびの値たちから微分値を計算しなくちゃいけないことがよくあります。
f:id:kouya17:20170121093754p:plain:w350
そういう時に使われる方法の1つとして「差分法」があります。

今回使う3つの差分法

今回は前進差分、後退差分、中心差分を使います。それぞれ式で表すと以下の通りです。

前進差分 {f'(x_0) \simeq \frac{f(x_0 +h)-f(x_0)}{h}}

後退差分 {f'(x_0) \simeq \frac{f(x_0 )-f(x_0-h)}{h}}

中心差分 {f'(x_0) \simeq \frac{f(x_0 +h)-f(x_0 -h)}{2h}}

式は以下のサイトを参考にしました。
https://www.sci.hokudai.ac.jp/~inaz/lecture/butsurisuugaku2/html/model/node41.html

とりあえずsin(x)を差分法使って微分してみる

3つの差分法を使ってsin(x)を微分してみましょう。また、解析解cos(x)と比較して、誤差も出してみます。

今回作成したfortranプログラムを以下に載せます。このプログラムをコンパイル・実行すると、結果ファイル"centd.dat","ford.dat","backd.dat"が作られます。

結果

区間0から2πについて、分割数nを10として計算した結果を示します。
それぞれグラフの紫線が正しい解(cos(x))、緑線が差分によって得られた数値解です。

前進差分の結果(ford.dat)

f:id:kouya17:20170122000221p:plain:w400
正しい解と比べてx軸方向負の向きにずれちゃってますね。

後退差分の結果(backd.dat)

f:id:kouya17:20170122000240p:plain:w400
今度はx軸方向正の方向に。

中心差分の結果(cend.dat)

f:id:kouya17:20170122000253p:plain:w400
上2つよりいい感じ。

誤差の比較

それぞれ解析解との絶対誤差を比較します。
紫線が前進差分による誤差、緑線が後退差分による誤差、水色線が中心差分による誤差です。
f:id:kouya17:20170122000026p:plain:w400
この中だと中心差分による誤差が一番小さいですね。

まとめ

3つの差分法を使ってsin(x)を微分してみました。
解析解との誤差を比較してみたところ、この3つの中では中心差分による誤差が一番小さくなりました。

シリーズ第2回「差分法[1]sin(x)を微分する」は以上です。
次回は分割数nを変化させて、結果がどう変わるのかみてみます。

広告を非表示にする