2020-05-16

コンパクト差分

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

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 の場合の計算結果（マーカー）と厳密解（実線）をあらわしている．計算結果は厳密解とほぼ一致しており，正しく計算できていることが確認できる．

f:id:kimmy123:20200516214418p:plain — 図1 1階微分

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

例えば，

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

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

f:id:kimmy123:20200516215130p:plain — 図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次精度片側差分

2020-05-16

陽的差分法

数値計算

テキスト「空力音の直接数値計算を支える技術」（暗黙の型宣言著）の補足課題を作ってみました．

さっそく第一回目の課題．

テキスト3.3まで読んだらやってみよう！

課題

２次，4次精度中心差分法の誤差を計算する．

ヒント：

テキストList3.1のプログラムに誤差計算のコードを追加する．
格子点数を増やした時の誤差の減り方が，テキスト図3.11に示されているように， 2次精度はslope2，4次精度はslope4にほぼ沿うことを確認する．

解答例

テキストList3.1の「!ここで誤差を計算する」に以下を追加．

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

格子点数Nを8から1024まで変化させた場合の2次精度中心差分の最小平均誤差は，表1のE2-MinAveErrのとおり．

表1. 2次，4次精度中心差分誤差

N	Δx	E2-MinAveErr	E4-MinAveErr
8	0.785398	0.060164	0.007113
32	0.19635	0.00407	3.13E-05
128	0.049087	0.000256	1.23E-07
512	0.012272	1.6E-05	4.81E-10
1024	0.006136	3.99E-06	3.01E-11

これをグラフにすると，

f:id:kimmy123:20200516214921p:plain — 図1 陽的差分法誤差

E2-MinAveErrの近似曲線を引くと， $y=0.1013x^{1.9878}$ であり，ほぼ2次精度（正確には1.9878精度）であることが確認できる．

slope2のグラフは， $y=ax^{2}$ のグラフを引けばOK．(aは任意．0.5くらいにしておくと近似曲線と重ならなくて見やすい．)

FiniteDifferenceのブロックを4次精度に書き換える．

    FiniteDifference : block
        integer :: ip1,i,im1,ip2,im2
        do i=1,N
            ip1 = modulo((i+1)-1,N)+1
            im1 = modulo((i-1)-1,N)+1
            ip2 = modulo((i+2)-1,N)+1
            im2 = modulo((i-2)-1,N)+1
            f%dx(i) = (-f%var(ip2) + 8*f%var(ip1) -8*f%var(im1) + f%var(im2))/(12.d0*dx)
        end do
    end block FiniteDifference

4次精度中心差分式はテキスト式3.18参照．

誤差を計算すると，表1のE4-MinAveErrのとおり．近似曲線は $y=0.198x^{3.9816}$ なのでほぼ4次精度．

Tips

周期境界条件ではmoduloを使うと便利

周期境界条件というと，普通はif文を使って境界での参照点を設定するけれど，本テキストではmoduloをうまく使ってif文をいちいち書かなくて済むようにしている．

modulo((i+s)-1,N)のsは，いくつ先の格子なのかというのをあらわしている．プログラムを見やすくするためにsを明記している．

ちなみにmodulo((i+s)-1,N)+1なんてしないで，このほうが簡単でしょ？と思いつつ，

ip1 = modulo(i+1,N)
im1 = modulo(i-1,N)

と書き換えて，ip1とim1を出力させてみると・・・ ※簡単のためN = 8

N	ip1	im1
0	1	2
1	2	3
2	3	4
3	4	5
4	5	6
5	6	7
6	7	0
7	8	1

im1=0が配列外参照になってしまう．（fortran配列はデフォルトで1はじまり）

これを防ぐために，modulo内で-1して，あとから+1している．

誤差は格子点数Nで平均化する！

前述のCmpErrorでerr/Nとせず，単にerrだけにしてしまうと，格子全体での累積誤差になる．

格子点数を増やすほど誤差も多く足されるので，誤差は減らない．これはテキストの「ノルム」のl1ノルムで述べられているとおり．

Kimmy's Lab.

フルート好きのリケジョkimmyによる裏研究所

コンパクト差分

課題

解答例

陽的差分法

課題

解答例

Tips