フルートを科学する ～ 室温とピッチ ～
熱くてムシムシする～
こんなときにフルート吹いてると，なんだか音が鳴りにくい(><)
頭もボーっと・・・熱中症！？
防音室の空気が入れ替わって，快適な温度と湿度にしてから吹いてみると，いつもどおり鳴る．
温度や湿度によって吹きやすさが変わるのは，なぜだろう？
経験上，ピッチはすごくズレる．あまりにズレすぎて鳴りにくいのかな．
温度や湿度によってピッチがズレるのは，音速が変わるからだろう．
湿度はちょっと話がややこしいので，温度によってどれくらいピッチが変わるのか計算してみた．
温度と音速
音速は以下の式であらわされる．
空気の場合， kg/mol．
温度がのとき，絶対温度は，
℃．
式(1)から，温度と音速は以下のようになる．
表1 温度と音速
|
温度 |
15 | 19 | 20 | 21 | 25 |
|---|---|---|---|---|---|
|
音速 |
340.7 | 343.0 | 343.6 | 344.2 | 346.5 |
これがどのくらいピッチに影響するのか？
温度とピッチ
周波数と音速の関係は以下の式であらわされる．
ここでは波長．これは，ほぼ管の長さ
で決まるはず．（開口端補正や音孔の影響は無視）
フルートは両側開端なので，管の長さは，ほぼ半波長分の長さ．(
)
式(2)は，
20℃の環境で，いつものラの音，つまりA5(884Hz)の音でチューニングしたとすると，
低いラの音(A4: 442Hz)の管長は388 mm程度．（A4の倍音がきっちりA5と仮定）
音速が表1のように変化すると，A4の周波数は以下のように変化する．
表2 温度と周波数
|
温度 |
15 | 19 | 20 | 21 | 25 |
|---|---|---|---|---|---|
|
A4周波数 |
438.2 | 441.2 | 442.0 | 442.8 | 445.8 |
|
A5周波数 |
876.4 | 882.5 | 884.0 | 885.5 | 891.5 |
|
A6周波数 |
1314.6 | 1323.7 | 1326.0 | 1328.3 | 1337.3 |
温度が1℃上下すると，周波数は0.75Hzくらい上下する．
このくらいなら許容範囲？でも5℃上下すると3.7Hzくらい上下する．これは聴いてわかる？
同様にして，中音と高音のラについても周波数を求めてみると，高音ほど温度による周波数変化は大きい．
人間の耳がこれをどう感じるかはわからない．
もしかしたら高いラのピッチのズレは，低いラのズレより気にならないかもしれない．
人間の耳の感度はさておき．
このピッチの上がり下がりを補正するために，管をどれくらい抜き差ししないといけないか見積もってみた．
ふつうはA5でチューニングすると思うので，管をだけ抜き差しして，
にすることになる．これを満たすを計算すると以下のようになる．
表3 管の抜き差し必要量
|
温度 |
15 | 19 | 20 | 21 | 25 |
|---|---|---|---|---|---|
|
管長調整量 |
-1.7 | -0.3 | 0.0 | 0.3 | 1.7 |
温度が1℃上がる（下がる）と，0.3mmくらい抜く（差し込む）必要がある．
5℃変化すると1.7mmくらい．
今回の計算は，開口端補正などを無視してシンプルに計算しているので，実際はこのとおりではありません！
私の場合，夏と冬で練習環境の温度差は10℃くらいあり，3mmくらい管を抜き差ししているので，だいたいこの計算どおりか．
日本の夏場は，湿度も上がるので，音速がやや速くなり，さらに管を抜く必要が出てくると思われる．それが何ミリかは・・・興味ある人は見積もって！
フルートを科学する ～楽器の支え～
最近、フルートを吹いていて、左手人差し指の付け根が痛くなることがあり、 気になって考えてみた。
フルート演奏時の力・モーメントのつり合い
フルートは，主に，唇，左手人差し指付け根、右手親指で支える（３点支持）．
それにプラスして，補助的に，右手小指も楽器を支える役割をしている．
唇は「支点」としての役割をしていると思われる．
模式的にあらわすと以下の図のようになる．
ここで，Rは唇による支点反力をあらわす．
楽器を支えるためには，力とモーメントがつり合っていないといけない。
- 力のつり合い
- 唇まわりのモーメントのつり合い
a:b:c=3:2:1くらいとすると，
式②から，
つまり，右小指の力がゼロなら()，
左人差し指で押す力()は、右親指で押す力(
)の2倍くらいでいいはず．
しかし，右小指に力が入ってしまうと(>0)，その力の2倍くらい，左人差し指で押す力を強くしないといけなくなる．
右手は，支点（唇）から離れている分，少しの力で大きく楽器を前方に回転させるモーメントをかけてしまう．
それに対抗するために，左手の支えを強くする必要があり，左人差し指付け根に負荷をかけてしまう．
右小指に力が入りすぎないようにしないと．
力とモーメントのつり合いがとれるためのは，さまざまな値をとりうる．
いろいろ試して，自分にとって一番いい指の位置や力加減を探していこう．
2次元圧縮性流体計算（無反射境界条件）
テキスト「空力音の直接数値計算を支える技術」（暗黙の型宣言著）の補足課題7
全部読んだらやってみよう！その2
課題
補足課題6で作成した圧縮性流体計算プログラムに、無反射境界条件を適用する．
- ヒント
- 特性波の方向に注意．入ってくる特性波を0にして打ち消す．
解答例
プログラム全体は載せられないので，ポイントだけ解説します．
- main.f90
TVDRK : block ! : ! 略 ! : !--- 時間積分の前にNSCBCを適用 ! NSCBC x-dir bd = ([1, Nx]) do j = 1, Ny ! x=0側の境界． !入ってくる特性波（つまり流速が正の特性波）を0にして打ち消す． L11(1,j) = max((u(1,j) - C0),0.d0)*(d_p_dx(1,j) - dens(1,j)*C0*d_u_dx(1,j)) ! L11の特性波の流速は(u - C0) L12(1,j) = max( u(1,j) ,0.d0)*((C0**2)*d_dens_dx(1,j) - d_p_dx(1,j)) ! L12の特性波の流速はu L13(1,j) = max( u(1,j) ,0.d0)*d_v_dx(1,j) ! L13の特性波の流速はu L15(1,j) = max((u(1,j) + C0),0.d0)*(d_p_dx(1,j) + dens(1,j)*C0*d_u_dx(1,j)) ! L15の特性波の流速は(u + C0) ! x=Lx側の境界． ! 入ってくる特性波（つまり流速が負の特性波）を0にして打ち消す． L11(2,j) = min((u(Nx,j) - C0),0.d0)*(d_p_dx(Nx,j) - dens(Nx,j)*C0*d_u_dx(Nx,j)) L12(2,j) = min( u(Nx,j) ,0.d0)*((C0**2)*d_dens_dx(Nx,j) - d_p_dx(Nx,j)) L13(2,j) = min( u(Nx,j) ,0.d0)*d_v_dx(Nx,j) L15(2,j) = min((u(Nx,j) + C0),0.d0)*(d_p_dx(Nx,j) + dens(Nx,j)*C0*d_u_dx(Nx,j)) end do do nb = 1, 2 d11 = (L12(nb,:) + 0.5d0*(L15(nb,:) + L11(nb,:)))/C0**2 d12 = (L15(nb,:) + L11(nb,:))*0.5d0 d13 = (L15(nb,:) - L11(nb,:))*0.5d0/(dens(bd(nb),:)*C0) d14 = L13(nb,:) ! NSCBCを適用した支配方程式(テキストの式(6.1)~(6.5))と， ! もともとの支配方程式(テキストの式(2.1)~(2.5))を比較し， ! 境界の値にNSCBCを適用する． d_dens_u_dx(bd(nb),:) = d11(:) d_dens_u_u_dx(bd(nb),:) = u(bd(nb),:)*d11(:) + dens(bd(nb),:)*d13(:) d_p_dx(bd(nb), :) = 0.d0 d_visyx_dx(bd(nb), :) = 0.d0 d_dens_u_v_dx(bd(nb),:) = v(bd(nb),:)*d11(:) + dens(bd(nb),:)*d14(:) d_dens_u_TE_dx(bd(nb),:)= 0.5d0*(u(bd(nb),:)**2 + v(bd(nb),:)**2)*d11(:) + d12(:)/(Gamma-1.d0) & + dens(bd(nb),:)*u(bd(nb),:)*d13(:) + dens(bd(nb),:)*v(bd(nb),:)*d14(:) d_u_p_dx(bd(nb),:) = 0.d0 d_qx_dx(bd(nb),:) = 0.d0 end do ! NSCBC y-dir ! : ! 略 ! : !--- Time Integration ! : ! 略 ! : end block TVDRK
計算結果
NSCBCを適用しても，壁からの非物理的な波の反射を完全には抑えられない． 特にdilatation分布を見ると，そのことがよくわかる．
局所境界条件であるNSCBCだけでは反射を完全に抑えられないので，適宜，領域境界条件を適用する． これはテキスト6.3でも述べられていることだけど，実際にやってみるとそのことがよく理解できる．
流速分布
図1 流速分布 渦度分布
図2 渦度分布 Dilatation分布
t>3/4で，四隅から波が反射して返ってきている．
図3 Dilatation分布
2次元圧縮性流体計算（単極渦）
テキスト「空力音の直接数値計算を支える技術」（暗黙の型宣言著）の補足課題6
全部読んだらやってみよう！その1
課題
いよいよ圧縮性流体計算！
8.2章に書かれている例題（超音速一様流で移流する単極渦まわりの流れ場）に取り組もう．
- ヒント
解答例
空間微分，時間積分はこれまでの補足課題と同様に以下の条件とします．
プログラム全体は載せられないので，ポイントだけ解説します．
- main.f90
do n = 1, Nt TVDRK : block integer :: i, j dens_u = dens*u dens_v = dens*v dens_TE = dens*TE ! TVDRKで時間進行させる保存量をキャッシュに入れる call cacheTVDRK_2D_NS(dens, dens_u, dens_v, dens_TE) do step=1,3 do j=1, Ny ! ひたすら空間微分 call compute1stDerivative(d_dens_u_dx(:,j), dens(:,j)*u(:,j), dx, Gx) call compute1stDerivative(d_dens_u_u_dx(:,j), dens(:,j)*u(:,j)*u(:,j), dx, Gx) call compute1stDerivative(d_p_dx(:,j), p(:,j), dx, Gx) call compute1stDerivative(d_dens_u_v_dx(:,j), dens(:,j)*u(:,j)*v(:,j), dx, Gx) call compute1stDerivative(d_dens_u_TE_dx(:,j), dens(:,j)*u(:,j)*TE(:,j), dx, Gx) call compute1stDerivative(d_u_p_dx(:,j), u(:,j)*p(:,j), dx, Gx) call compute1stDerivative(d_temp_dx(:,j), temp(:,j), dx, Gx) kp(:,j) = mu(:,j)*Cp/Pr qx(:,j) = -kp(:,j)*d_temp_dx(:,j) call compute1stDerivative(d_qx_dx(:,j), qx(:,j), dx, Gx) call compute1stDerivative(d_u_dx(:,j), u(:,j), dx, Gx) call compute1stDerivative(d_v_dx(:,j), v(:,j), dx, Gx) end do do i=1, Nx call compute1stDerivative(d_dens_v_dy(i,:), dens(i,:)*v(i,:), dy, Gy) ! : ! 省略 ! : end do ! 粘性応力は，x方向とy方向両方の微分から決まるため，ひととおりx方向とy方向の微分値を求めた後に計算． visxx = mu*(2.d0*d_u_dx - (2.d0/3.d0)*(d_u_dx + d_v_dy)) visxy = mu*(d_u_dy + d_v_dx) visyx = mu*(d_v_dx + d_u_dy) visyy = mu*(2.d0*d_v_dy - (2.d0/3.d0)*(d_u_dx + d_v_dy)) ! 粘性応力の微分．もっと厳密にやるなら，テキスト8.1章「2階微分の計算法」に書かれているやり方でやりましょう． do j=1, Ny call compute1stDerivative(d_visxx_dx(:,j), visxx(:,j), dx, Gx) call compute1stDerivative(d_visyx_dx(:,j), visyx(:,j), dx, Gx) call compute1stDerivative(d_U_visx_dx(:,j), u(:,j)*visxx(:,j) + v(:,j)*visyx(:,j), dx, Gx) end do do i=1, Nx ! : ! 省略 ! : end do ! 時間積分 call integrateTVDRK_2D_NS(dens, -(d_dens_u_dx + d_dens_v_dy), dt, step, 1) !OK call integrateTVDRK_2D_NS(dens_u, & -(d_dens_u_u_dx + d_dens_v_u_dy) - d_p_dx + d_visxx_dx + d_visxy_dy, & dt, step, 2) call integrateTVDRK_2D_NS(dens_v, & -(d_dens_u_v_dx + d_dens_v_v_dy) - d_p_dy + d_visyx_dx + d_visyy_dy, & dt, step, 3) call integrateTVDRK_2D_NS(dens_TE, & -(d_dens_u_TE_dx + d_dens_v_TE_dy + d_u_p_dx + d_v_p_dy) & + d_U_visx_dx + d_U_visy_dy - (d_qx_dx + d_qy_dy), dt, step, 4) ! 各物理量を計算．これを次の時刻で空間微分する． u = dens_u/dens v = dens_v/dens TE = dens_TE/dens dens_e = dens_TE - 0.5d0*(dens*u*u + dens*v*v) e = dens_e/dens p = dens_e*(Gamma-1.d0) temp = e/Cv mu = Mu0*((Temp0 + Temp1)/(temp + Temp1))*((temp/Temp0)**(3/2)) enddo end block TVDRK end do
- module_TimeIntegration.f90
! TVDRKで時間進行させる保存量をキャッシュに入れる subroutine cacheTVDRK_2D_NS(dens, dens_u, dens_v, dens_TE) real(8), intent(in) :: dens(:,:), dens_u(:,:), dens_v(:,:), dens_TE(:,:) dens2(:,:) = dens(:,:) dens_u2(:,:) = dens_u(:,:) dens_v2(:,:) = dens_v(:,:) dens_TE2(:,:) = dens_TE(:,:) end subroutine cacheTVDRK_2D_NS ! TVDRKでの時間積分は，こんなふうに書いてみた． ! けっきょくどの物理量も式①で積分することは同じなので，コードの繰り返しを避けるため，被積分配列をcase文で選択． subroutine integrateTVDRK_2D_NS(f, d_f_dt, dt, step, n) implicit none real(8), intent(inout) :: f(:,:) real(8), intent(in ) :: d_f_dt(:,:) real(8), intent(in ) :: dt integer, intent(in ) :: step integer, intent(in ) :: n real(8), parameter :: a(3) = [0d0, 3d0/4d0, 1d0/3d0] real(8), parameter :: b(3) = [1d0, 1d0/4d0, 2d0/3d0] select case (n) case (1) fn2 = dens2 case (2) fn2 = dens_u2 case (3) fn2 = dens_v2 case (4) fn2 = dens_TE2 end select f(:,:) = a(step)*fn2(:,:) + b(step)*(f(:,:)-dt*d_f_dt(:,:)) ! 式① end subroutine integrateTVDRK_2D_NS
計算結果
壁面に無反射境界条件を適用していないので，壁付近まで渦が広がっていったところ(t=4/8くらい)で，壁から非物理的な波が反射して返ってきます． その影響でt=5/8を過ぎると発散．
流速分布
図1 流速分布 渦度分布
図2 渦度分布 Dilatation分布
図3 Dilatation分布
この課題では，ここまでできればOK．次の課題で無反射境界条件に取り組もう！
２次元拡散方程式
テキスト「空力音の直接数値計算を支える技術」（暗黙の型宣言著）の補足課題5
7章まで読んだらやってみよう！その3
本課題も別冊テキスト「流体計算で覚えるPython3」（暗黙の型宣言著）に掲載されているテスト問題を使用します．
課題
2次元拡散方程式をFortranで解いてみよう．
- 支配方程式は，２次元拡散方程式．
- 初期条件は，2次元ベル型関数．（別冊テキスト7.7節と同様）
- ヒント
- 補足課題4「2次元移流方程式」で作成したプログラム一式をベースに，支配方程式を拡散方程式に変更すればOK．
解答例
プログラム全体は載せられないので，ポイントだけ解説します．
拡散方程式では，空間微分が1階微分ではなく２階微分になるので，これを計算するコードに修正． さらに，別冊テキストの離散化式(7.9)をTVD Runge-Kutta法で時間積分するコードを追加．
２次元移流方程式から２次元拡散方程式への変更はこれくらい．一瞬でできてしまう．
- main.f90
do n = 1, Nt TVDRK : block integer :: i, j call cacheTVDRK_2D(f%var) do step=1,3 do j=2, Ny-1 call compute2ndDerivative(f%dx2(:,j), f%var(:,j), dx, Gx) ! x方向の２階微分を計算 end do do i=2, Nx-1 call compute2ndDerivative(f%dy2(i,:), f%var(i,:), dy, Gy) ! y方向の２階微分を計算 end do ! 境界でゼロを強制 f%var(1 ,: ) = 0d0 f%var(Nx,: ) = 0d0 f%var(: ,1 ) = 0d0 f%var(: ,Ny) = 0d0 call integrateTVDRK_2D_DiffusionEq(f%var, f%dx2, f%dy2, dt, step, mu) ! TVDRK法で拡散方程式の時間積分するサブルーチンを呼び出す enddo end block TVDRK end do
ここまでの補足課題をこなしていれば，サブルーチン「compute2ndDerivative」は module_CompactScheme.f90 で既に定義してあるはず． もし定義していなければ，テキストを参考に定義してください．
- module_TimeIntegration.f90
subroutine integrateTVDRK_2D_DiffusionEq(f, d_f_dx2, d_f_dy2, dt, step, mu) implicit none real(8), intent(inout) :: f(:,:) real(8), intent(in ) :: d_f_dx2(:,:), d_f_dy2(:,:) real(8), intent(in ) :: dt integer, intent(in ) :: step real(8), intent(in ) :: mu real(8), parameter :: a(3) = [0d0, 3d0/4d0, 1d0/3d0] real(8), parameter :: b(3) = [1d0, 1d0/4d0, 2d0/3d0] integer :: Nx, Ny integer :: i, j Nx = size(f,1) Ny = size(f,2) do j=1, Ny do i=1, Nx f(i,j) = a(step)*fn2(i,j) + b(step)*(f(i,j) + mu*dt*(d_f_dx2(i,j) + d_f_dy2(i,j))) ! TVDRK法で拡散方程式の時間積分 end do end do end subroutine integrateTVDRK_2D_DiffusionEq
拡散数0.25の計算結果は図1のとおり．時間とともに円錐が拡散している． もっとしっかりやるなら，厳密解（２次元ガウス分布）との誤差も計算するとよいでしょう．
２次元移流方程式
テキスト「空力音の直接数値計算を支える技術」（暗黙の型宣言著）の補足課題4
7章まで読んだらやってみよう！その2
今回は２次元の課題にトライ．
8章に2次元のテスト問題が書かれているけれど，いきなりナビエ・ストークス方程式を解き始めるよりは，もう少し単純な課題で２次元計算の基礎を知っておいたほうがいいと思う．
なお，本課題では別冊テキスト「流体計算で覚えるPython3」（暗黙の型宣言著）に掲載されているテスト問題を使用します．
課題
「流体計算で覚えるPython3」７章に書かれている回転円錐問題をFortranで解いてみよう．
- 支配方程式は，２次元移流方程式．（7.1節に掲載）
- 初期条件は，2次元ベル型関数．（7.5節に掲載）
- ヒント
- 補足課題3「時間積分」で作成したプログラム一式をベースに２次元化すればOK．
解答例
プログラム全体は載せられないので，ポイントだけ解説します．
まずはmain.f90のInitializeブロックで初期条件を定義． 移流速度c_xはy，c_yはxで決まる，という点に注意．うっかり逆にしないこと！
Initialize : block : c_x(i,j) = -2.d0*PI*y(j) c_y(i,j) = 2.d0*PI*x(i) : call fileout_2dim(file_init, x, y, c_x, c_y, f%var) !初期条件出力 end block Initialize
初期条件を定義できたら，いったん出力してみて確認．
図1のt/T=0のようになっていればOK．
次に，２次元移流方程式を解くためのコードを書く．ここでは重要なポイントだけ紹介．
- main.f90
df/dx, df/dyをそれぞれ計算し，fを求める．というのを繰り返し，時間積分していく．
do n = 1, Nt TVDRK : block integer :: i, j call cacheTVDRK_2D(f%var) do step=1,3 do j=2, Ny-1 call compute1stDerivative(f%dx(:,j), f%var(:,j), dx, Gx) end do do i=2, Nx-1 call compute1stDerivative(f%dy(i,:), f%var(i,:), dy, Gy) end do call integrateTVDRK_2D_Adv(f%var, f%dx, f%dy, dt, step, c_x, c_y) ! 境界でゼロを強制 f%var(1 ,: ) = 0d0 f%var(Nx,: ) = 0d0 f%var(: ,1 ) = 0d0 f%var(: ,Ny) = 0d0 enddo end block TVDRK end do
- module_TimeIntegration.f90
２次元移流方程式を時間方向に離散化して，次ステップのfを求める． ここではTVD Runge-Kutta法で時間積分．
subroutine integrateTVDRK_2D_Adv(f, d_f_dx, d_f_dy, dt, step, conv_x, conv_y) implicit none real(8), intent(inout) :: f(:,:) real(8), intent(in ) :: d_f_dx(:,:), d_f_dy(:,:) real(8), intent(in ) :: dt integer, intent(in ) :: step real(8), intent(in ) :: conv_x(:,:), conv_y(:,:) real(8), parameter :: a(3) = [0d0, 3d0/4d0, 1d0/3d0] real(8), parameter :: b(3) = [1d0, 1d0/4d0, 2d0/3d0] integer :: Nx, Ny integer :: i, j Nx = size(f,1) Ny = size(f,2) do j=1, Ny do i=1, Nx f(i,j) = a(step)*fn2(i,j) + b(step)*(f(i,j)-conv_x(i,j)*dt*d_f_dx(i,j)-conv_y(i,j)*dt*d_f_dy(i,j)) end do end do end subroutine integrateTVDRK_2D
ひととおりプログラムができたら，円錐が1/4周，2/4周，3/4周，1周したときの計算結果を出力し，正しく反時計回りにまわっていることを確認． 可視化ソフトは「AVえせ」というフリーウェアを使用．# あやしいソフトじゃないよ！
格子点数やクーラン数を変えてみて，1周したときの計算結果を可視化して確認する．
まず，総格子点数を変えてみた場合．可視化結果にほとんどちがいはみられない．（格子点数の分だけベクトルが表示される）
次に，クーラン数(CFL)を変えてみた場合．こちらも，ほとんどちがいはみられない．
1周したときの領域全体のエラーを出力する．
print *, "CFL:", CFL, "error:", sum(abs(f%var-f0%var))/(Nx*Ny) ! f%varは計算結果．f0%varは初期値．
格子点数とクーラン数を変えたときのエラーをグラフにすると，
クーラン数よりも格子点数によってエラーが決まっている．
時間積分
テキスト「空力音の直接数値計算を支える技術」（暗黙の型宣言著）の補足課題3
7章まで読んだらやってみよう！
課題
(1) 時間積分にEuler法を用いて移流方程式を計算する．
- 初期値はテキスト式(7.11)とし，波が境界を出ない時間まで計算する．
- 境界では勾配ゼロ(df/dx=0)の境界条件を与える．
- 誤差も計算して，実効的な精度を確認する．
- テキストの以下プログラム参照．
- List7.2(移流方程式メインルーチン)
- List7.3(Euler法による時間積分)
- 空間微分は，補足課題2を参考に，以下を用いる．
- グローバルスキーム・・・6次精度
- 境界1点内側・・・5次精度
- 境界・・・4次精度片側差分
(2) 同様に，TVD Runge-Kutta法（以下，TVDRK）を用いて移流方程式を計算する．
- List7.14(TVDRK法による時間積分)．
解答例
まずはEuler法．テキストを写経して，以下のファイルを作成する．
- main.f90(List7.3参照)
- module_TimeIntegration.f90(List7.4参照)
main.f90を以下のように修正する．（修正箇所は「!---」, ただのコメントは「!」）
real(8), parameter :: EndT = 10d0 !--- 計算時間を，波が境界に到達しない時間に設定 : !--- 境界条件用パラメータを追加 real(8), parameter :: df1 = 0d0 real(8), parameter :: dfN = 0d0 !--- 3重対角行列を指定 type(TriDiagonalMatrix) :: G : Initialize : block : !--- 両端で勾配ゼロになるようオフセットする(3.8章参照) f%var(1) = f%var(2) f%var(Nx)= f%var(Nx-1) end block Initialize : do n = 0, Nt-1 Euler : block call compute1stDerivative(f%dx, f%var, dx, G, df1, dfN) !--- 境界条件用パラメータdf1, dfNを追加 call integrateEuler(f%var, f%dx, dt, conv) end block Euler end do !--- 誤差計算・出力用ブロックを追加 ResultOut : block integer, parameter :: Nx_err = nint(wavelength/dx) ! 波長を誤差計算範囲とする integer :: i, i_wavecenter_f, i_wavecenter_f0 real :: err err = 0d0 i_wavecenter_f0= nint(Nx/2d0) + 1 ! 初期の波f0の中心位置 i_wavecenter_f = i_wavecenter_f0 + nint(conv*dt*(Nt-1)/dx) ! 伝播した波fの中心位置 open(1, file='FileOut.txt', status='replace') do i=1,Nx write(1,*) dble(i-1)*dx - Lx/2d0, f0%var(i), f%var(i) ! x, f0, f を出力 enddo do i=-Nx_err, Nx_err err = err + abs(f%var(i_wavecenter_f+i)-f0%var(i_wavecenter_f0+i)) ! 初期の波と伝播した波の振幅差を計算 enddo print *, "MinAveError:", err/(Nx_err*2+1) ! 累積誤差を格子点数で割って，最小平均誤差を出力 close(1) end block ResultOut : contains subroutine compute1stDerivative(d_f_dx, f, dx, G, GraditntB1, GraditntBN) !--- 境界条件用パラメータを追加 use TypeTriDiagonalMatrix !--- 3重対角行列を指定 implicit none : ! 境界条件用パラメータを追加 real(8), intent(in ), optional :: GraditntB1 real(8), intent(in ), optional :: GraditntBN call construct1stDerivativeVector(d_f_dx, f, dx, G, GraditntB1, GraditntBN) !--- 境界条件用パラメータを追加 call solve(G, d_f_dx) end subroutine compute1stDerivative end program main
- module_CompactScheme.f90(List3.4参照)
- TypeTriDiagonalMatrix.f90(List4.2参照)
以上でEuler法のプログラムは完成．
パラメータ「CFL」をさまざまに変化させて，誤差の変化を確認する．確認結果は後述．
次に，TVDRK法．main.f90のEuler法用ブロックを，TVDRK法用ブロックに変更．
- main.f90
:
do n = 0, Nt-1
TVDRK : block
call cacheTVDRK(f%var) ! fnの値を保持
do step=1,3
call compute1stDerivative(f%dx, f%var, dx, G, df1, dfN)
call integrateTVDRK(f%var, f%dx, dt, step, conv) ! TVDRKのサブルーチンをcall
enddo
end block TVDRK
end do
:
さらに，module_TimeIntegration.f90にテキストList7.14のサブルーチンcacheTVDRKとintegrateTVDRKを追加する．
テキストにも書いてあるように，Euler法からTVDRK法への変更は簡単♪
TVDRK法は，Euler法から得られた結果と現時間ステップの結果に重みづけして，次の時間ステップの値を決めている，と解釈できる．
パラメータ「CFL」をさまざまに変化させて，誤差の変化を確認した結果は図1のとおり．
Euler法はslope 2（CFLを1/10にすると誤差が1/100）よりも傾きが浅い．近似直線を引いてみると，約1.1次精度であった．
TVDRK法(Nx = 400)は，0.1<CFL<0.6でほぼslope3に沿っているので約3次精度． 0.6<CFLでslope3より誤差が大きいのは，格子が荒すぎる影響．
CFL<0.1で誤差の減少がサチっているのは？
Nx=400では0.00001付近でサチっているのに対し，格子を2倍荒くしたNx=200では誤差が0.001付近でサチっている．
CFL=0.002では，Nx=400の誤差は，Nx=200の誤差の約倍であった．
空間微分のグローバルスキームが6次精度なので，格子を2倍にすることで誤差が
倍に減少している．
つまり，誤差の減少がサチるのは，空間微分の精度の影響．
