陽的差分法

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

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

テキスト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

これをグラフにすると，

図1 陽的差分法誤差

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

slope2のグラフは，のグラフを引けば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のとおり．近似曲線はなのでほぼ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ノルムで述べられているとおり．