コンパクト差分

テキスト「空力音の直接数値計算を支える技術」（暗黙の型宣言著）の補足課題２

4章まで読んだらやってみよう！

課題

(1) コンパクト差分を用いて正弦波の1階微分を計算するプログラムを作成する．
※誤差も計算して，実効的な精度を確認する．
※連立方程式はLU分解(Crout法)を用いる．
※テキストの以下プログラム参照．

• List3.2（コンパクト差分main）
• List3.4（1階微分に対するコンパクト差分用module）
• List4.2（3重対角行列定義）
• List4.3（LU分解）

(2) 同様に，2階微分を計算するプログラムを作成する．
※テキストのList3.11(2階微分に対するコンパクト差分用module)参照．

解答例

まずは1階微分．テキストを写経して，以下３つのファイルを作成する．

• main.f90(List3.2参照)
• module_CompactScheme.f90(List3.4参照)
• TypeTriDiagonalMatrix.f90(List4.2参照)

コードに以下を追加する．

• List4.2(module TypeTriDiagonalMatrix)の「ここに解法のサブルーチンを書く」に，List4.3(subroutine decompose, solve)を追加．

• List3.2(program main)の「ここで誤差を計算する」に以下を追加．
  ※厳密解cos(x)に対する最小平均誤差を計算している．

    CmpError : block
        integer :: i
        real(8) :: x, err
        err = 0d0
        do i=1,Nx
            x = dble(i-1)*dx
            err = err + abs(f%dx(i)-cos(x))
            !print * , f%dx(i)
        end do
        print * , "MinAveError:", err/Nx
    end block CmpError   

以上で1階微分計算用プログラムは完成．
動作確認のため，上記コードでコメントアウトしたprint分をコメントインしてみて，正しく1階微分が計算できているか確認してみる．

図1は，Nx=8 の場合の計算結果（マーカー）と厳密解（実線）をあらわしている．計算結果は厳密解とほぼ一致しており，正しく計算できていることが確認できる．

図1 1階微分

格子分割数(Nx)やスキーム(Interior1st, Boundary1st, NearBdry1st)を変えて，誤差の変化を確認してみる．

例えば，

• グローバルスキーム・・・6次精度
• 境界1点内側・・・5次精度
• 境界・・・4次精度片側差分

を用いた場合の，格子間隔Δxと誤差は，図2のC6N5B4T1のとおり．

図1 1階微分誤差

これに近似直線を引くと，実質的な誤差は4.9178なので，約5次精度であることがわかる．

境界1点内側を4次精度(C6N4B4T1)にしてみても，精度はほとんど変わらない．

グローバルスキームは6次精度だが，境界が4次もしくは5次精度なので，全体としては約5次精度に落ちる．

次に2階微分．プログラムの修正箇所は以下．

• module_CompactScheme.f90にList3.11のソースを追加．
• main.f90で呼び出すモジュールを2階微分用のモジュールに変更．

        call construct2ndDerivativeMatrix(G)
        call construct2ndDerivativeVector(f%dx, f%var, dx, G)

• main.f90の誤差計算ブロック内の厳密解を-sin(x)に変更．

err = err + abs(f%dx(i)-(-sin(x)))

以下のスキームを用いた場合の誤差は図3のとおり．実質，約4次精度．

• グローバルスキーム・・・6次精度
• 境界1点内側・・・4次精度
• 境界・・・3次精度片側差分

図3 2階微分誤差