プログラミング言語

【2次精度の片側差分】 1階微分と2階微分を離散化する。

こんにちは(@t_kun_kamakiri)。

当ブログを読んでくれた方から質問をもらったので、記事にして回答したいと思います。

基本的に質問に対しては、お答えするスタンスです(‘◇’)ゞ

☟頂いた質問内容はこちら・・・

👆こんな感じの質問です。

カマキリ

どの記事を読んだ際の質問か?

記事はこちら

 

カマキリ

質問内容は?

「乱流の数値シミュレーション」のp.28(2.24)、(2.25)の導出方法がわからないから教えてほしいという内容です。
ちなみに「乱流の数値シミュレーション」とは☟こちらの書籍です。

乱流の数値シミュレーション 改訂版

乱流の数値シミュレーション 改訂版

梶島 岳夫
4,070円(09/28 02:30時点)
発売日: 2017/01/27
Amazonの情報を掲載しています

こちらの書籍は乱流モデルについて、標準的な難易度でわかりやすく書かれています(‘ω’)
宣伝はさておき・・・・

 

本日の記事の内容はこちら☟

記事の内容はこちら
「乱流の数値シミュレーション」

  • p.28(2.24)、(2.25)
  • p.29(2.28)、(2.29)

を導出する。

 

p.28(2.24)、(2.25)は何かというと・・・・☟これです。

\(\frac{df}{dx}\)のような一階微分を、等間隔格子での片側差分でどう表されるのかという話です。
結果を示しておくと・・・

(2.24)

\begin{align*}
f^{\prime}=\frac{f_{j-2}-4f_{j-1}+3f_{j}}{\Delta}+\frac{\Delta^2}{3}f^{(3)}_{j}+\omicron[\Delta^3]\tag{2.24}
\end{align*}

(2.25)

\begin{align*}
f^{\prime}=\frac{-3f_{j}+4f_{j+1}-f_{j+2}}{\Delta}-\frac{\Delta^2}{3}f^{(3)}_{j}+\omicron[\Delta^3]\tag{2.25}
\end{align*}

\(f^{(3)}_{j}\)は\(\frac{d^3 f}{dx^3}\)の意味ですが、数値計算ではこの項は使わないので(2.24)(2.25)はテーラー展開の2次の項までしか使っていないため、2次精度の離散化です。
 
※ここでの\(\Delta\)は、\(\Delta x\)の意味です。

※余談ですが、片側差分だと下記のような「前進差分」や「後退差分」があります。

前進差分

\begin{align*}\frac{df}{dx}=\frac{f(x+\Delta x)-f(x)}{\Delta x}\end{align*}

後退差分

\begin{align*}\frac{df}{dx}=\frac{f(x)-f(x-\Delta x)}{\Delta x}\end{align*}


 

質問内容から、p.29の(2.28),(2.29)の導出方法もわからないと思うので、合わせて解説をしたいと思います。

p.29(2.28)、(2.29)は何かというと・・・・☟これです。

ちなみに、p.29の(2.28),(2.29)は\(\frac{d^2 f}{dx^2}\)のような二階微分を、等間隔格子での片側差分でどう表されるのかという話です。
結果を示しておくと・・・

(2.28)

\begin{align*}
f^{\prime\prime}=\frac{-f_{j-3}+4f_{j-2}-5f_{j-1}+2f_{j}}{\Delta^2}+\frac{11\Delta^2}{12}f^{(4)}_{j}+\omicron[\Delta^3]\tag{2.28}
\end{align*}

(2.29)

\begin{align*}
f^{\prime\prime}=\frac{2f_{j}-5f_{j+1}+4f_{j+2}-f_{j+3}}{\Delta^2}+\frac{11\Delta^2}{12}f^{(4)}_{j}+\omicron[\Delta^3]\tag{2.29}
\end{align*}

\(f^{(4)}_{j}\)は\(\frac{d^4 f}{dx^4}\)の意味ですが、数値計算ではこの項は使わないので(2.28)(2.29)はテーラー展開の2次の項までしか使っていないため、2次精度の離散化です。
 
※ここでの\(\Delta\)は、\(\Delta x\)の意味です。

早速、導出をしていこうと思いますが、(2.24)(2.25)は基本的に同じ導出の仕方なので(2.24)のみ導出します。
また、(2.28)、(2.29)も同じ導出方法なので(2.28)のみ示します。

(2.24)(2.25)(2.28)(2.29)は全て同じ要領で導出できるのですが、それぞれ違う方法で導出してみます。

  • (2.24)(2.25)→連立方程式を代数的に解く
  • (2.28)(2.29)→連立方程式を行列計算で解く

連立方程式を行列計算で解く方法の方が便利です。

では導出をしていきましょう(^^)/

スポンサーリンク

3点片側 後退差分式(2.24)式の導出

 

カマキリ

まず、簡単なやつからやります。

 

(2.24)

\begin{align*}
f^{\prime}=\frac{f_{j-2}-4f_{j-1}+3f_{j}}{\Delta}+\frac{\Delta^2}{3}f^{(3)}_{j}+\omicron[\Delta^3]\tag{2.24}
\end{align*}

(2.24)式が意味しているのは、以下のような図の一階微分\(\frac{df}{dx}\)を求めることです。

自身の位置\(j\)より後ろ側(\(j-1,j-2\))の値を使って一階微分\(\frac{df}{dx}\)を表したいのです。

精度として保証できそうな次数までテーラー展開をしてすることで求めることができます。
今回は、3次の項までテーラー展開を行い、

出てきた以下の(1)式~(3)式を式変形することで、1階微分を差分で表現します。

\begin{align*}
f(x-\Delta x)=f(x)-\frac{df}{dx}\Delta x+\frac{1}{2}\frac{d^2f}{dx^2}\Delta x^2-\frac{1}{6}\frac{d^3f}{dx^3}\Delta x^3+\omicron(\Delta x^4)\tag{1}
\end{align*}

\begin{align*}
f(x-2\Delta x)=f(x)-2\frac{df}{dx}\Delta x+2\frac{d^2f}{dx^2}\Delta x^2-\frac{4}{3}\frac{d^3f}{dx^3}\Delta x^3+\omicron(\Delta x^4)\tag{2}
\end{align*}

以上で、2次の項までテーラー展開をした結果が用意できました。

ここまで来たら下記の連立方程式を解くことで「未知の値」を求められます。

 

 

つまり、\(\frac{d f}{dx},\frac{d^2 f}{d x^2}\)を\(f(x),f(x-\Delta x), f(x-2\Delta x)\)で表せばよいという事です。

ここまで来れば(2.24)を求めるのは簡単で、\(\frac{df}{dx}\)を求めたいのですから、
(1)×4-(2)として、\(\frac{d^2 f }{dx^2}\)を消去すれば良いです。

(1)×4-(2)

\begin{align*}
4f(x-2\Delta x)-f(x-2\Delta x)=3f(x)-2\frac{df}{dx}\Delta x+2\frac{d^2f}{dx^2}\Delta x^2+\frac{2}{3}\frac{d^3f}{dx^3}\Delta x^3+\omicron(\Delta x^4)
\end{align*}

よって、

\begin{align*}
\frac{df}{dx}=\frac{f(x-2\Delta x)-4f(x-\Delta x)+3f(x)}{\Delta x}
+\frac{\Delta x^2}{3}\frac{d^3f}{dx^3}+\omicron(\Delta x^4)\tag{3}
\end{align*}

今は、\(f(x)\)が連続であるとしてテーラー展開しましたが、実際の数値計算では格子点にしか値がないため、とびとびの値をとることになります。

👆このように考えると(3)式は、

\begin{align*}
\frac{df}{dx}=\frac{f_{j-2}-4f_{j-1}+3f_{j}}{\Delta x}+\frac{\Delta x^2}{3}\frac{d^3f}{dx^3}+\omicron(\Delta x^4)\tag{4}
\end{align*}

これが(2.24)式ということになります。

\begin{align*}
f^{\prime}=\frac{f_{j-2}-4f_{j-1}+3f_{j}}{\Delta}+\frac{\Delta^2}{3}f^{(3)}_{j}+\omicron[\Delta^3]\tag{2.24}
\end{align*}

この離散化はテーラー展開において、3次以降は無視するため2次精度の片側差分となります。

そして、自身の位置\(j\)とその後ろ側(\(j-1,j-2\))の計3点を使っているため、これは3点式の2次精度の後退差分と呼びます。

3点片側 前進差分式(2.25)の導出

 

(2.25)を導出するには、(2.24)と同じ導出過程をたどれば良いです。

今度は、自身の位置\(j\)より前側(\(j+1,j+2\))の値を使って\(\frac{df}{dx}\)を表わすことを考えます。

導出過程は、読者の課題にしておきます。

☟結果は以下のようになるので、試してみてください。

(2.25)

\begin{align*}
f^{\prime}=\frac{-3f_{j}+4f_{j+1}-f_{j+2}}{\Delta}-\frac{\Delta^2}{3}f^{(3)}_{j}+\omicron[\Delta^3]\tag{2.25}
\end{align*}

4点片側 後退差分式(2.28)式の導出

 

カマキリ

今度は、ちょっと複雑です!

(2.28)

\begin{align*}
f^{\prime\prime}=\frac{-f_{j-3}+4f_{j-2}-5f_{j-1}+2f_{j}}{\Delta^2}+\frac{11\Delta^2}{12}f^{(4)}_{j}+\omicron[\Delta^3]\tag{2.28}
\end{align*}

(2.28)式が意味しているのは、以下のような図の二階微分\(\frac{d^2f}{dx^2}\)を求めることです。

自身の位置\(j\)より後ろ側(\(j-1,j-2,j+3\))の値を使って二階微分\(\frac{d^2f}{dx^2}\)を表したいのです。

まず、4次の項までのテーラー展開を用意します。

\begin{align*}
f=f\tag{5}
\end{align*}

\begin{align*}
f(x-\Delta x)=f(x)-\frac{df}{dx}\Delta x+\frac{1}{2}\frac{d^2f}{dx^2}\Delta x^2-\frac{1}{6}\frac{d^3f}{dx^3}\Delta x^3+\frac{1}{24}\frac{d^4f}{dx^4}\Delta x^4+\omicron(\Delta x^5)\tag{6}
\end{align*}

\begin{align*}
f(x-2\Delta x)=f(x)-2\frac{df}{dx}\Delta x+2\frac{d^2f}{dx^2}\Delta x^2-\frac{4}{3}\frac{d^3f}{dx^3}\Delta x^3+\frac{2}{3}\frac{d^4f}{dx^4}\Delta x^4+\omicron(\Delta x^5)\tag{7}
\end{align*}

\begin{align*}
f(x-3\Delta x)=f(x)-3\frac{df}{dx}\Delta x+\frac{9}{2}\frac{d^2f}{dx^2}\Delta x^2-\frac{9}{2}\frac{d^3f}{dx^3}\Delta x^3+\frac{27}{8}\frac{d^4f}{dx^4}\Delta x^4+\omicron(\Delta x^5)\tag{8}
\end{align*}

以上で、4次の項までテーラー展開をした結果が用意できました。

ここまで来たら先ほどと同様に連立方程式を解けばよいという事になります。

さすがにこの連立方程式を代数的に解くのはめんどうなので、行列計算で処理したいと思います。

👆こちらを行列式にまとめて・・・・

\begin{align*}
\begin{pmatrix}
f\\
f(x-\Delta x)-\frac{\Delta x^4}{24}\frac{d^4f}{dx^4}\\
f(x-2\Delta x)-\frac{2\Delta x^4}{3}\frac{d^4f}{dx^4}\\
f(x-3\Delta x)-\frac{27\Delta x^4}{8}\frac{d^4f}{dx^4}
\end{pmatrix}
=\begin{pmatrix}
1 &0 &0 &1 \\
1 &-1 &\frac{1}{2} &-\frac{1}{6} \\
1 &-2 &2 &-\frac{4}{3} \\
1 &-3 &\frac{9}{2} & \frac{9}{2}
\end{pmatrix}
\begin{pmatrix}
f\\
\frac{df}{dx}\Delta x\\
\frac{d^2f}{dx^2}\Delta x^2\\
\frac{d^3f}{dx^3}\Delta x^3
\end{pmatrix}\tag{9}
\end{align*}

さらにまとめると・・・・行列計算にすることができます。

 

見やすさのために、右辺と左辺を逆にして・・・・

\begin{align*}
AX=Y\tag{10}
\end{align*}

行列\(X\)を求めたいので、逆行列\(A^{-1}\)を求めて(10)式の左から作用させて・・・

\begin{align*}
X=A^{-1}Y\tag{11}
\end{align*}

とすれば、行列\(X\)が求まるのがわかりますね(^^)

いったんわかりやすくするために文字で置き換えておきます。

\begin{align*}
\begin{pmatrix}
1 &0 &0 &0 \\
1 &-1 &\frac{1}{2} &-\frac{1}{6} \\
1 &-2 &2 &-\frac{4}{3} \\
1 &-3 &\frac{9}{2} &- \frac{9}{2}
\end{pmatrix}
\begin{pmatrix}
x_{1}\\
x_{2}\\
x_{3}\\
x_{4}
\end{pmatrix}
=
\begin{pmatrix}
y_{1}\\
y_{2}\\
y_{3}\\
y_{4}
\end{pmatrix}
\end{align*}

ではどうやって(11)式の\(X\)を求めるのか・・・・??

線形代数の知識を使えばいろいろな解法がありますが、一番簡単な方法として「拡大係数行列」を使って計算していく方法です。

途中式はごちゃぐちゃしていますが、「スタート」と「ゴール」だけ確認できれば良いでしょう!

(2.24)式を導出したよりもスマートにすべての変数を求めることができていますよね。

それで・・・

とにかく今は求めたい式というのが3行目の、

\begin{align*}
x_{3}=2y_{1}-5y_{2}+4y_{3}-y_{4}\tag{12}
\end{align*}

なので、

(12)式の変数に置いた式を元に戻して・・・

\begin{align*}
\frac{d^2f}{dx^2}\Delta x^2 &=2f-5\bigg(f(x-\Delta x)-\frac{\Delta x^4}{24}\frac{d^4f}{dx^4}\bigg)\\
&+4\bigg(f(x-2\Delta x)-\frac{2\Delta x^4}{3}\frac{d^4f}{dx^4}\bigg)\\
&-\bigg(f(x-3\Delta x)-\frac{27\Delta x^4}{8}\frac{d^4f}{dx^4}\bigg)\tag{13}
\end{align*}

(13)式を整理すると・・・

\begin{align*}
\frac{d^2f}{dx^2}=\frac{2f-5f(x-\Delta x)+4f(x-2\Delta x)-f(x-3\Delta x)}{\Delta x^2}-\frac{11\Delta^2}{12}\frac{d^4f}{dx^4}\tag{14}
\end{align*}

今は、\(f(x)\)が連続であるとしてテーラー展開しましたが、実際の数値計算では格子点にしか値がないため、とびとびの値をとることになります。

👆このように考えると(14)式は、

\begin{align*}
\frac{d^2f}{dx^2}=\frac{-f_{j-3}+4f_{j-2}-5f_{j-1}+2f_{j}}{\Delta^2}+\frac{11\Delta^2}{12}f^{(4)}_{j}+\omicron[\Delta^3]\tag{15}
\end{align*}

これが(2.28)式ということになります。

\begin{align*}
f^{\prime\prime}=\frac{-f_{j-3}+4f_{j-2}-5f_{j-1}+2f_{j}}{\Delta^2}+\frac{11\Delta^2}{12}f^{(4)}_{j}+\omicron[\Delta^3]\tag{2.28}
\end{align*}

この離散化はテーラー展開において、3次以降は無視するため2次精度の片側差分となります。

そして、自身の位置\(j\)とその後ろ側(\(j-1,j-2,j-3\))の計4点の値を使っているため、これは4点式の2次精度の後退差分と言えるでしょう!

念のためPythonで\(A\)の逆行列を求める

(2.28)式を手計算で求めましたが、ちょっと不安なので数値計算で本当に正しいかを見ておきます。

何を確認するかですが・・・\(A\)の逆行列を確認します。

Pythonのnumpyを使えばあっという間に逆行列が計算できます。

【結果】

分数の部分は割り算で出てきているので、わかりにくいですが逆行列\(A^{-1}\)が手計算の結果と同じですね。

どうやら計算ミスはなさそうです。

 

4点片側 前進差分式(2.29)式の導出

 

(2.29)を導出するには、(2.28)と同じ導出過程をたどれば良いです。

今度は、自身の位置\(j\)より前側(\(j+1,j+2,j+3\))の値を使って\(\frac{d^2f}{dx^3}\)を表わすことを考えます。

導出過程は、読者の課題にしておきます。

☟結果は以下のようになるので、試してみてください。

(2.29)

\begin{align*}
f^{\prime\prime}=\frac{2f_{j}-5f_{j+1}+4f_{j+2}-f_{j+3}}{\Delta^2}+\frac{11\Delta^2}{12}f^{(4)}_{j}+\omicron[\Delta^3]\tag{2.29}
\end{align*}

というわけで導出は以上です(^^)/

まとめ

 

記事の内容はこちら
「乱流の数値シミュレーション」

乱流の数値シミュレーション 改訂版

乱流の数値シミュレーション 改訂版

梶島 岳夫
4,070円(09/28 02:30時点)
発売日: 2017/01/27
Amazonの情報を掲載しています
  • p.28(2.24)、(2.25)
  • p.29(2.28)、(2.29)

を導出する。

 

\(\frac{df}{dx}\)のような一階微分を、等間隔格子での片側差分でどう表されるのかという話です。
結果を示しておくと・・・

(2.24)


\begin{align*}
f^{\prime}=\frac{f_{j-2}-4f_{j-1}+3f_{j}}{\Delta}+\frac{\Delta^2}{3}f^{(3)}_{j}+\omicron[\Delta^2]\tag{2.24}
\end{align*}

(2.25)


\begin{align*}
f^{\prime}=\frac{-3f_{j}+4f_{j+1}-f_{j+2}}{\Delta}-\frac{\Delta^2}{3}f^{(3)}_{j}+\omicron[\Delta^2]\tag{2.25}
\end{align*}

\(f^{(3)}_{j}\)は\(\frac{d^3 f}{dx^3}\)の意味ですが、数値計算ではこの項は使わないので(2.24)(2.25)はテーラー展開の2次の項までしか使っていないため、2次精度の離散化です。
 
※ここでの\(\Delta\)は、\(\Delta x\)の意味です。

p.29の(2.28),(2.29)は\(\frac{d^2 f}{dx^2}\)のような二階微分を、等間隔格子での片側差分でどう表されるのかという話です。
結果を示しておくと・・・
 

(2.28)


\begin{align*}
f^{\prime\prime}=\frac{-f_{j-3}+4f_{j-2}-5f_{j-1}+2f_{j}}{\Delta^2}+\frac{11\Delta^2}{12}f^{(4)}_{j}+\omicron[\Delta^3]\tag{2.28}
\end{align*}

(2.29)


\begin{align*}
f^{\prime\prime}=\frac{2f_{j}-5f_{j+1}+4f_{j+2}-f_{j+3}}{\Delta^2}+\frac{11\Delta^2}{12}f^{(4)}_{j}+\omicron[\Delta^3]\tag{2.29}
\end{align*}

\(f^{(4)}_{j}\)は\(\frac{d^4 f}{dx^4}\)の意味ですが、数値計算ではこの項は使わないので(2.28)(2.29)はテーラー展開の2次の項までしか使っていないため、2次精度の離散化です。
※ここでの\(\Delta\)は、\(\Delta x\)の意味です。

 

記事にして質問者さんに返答したらとても喜んでくれました(^^)/

このように、記事の書いても誰の役に立つかはわかりませんが、きっと誰かの役に立つと思って記事を書いていこうと思います_(._.)_

 

数値計算は偏微分方程式を離散化するがために、方程式の精度が悪くなったり、解が不安定になったりします。

色々と勉強することは多いです・・・

☟こちらの記事もどうぞ

【プロフィール】

カマキリ
(^^)

大学の専攻は物性理論で、Fortranを使って数値計算をしていました。
CAEを用いた流体解析は興味がありOpenFOAMを使って勉強しています。

プロフィール記事はこちら

 

大学学部レベルの物理の解説をします 大学初学者で物理にお困りの方にわかりやすく解説します。

このブログでは主に大学以上の物理を勉強して記事にわかりやすくまとめていきます。

  • ・解析力学
  • ・流体力学
  • ・熱力学
  • ・量子統計
  • ・CAE解析(流体解析)
  • note
    noteで内容は主に「プログラミング言語」の勉強の進捗を日々書いています。また、「現在勉強中の内容」「日々思ったこと」も日記代わりに書き記しています。
  • youtube
    youtubeではオープンソースの流体解析、構造解析、1DCAEの操作方法などを動画にしています。
    (音声はありません_(._.)_)
  • Qiita
    Qiitaではプログラミング言語の基本的な内容をまとめています。
関連記事もどうぞ