差分法って言葉、聞いたことありますか?
日常では多分聞かない言葉ですが、数値計算の分野では耳にする言葉です。
そして数値計算の分野では、有名な手法であり、私も中心差分法を使ってプログラミングをしてました。
今回は中心差分法について説明したいと思います。
微分と差分法
数値計算のプログラミングするということは何らかの式を入れると思います。
例えば、熱に関するなら、熱伝導方程式、流速ならナビエストークス方程式等、プログラムに組み込むと思います。
そして、紹介した式には微分が含まれています。
(大学で紹介される式には大体微分が含まれている気がする∂)
熱伝導方程式を例にすると、
ρc∂T/∂t = k(∂²T/∂x² +∂²T/∂y² + ∂²T/∂z²) + q′v
です。
そして残念ながら、上記の式をそのままプログラミング出来ません。
コンピュータは基本四則演算で表現しないといけないんです。
(パソコン…思ったほど優秀じゃないな)
上記の熱伝導方程式だと微分の部分を何とか四則演算で表現しないとプログラミング出来ないのです。
そこで使われる手法が差分法になります。
中心差分法
差分法にはいろいろな種類がありますが、今回は中心差分法を説明します。
簡潔に説明すると中心差分法は求めたい場所のひとつ前とひとつ後の値を使って、求めたい場所を決めます。
例えば、ƒ(i-1)に200℃、ƒ(i+1)に100℃という値を入れ、この2つの平均をf(i)に入れるという式ならf(i)は150℃になります。
実際に熱伝導方程式をプログラミングできるように変形していきます。
まず微分式を中心差分法(2次精度)を使って変形します。
dƒ/dx = (ƒ(i+1) – ƒ(i-1)) / 2Δx
d²ƒ/dx² = (ƒ(i+1) – 2ƒ(i) + ƒ(i-1)) / Δx²
分かりやすくするために定常計算、発熱なし、1次元の熱伝導方程式にします。
ρc∂T/∂t = k(∂²T/∂x² +∂²T/∂y² + ∂²T/∂z²) + q′v
∂T/∂t=0 , ∂²T/∂y² = 0 , ∂²T/∂z² = 0 , q′v = 0 より
0 = k(∂²T/∂x² ) … 一次元の熱伝導方程式
ここで中心差分法を適応すると
0 = k((T(i+1) – 2T(i) + T(i-1)) / Δx²)
T(i)を求めたいので
T(i) = 0.5 * (T(i+1)+ T(i-1))
これでプログラミングできる形になりました。
(do 文と組み合わせることでプログラミングできる)
境界条件の重要性
中心差分法について説明しましたが、何か疑問に思った人もいると思います。
それは一番端の値はどうやって求めるのか?
答えは求めることはできません。(中心差分法を使ってですが)
なので一番端の値はあらかじめ決めなければなりません。
それを境界条件といいます。
境界条件は重要なファクターであり、例えば、境界条件が1000℃なのか、500℃なのかで答えが変わってしまいます。
商用シミュレーションソフトを使ったことがある人は境界条件はよく考えて設定するようにと学んでるかもしれないですね。
それくらい重要です。
まとめ
中心差分法について説明してきましたが、他に片側差分法等もあります。
中心差分法と考え方は似ているので、中心差分法を理解したなら扱えると思います。
境界条件に関しては、商用シミュレーションソフトだろうが自作ソフトだろうがよく考えて設定しましょう。