.. sectionauthor:: 清水康行 `本テキストへの感想, 質問, 要望, バグ報告などはこちらへお願いします. `_ ======================== 高精度風上差分法 ======================== 変数の格子間分布の補完 ========================== いま， :numref:`2dcip_1` に示すように黒線で示す三角形分布の変数 :math:`f` の時刻 :math:`t` における分布 :math:`f(t)` が，速度 :math:`u(>0)` で :math:`x` 方向に移流する状態を想定します．なお，時刻 :math:`t` においては 三角形の頂点である :math:`A` 点は :math:`i` 番目の格子点に位置しているとします． :math:`\Delta t` 時間で移動する距離は :math:`u \Delta t` なので，図の赤の破線の分布に移動するはずで， ピークは赤丸の :math:`A'` 点に移動することになります． .. _2dcip_1: .. figure:: images/01/2dcip_1.png :width: 70% ::変数の移流(1) ところが，:numref:`2dcip_1` のように実際は :math:`A'` 点は計算格子上には無く, :math:`i` 番目の格子と :math:`i+1` 番目の 格子の間に来てしまうので，移動後の :math:`A'` 点を頂点とする三角形分布を離散値で表現することは出来ません． これを，現存の格子点上で表現しようとすると, :numref:`2dcip_2` の青線のような図形になってしまって， 結果的にピークの :math:`A'` 点は :math:`B_1` 点と :math:`B_2` 点のようにピークが減衰した形でしか表現出来ません． これは, 前章の「 :ref:`art_dif` 」で述べた数値拡散を図形的に説明したものとも言えます． .. _2dcip_2: .. figure:: images/01/2dcip_2.png :width: 70% ::変数の移流(2) 以上の説明は，時刻が :math:`t` から :math:`t+\Delta t` に進む 時の説明でしたが，これを，時刻が :math:`t-\Delta t` から :math:`t` に進むとし場合に置き換えて見ましょう． :numref:`2dcip_3` に示すように, 時刻 :math:`t` における :math:`f` の分布が正しく :math:`A` を 頂点とえうる三角形分布で表現される ためには，時刻 :math:`t-\Delta t` における :math:`A'` を頂点とする三角形分布が正しく表現されている必要があります． しかしながら, :math:`A'` 点は格子点上に無いため，離散値でこの形を表現することは不可能です． そこで，時刻 :math:`t-\Delta t` における分布形は離散値ではなく，何等かの連続関数である必要があります． これが，次に説明するCIP法の起源となります． .. _2dcip_3: .. figure:: images/01/2dcip_3.png :width: 70% ::変数の移流(3) .. _label1: `CIP法 `_ の基本概念 ================================================================== 1次元移流方程式を再記します． .. math:: :label: kiso_11 \cfrac{\partial f}{\partial t}+u\cfrac{\partial f}{\partial x}=0 なお，前章の :eq:`kiso_1` では移流速度を :math:`c` としていましたが， ここでは :math:`u` としてあります．(特に意味も理由も無いのですが，流体の計算では移流速度は流体の流速 が用いられ，流速は :math:`u, v, w` などが用いられる場合が多いのでここでは :math:`u` としました．) :eq:`kiso_11` 式の1次元移流方程式を数値的に計算するということは, 要するに :math:`f(x,t+\Delta t)` をどう推定するかということです. 前節で説明したように，通常の風上差分ではこれを :math:`f(x-\Delta x,t)` と :math:`f(x,t)` を使って計算します．ただし，前記のようにこの方法は計算は安定ですが， 計算結果が拡散してしまって（数値拡散）正確ではありません． そこで，ここでは一次精度風上差分に代わる高精度の差分法である矢部ら [Ref:1]_ によって開発された `CIP法 `_ を紹介します． そもそも :eq:`kiso_11` 式の意味するところは 変数 :math:`f` の分布形を速度 :math:`u` で伝えているのですから, :math:`\Delta t` 時間の間には :math:`f` の分布形が伝わる距離は :math:`u \Delta t` のはずです. 即ち, 将来の :math:`f=f(x, t+\Delta t)` は 現在の :math:`f=f(x-u\Delta t, t)` がそのまま伝わっているはずです( :numref:`trans` 参照). .. _trans: .. figure:: images/01/fig1-2.png :width: 50% ::math:`f(x,t+\Delta t)` の伝わり方 将来の :math:`f(x,t+\Delta t)` を推定する代わりに, 現在の :math:`f(x-u\Delta t, t)` が推定できれば, その値をそのまま用いることが可能となるわけです. これを式で書けば, .. math:: :label: yosoku f(x,t+\Delta t) \simeq f(x-u\Delta t, t) となります. 格子間の変数分布形推定方法 ============================= 3次関数の導入 ------------------- :math:`f(x-u\Delta t, t)` を推定する方法はいろいろあると思いますが, 問題は :math:`f` は離散量として与えられているので, :math:`f` の既知量は :numref:`trans` の●点であり, この2つの●の点から○の :math:`f(x-u\Delta t, t)` を推定する必要があるということになります. 最も単純には2点間を直線で結んで比例配分するという方法でしょうが, 2点の間で :math:`f` が急激に変化している場合や, 1度増加して減少する場合などを 考慮できるためには3次関数で表すのが妥当です. .. _moshiki: .. figure:: images/01/fig1-3.png :width: 50% : CIP法の説明模式図 :numref:`moshiki` に示すように :math:`x` 方向に :math:`i` , :math:`t` 方向に :math:`n` で離散量を表し, A点およびB点における :math:`f` を :math:`f(i-1,n)` および :math:`f(i,n)` とします. AB間における :math:`f` の分布を3次関数で, :math:`X` をBからAに向かった距離として表すと .. math:: :label: fvalue F=a_1 X^3 + a_2 X^2 + a_3 X + a_4 となり, :math:`a_1 \sim a_4` の4個のパラメータを含む3次式となります. ここで注意して頂きたいのは，:math:`f` は各格子点で定義される離散量であるのに 対して， :math:`F` は各格子点の :math:`f` を格子の間を補完する連続関数 (3次式)であるという点です． :eq:`fvalue` 式中の4個の パラメータを決定するためには4つの条件が必要となります. このうち2つは, .. math:: :label: bc0 \mbox{B点で:} & F_{(X=0)}=f(i,n) \\ \mbox{A点で:} & F_{(X=\Delta x)} = f(i-1,n) の条件で決まりますが, あと2つはどうしたら良いでしょうか? これもいろいろ考えられます(たとえば, :math:`f(i-2,n)` や :math:`f(i+1,n)` などを用いる) がCIP法では, :math:`f` の :math:`x` 軸方向の勾配を用います. つまり, .. math:: :label: bc1 \mbox{B点で:} & \displaystyle{{\partial F}\over{\partial x}}_{(X=0)}= {{\partial f}\over{\partial x}} (i,n) \\ \mbox{A点で:} & \displaystyle{{\partial F}\over{\partial x}}_{(X=\Delta x)} = {{\partial f}\over{\partial x}}(i-1,n) :eq:`bc0` :math:`\sim` :eq:`bc1` 式によって :math:`i` と :math:`i-1` の情報のみでAB間の :math:`f` の分布 :math:`F` が推定できることになり, これもCIP法の大きなメリットの一つです. :eq:`fvalue` 式およびこれを :math:`x` で偏微分した, .. math:: :label: fprof F_x= {{\partial F}\over{\partial x}} =3 a_1 X^2 + 2 a_2 X +a_3 に :eq:`bc0` :math:`\sim` :eq:`bc1` 式の4つの条件を適用して, :math:`a_1 \sim a_4` のパラメータを決定して :math:`F` を求めれば良いのですが, これはちょっと面倒そうなので :eq:`bc0` :math:`\sim` :eq:`bc1` 式の内, B点の条件( :math:`X=0` )を自動的に満たす .. math:: :label: fxnew F(X)=\left[ (aX + b)X + f_x(i,n) \right] X + f(i,n) をあらためて用いることとします. ただし, .. math:: :label: fnew1 f_x(i,n) = \displaystyle{{\partial f}\over{\partial x}} (i,n) です. :eq:`fxnew` 式を書き直して, .. math:: :label: fnew2 F(X)=aX^3+bX^2+f_x(i,n)X+f(i,n) となり, これを :math:`x` で偏微分すると, .. math:: :label: fnew3 F_x(X)=3aX^2+2bX+f_x(i,n) となります． :eq:`fnew2` 式と :eq:`fnew3` 式は :eq:`bc0` 式と :eq:`bc1` 式のうちB点の条件( :math:`X=0` )を自動的に満たしています． あとは，:eq:`bc0` 式と :eq:`bc1` 式のうちB点の条件( :math:`X=\Delta x` ) を :eq:`fnew2` 式と :eq:`fnew3` 式に適用して, .. math:: :label: ab1 a\Delta x^3 + b \Delta x^2 + f_x(i,n) \Delta x + f(i,n) = f(i-1,n) \\ 3a\Delta x^2 + 2b\Delta x + f_x(i,n) = f_x(i-1,n) これより, .. math:: :label: ab2 a=\displaystyle{{ \left[(f_x(i-1,n)+f_x(i,n) \right] \Delta x + 2 \left[f(i,n)-f(i-1,n) \right] }\over{\Delta x^3}} \\ b=\displaystyle{{ 3\left[(f(i-1,n)-f(i,n) \right]- \left[f_x(i-1,n)+2f_x(i,n) \right] \Delta x }\over{\Delta x^2}} が求まります. 演習問題(1) ^^^^^^^^^^^^^^ .. list-table:: * - :eq:`ab2` 式の誘導をして下さい． CIP法による格子間変数分布式の推定法のまとめ -------------------------------------------------- 以上をまとめると, :eq:`kiso_11` 式にCIP法を適用して :math:`f(x,t+\Delta t)` を 求める手順は以下となります． 1) :eq:`ab2` 式で :math:`a, b` を求める. 2) :eq:`fxnew` 式に :math:`a, b` を代入して, :math:`X=u\Delta t` とすることで, :math:`F(X=u\Delta t)=f(x-u\Delta t, t)` を求める． 3) 同時に :math:`F_x(X=u\Delta t)=f_x(x-u\Delta t, t)` は :eq:`fnew3` 式で求める． 4) :eq:`yosoku` 式で :math:`f(x,t+\Delta t)` を求める． これを :math:`i=1 \sim N_i` で実行し, 終了後に時間更新する. Pythonコードと計算例 ------------------------ CIP法による1次元移流方程式の数値解法のPythonコードの例を https://github.com/YasuShimizu/Python-Code/blob/main/cip.py に示します． この計算を実行すると，以下のような結果が得られます． ただし, :math:`x_b` = 10, :math:`x_p` = 10, :math:`f_p` = 0.5, :math:`u` = 0.5 としてあります． .. _cip_tri: .. figure:: images/01/cip.gif :width: 90% : 1次元移流方程式(CIP法) このようにCIP法を用いると初期条件の三角形の分布系が崩れることなく，一定速度で 下流側に進んでいく様子が計算されていることが分かります． Appendixに掲載している理論解( :numref:`adv_anim` )と比較してもほぼ同じ結果 が得られています． 移流方向が変化する移流方程式 ================================== 前記の :ref:`label1` では移流方向が一定である場合の計算法について述べましたが，ここでは移流方向が変化する 場合にいて述べます． 1次元の移流方程式 :eq:`kiso_11` 式において，これまでの説明では移流速度 :math:`u` は一定 として扱って来ましたが， ここでは :math:`u` は一定ではなく，時間・空間的に変化する可変量とします．一般に流体では :math:`u` は流れの向きに よって正だったり負だったりしますので， :math:`u>0` でも :math:`u<0` でも対応可能なアルゴリズムを用意 する必要があります． *u<0* の場合 ---------------- .. _moshiki_u-: .. figure:: images/01/fig1-4.png :width: 50% : :math:`u<0` の場合の模式図 上図のように :math:`u<0` の場合 .. math:: :label: X_nega X=-u\Delta t とおいて，:math:`i` と :math:`i+1` の間の :math:`f` の分布を次式で表します. ただし :math:`X>0` です． .. math:: :label: X_nega_1 F(X)=\left[ (aX + b)X + f_x(i) \right] X + f(i) また, :math:`F(X)` の微分量は, .. math:: :label: X_nega_2 F'(X)=3aX^2+2bX+f_x(i) ここで, :math:`f_x` は .. math:: :label: X_nega_3 f_x = {{\partial f}\over{\partial x}} のことです. :math:`F` および :math:`F'` は .. math:: :label: X_nega_4 X=0 \mbox{で:} & F(X)=f(i) \mbox{および} F'(X)=f_x(i) \\ X=\Delta x \mbox{で:} & F(X)=f(i+1) \mbox{および} F'(X)=f'(i+1) の条件より, .. math:: :label: X_nega_5 a=\displaystyle{{ \left[(f_x(i+1)+f_x(i) \right] \Delta x + 2 \left[f(i)-f(i+1) \right] }\over{\Delta x^3}} \\ b=\displaystyle{{ 3\left[(f(i+1)-f(i) \right]- \left[f_x(i+1)+2f_x(i) \right] \Delta x }\over{\Delta x^2}} が得られます. 演習問題(2) ^^^^^^^^^^^^^^^^ .. list-table:: * - :eq:`X_nega_5` 式の誘導をして下さい． *u>0* の場合 ----------------- .. figure:: images/01/fig1-5.png :width: 50% : :math:`u>0` の場合の模式図 上図のように :math:`u>0` の場合も .. math:: :label: X_posi_1 X = - u \Delta t と起きます. 前節では :math:`X = u \Delta t (>0)` としましたが, ここでは :eq:`X_posi_1` 式のように :math:`X<0` となるように設定します. この理由はあとで分かりますが, :math:`u<0` でも :math:`u>0` でも同じアルゴリズムを使用可能とするためのものです. :math:`u>0` の場合には, :math:`i` と :math:`i-1` の間の :math:`f` の分布を次式で表します. これは :eq:`X_nega_1` 式と全く同じ形です. .. math:: :label: X_posi_2 F(X)=\left[ (aX + b)X + f_x(i) \right] X + f(i) また, :math:`F(X)` の微分量は, .. math:: :label: X_posi_3 F'(X)=3aX^2+2bX+f_x(i) で表されます. :math:`F` および :math:`F'` は .. math:: :label: X_posi_4 X=0 \mbox{で: } & F(X)=f(i) \mbox{および} F'(X)=f_x(i) \\ X= -\Delta x \mbox{で: } & F(X)=f(i-1) \mbox{および} F'(X)=f'(i-1) の条件より, .. math:: :label: X_posi_5 a=\displaystyle{{ \left[(f_x(i-1)+f_x(i) \right] \Delta x - 2 \left[f(i)-f(i-1) \right] }\over{\Delta x^3}} \\ b=\displaystyle{{ 3\left[(f(i-1)-f(i) \right]+ \left[f_x(i-1)+2f_x(i) \right] \Delta x }\over{\Delta x^2}} が得られます. 演習問題(3) ^^^^^^^^^^^^^^^ .. list-table:: * - :eq:`X_posi_5` 式の誘導をして下さい． *u>0* でも *u<0* でも使えるようにするには? -------------------------------------------- :eq:`X_nega_5` 式が :math:`u<0` の場合で, :eq:`X_posi_5` 式が :math:`u>0` の場合ですが, これらのペアを比べてみると, 大変良く似ています. そこで， .. math:: :label: sig \mbox {sign}(Y) & =1 & ,Y>0 \\ \mbox {sign}(Y) & =-1 & ,Y<0 で定義されるsign関数(シグナル関数もしくは符号関数とも呼ばれます)を利用して， .. math:: :label: sig_1 i_m = i - \mbox{sign}(u) という :math:`i_m` を定義すると, .. math:: :label: sig_2 u>0 \mbox{のとき} i_m & = i-1 \\ u<0 \mbox{のとき} i_m & = i+1 となります. したがって, :eq:`X_nega_5` 式および :eq:`X_posi_5` 式の代わりに, .. math:: :label: ab a=\displaystyle{{ \left[(f_x(i_m)+f_x(i) \right] \Delta x \cdot \mbox{sign}(u)- 2 \left[f(i)-f(i_m) \right] }\over{\Delta x^3 \cdot \mbox{sign}(u)}}\\ b=\displaystyle{{ 3\left[(f(i_m)-f(i) \right]+ \left[f_x(i_m)+2f_x(i) \right] \Delta x \cdot \mbox{sign}(u)}\over{\Delta x^2}} とすることで, :math:`u<0` の場合も :math:`u>0` の場合も同じ式を用いることが出来ます． 例題 ^^^^^^^ .. list-table:: * - 移流速度 :math:`u` が次式で表されるような正負に単振動する場合の1次元移流方程式をCIP法で計算してみましょう． * - 初期条件は :numref:`triangle` と同じとします． .. math:: :label: uprof u=u_0 \sin \left( \cfrac{2\pi}{T_l} t \right) ただし, :math:`u_0` は振動流速の最大値，:math:`T_l` は振動流速の周期である． Pythonコードと計算例 ------------------------- 移流方向が正負に振動する場合の移流方程式の CIP法に数値解法のPythonコードの例を https://github.com/YasuShimizu/Python-Code/blob/main/cip_pm.py に示します． このコードで，:math:`u_0` =2.0, :math:`T_l` =100 :math:`x_b` = 10, :math:`x_p` = 20, :math:`f_p` = 0.5とした場合の計算結果を :numref:`pm` に示します． .. _pm: .. figure:: images/01/cip_pm.gif :width: 50% : 移流速度が正負振動する場合の計算結果 Source項を含む移流方程式 ================================= ここでは, 流体の数値計算を行う場合に不可欠となる分離解法の説明を行います. 流体の方程式の多くは, :eq:`kiso_11` 式のように右辺がゼロの場合は稀で, 多くは 右辺に値を持ちます. これを一般的に次式で表します. .. math:: :label: sr1 \cfrac{\partial f}{\partial t}+u\cfrac{\partial f}{\partial x}=G :math:`G` は流れの運動方程式の場合, 水面勾配, 河床勾配, 摩擦など様々 外力項(Source項)となりますが, 要するに移流項以外の成分のことです. :eq:`sr1` 式を便宜的に以下の2つの式に分けます. .. math:: :label: bunri1 {{\partial f}\over{\partial t}} = G .. math:: :label: bunri2 {{\partial f}\over{\partial t}}+u{{\partial f}\over{\partial x}} = 0 :math:`f(t)` および :math:`f(t+\Delta t)` をそれぞれ :math:`f^n` および :math:`f^{n+1}` とし, :eq:`sr1` 式で :math:`f^n \rightarrow f^{n+1}` を計算するのを, 一旦 :eq:`bunri1` 式で :math:`f^n \rightarrow \widetilde{f}` の計算を行い, さらに :eq:`bunri2` 式で :math:`\widetilde{f} \rightarrow f^{n+1}` を計算するというように 2段階で計算します. これを分離解法と言います. このように分離 することによって初めて :eq:`bunri2` 式の形でCIP法が使えるようになります. この2段階はNon-Advection Phaseおよび, Advection Phaseと呼ばれます. CIP法では, :math:`f` の時間変化の他に :eq:`ab` 式に表れる :math:`f_x` 即ち空間微分量も 常に必要となるため, 微分量も同時に, Non-Advecion Phaseで :math:`\displaystyle{{{\partial f}\over{\partial x}}^n \rightarrow \widetilde{{\partial f}\over{\partial x}}}`, Advecion Phaseで :math:`\displaystyle{\widetilde{{\partial f}\over{\partial x}} \rightarrow {{\partial f}\over{\partial x}}^{n+1}}` というように分離して計算します． 以上を整理すると以下のようになります. .. _bunri_hyo: .. list-table:: 分離解法 :header-rows: 1 * - Phase名 - 使用する数式 - :math:`f` の更新 - :math:`\cfrac{\partial f}{\partial x}` の更新 * - **[A]** \ Non-advection Phase - :math:`\displaystyle{{\partial f}\over{\partial t}}=G` - **(a)** \ :math:`f^n \rightarrow \widetilde{f}` - **(b)** \ :math:`\displaystyle{{{\partial f}\over{\partial x}}^n \rightarrow \widetilde{{\partial f}\over{\partial x}}}` * - **[B]** \ Advection Phase - :math:`\displaystyle{{{\partial f}\over{\partial t}}+u{{\partial f}\over{\partial x}}}=0` - **(c)** \ :math:`\widetilde{f} \rightarrow f^{n+1}` - **(d)** \ :math:`\displaystyle{\widetilde{{\partial f}\over{\partial x}} \rightarrow {{\partial f}\over{\partial x}}^{n+1}}` 以上のように分離した方程式を用いてそれぞれの段階での計算方法を以下に 示します. なお以下の説明においては移流項以外を差分表示する必要が ある場合は, 出来るだけ精度の高い中央差分を用いることととします. Non-Advection Phase ------------------------------- :eq:`bunri1` 式より .. math:: :label: na1 {{\widetilde{f}-f^n}\over{\Delta t}} = G^n .. math:: :label: na2 \widetilde{f} = f^n +G^n \Delta t 計算点を :math:`i` とすると, .. math:: :label: na3 \widetilde{f}(i) = f^n(i) +G^n(i) \Delta t でNon-Advection Phaseの :math:`f` が更新されます( :math:`f^n \rightarrow \widetilde{f}` ). :eq:`na2` の両辺を :math:`x` で微分して, .. math:: :label: na4 \widetilde{{\partial f}\over{\partial x}} ={{\partial f}\over{\partial x}}^n +{{\partial G}\over{\partial x}}^n \Delta t 上式の :math:`\displaystyle{{{\partial G}\over{\partial x}}^n}` を中央差分で表すと, .. math:: :label: na5 \widetilde{{\partial f}\over{\partial x}} ={{\partial f}\over{\partial x}}^n +{{G^n(i+1)-G^n(i-1)}\over{2\Delta x}} \Delta t ここで, :eq:`na3` 式より, .. math:: :label: na6 G^n(i+1)={{\widetilde{f}(i+1)-f^n(i+1)}\over{\Delta t}} \\ G(i-1)={{\widetilde{f}(i-1)-f^n(i-1)}\over{\Delta t}} なので, :eq:`na6` を :eq:`na5` に代入して, .. math:: :label: na8 \widetilde{{\partial f}\over{\partial x}}(i)={{\partial f}\over{\partial x}}^n(i)+{{1}\over{2\Delta x}} \left[\widetilde{f}(i+1)-f^n(i+1)-\widetilde{f}(i-1)+f^n(i-1)\right] となります. したがって, Non-Advection Phaseにおける変数 :math:`f` および その空間微分量 :math:`\displaystyle{{{\partial f}\over{\partial x}}}` の更新は, 表 :numref:`bunri_hyo` 中の(a)の部分は :eq:`na3` 式で,表 :numref:`bunri_hyo` 中の(b)の部分は :eq:`na8` 式で求めることが出来ます. Advection Phase --------------------------- 変数 :math:`f`の更新は表 :numref:`bunri_hyo` 中の(c)の部分で行いますが， これは今まで解説してきたCIP法をそのまま適用できます. 即ち :eq:`X_nega_1` 式または :eq:`X_posi_2` 式(同じものですが)を適用して, .. math:: :label: ir1 {f(i)}^{n+1}=F(X)= \left[(aX+b)X+ \widetilde{{\partial f}\over{\partial x}}(i) \right] X +\widetilde{f}(i) ただし, :math:`a` および :math:`b`は :eq:`ab` 式で与えられます. 次にAdvection Phaseにおける空間微分量 :math:`\displaystyle{{\partial f}\over{\partial x}}` の更新法を考えて見ます[表 :numref:`bunri_hyo` 中の(d)の部分]. :eq:`bunri2` 式の両辺を :math:`x` で偏微分します. .. math:: :label: ad11 {{\partial \over \partial x}}\left( {{\partial f}\over{\partial t}} \right) + {{\partial \over \partial x}} \left( u {{\partial f}\over{\partial x}} \right) =0 ここで :math:`\displaystyle{{\partial f}\over{\partial x}}` を :math:`f'` と 置いて第2項目を展開すると, .. math:: :label: ad12 {{\partial f'}\over{\partial t}} + u{{\partial f'}\over{\partial x}} = -{{\partial u}\over{\partial x}} f' となります.いまここで考えているのは, Advection Phaseでの :math:`\displaystyle{{\partial f}\over{\partial x}}` の更新ですので, この式を用いて, :math:`\widetilde{f'} \rightarrow {f'}^{n+1} \left( \displaystyle{ \widetilde{{\partial f}\over{\partial x}} \rightarrow {{\partial f}\over{\partial x}}^{n+1}} \right)` を行えば良いことになります. ここで，:eq:`ad11` 式と :eq:`sr1` 式を比べてみましょう. :eq:`ad12` 式は :eq:`sr1` 式で :math:`f` を :math:`f'` に, :math:`G` を :math:`\displaystyle{-{{\partial u}\over{\partial x}} f'}` としたものに一致しています. そこで :eq:`sr1` 式でやったのと同じように :eq:`ad11` 式に対して分離解法を適用します. 即ち, :eq:`ad11` を .. math:: :label: br_a1 {{\partial f'}\over{\partial t}} + u{{\partial f'}\over{\partial x}} = 0 .. math:: :label: br_a2 {{\partial f'}\over{\partial t}} = -{{\partial u}\over{\partial x}} f' のように分離して, :math:`\widetilde{f'}` と :math:`{f'}^{n+1}` の間に中間量として :math:`\widehat{f'}` を考え, :eq:`br_a1` 式で :math:`\widetilde{f'} \rightarrow \widehat{f'}` を, :eq:`br_a2` 式で :math:`\widehat{f'} \rightarrow {f'}^{n+1}` を計算することにします. まず :eq:`br_a1` 式ですが, これに再度CIP法を適用するのも一つの方法ですが, CIP法で用いられる :math:`F(X)` の空間微分量である :eq:`X_nega_2` 式または :eq:`X_posi_3` 式(どちらも同じ形です)の :math:`F'(X)` を用いる方が賢明のようです. 即ち, .. math:: :label: br_a3 \widehat{f'}(i) = \widetilde{F'}(X) として, .. math:: :label: br_a4 \widehat{f'}(i) = 3aX^2 + 2bX + \widetilde{f'}(i) で計算します. :eq:`br_a2` 式のほうは, .. math:: :label: br_a5 {{{f'}^{n+1}(i) - \widehat{f'}(i)}\over{\Delta t}} = - {{u(i+1)-u(i-1))}\over{2\Delta x}} \widehat{f'}(i) と中央差分を適用し, これより .. math:: :label: br_a6 {f'}^{n+1}(i) = \widehat{f'}(i) -{{u(i+1)-u(i-1))}\over{2\Delta x}} \widehat{f'}(i) \Delta t で計算します. 演習問題(4) ^^^^^^^^^^^^^^ .. list-table:: * - :eq:`sr1` 式の右辺 :math:`G` が :math:`\displaystyle{D{{\partial^2 f}\over{\partial x^2}}}` のとき, .. math:: \displaystyle{{{\partial f}\over{\partial t}}+u{{\partial f}\over{\partial x}}= D{{\partial^2 f}\over{\partial x^2}}} .. list-table:: * - のような移流 :math:`\cdot` 拡散方程式となります. これを分離解法で計算する方法について述べよ. Pythonコードと計算例 ------------------------ 移流拡散方程式の計算用Pythonコードの例を https://github.com/YasuShimizu/Python-Code/blob/main/cip_pm_diffusion.py に示します． このコードで拡散係数 :math:`D=0.5` , :math:`u_0` =2.0, :math:`T_l` =100 と した場合の計算結果を :numref:`pm_dif` に示します．左右に移動しながら最初は三角形の分布形状が 徐々に拡散していく様子が表わされています． .. _pm_dif: .. figure:: images/01/cip_pm_dif.gif :width: 50% : 移流速度が正負振動する場合の移流拡散方程式計算結果 なお，移流項の高精度差分法はCIP法以外にも数多く提案されており，そのうちいくつかを，Appendix ( :ref:`others`) で紹介しておきます．