目 次
第 1 章 常微分方程式 ut= au の差分解法 3 1.1 一階微分の差分 (finite difference) 近似 . . . . 3 1.2 ut= au 差分方程式とその解法 . . . . 3 1.3 振動方程式・前方差分近似における Von Neumann の安定性解析 . . . . 5 1.4 振動方程式・中央差分近似における Von Neumann の安定性解析 . . . . 6 1.5 振動方程式・中央差分近似での数値計算 . . . . 6 1.6 時間に中央差分を用いる解法 (Leap Frog) . . . . 7 1.7 Leap Frog における，物理モードと計算機モード . . . . 8 1.8 離散化された微分の打ち切り誤差 . . . . 9 第 2 章 拡散方程式 Ut= aUxxの差分解法 10 2.1 拡散方程式の差分方程式 . . . . 10 2.2 拡散方程式についての Neumann の安定性解析 . . . . 10 2.3 前方差分による，拡散方程式の数値計算 . . . . 11 第 3 章 移流方程式 ut=−cux 13 3.1 移流方程式の差分式 . . . . 13 3.2 移流方程式の Neumann の安定性解析 . . . . 133.3 移流方程式の境界条件（radiative boundary condition と cyclic boundary condition) . . . . 14
3.4 移流方程式の数値計算 . . . . 15 3.5 CFL condition の物理的解釈 . . . . 15 3.6 時間に前方・風上（空間後方）差分の影響領域 . . . . 15 3.7 風上差分における人工粘性 . . . . 16 3.8 分散性 −空間差分の影響− . . . . 16 3.9 移流拡散方程式 . . . . 17 第 4 章 波動方程式 utt= c2uxx 19 4.1 波動方程式の解釈 . . . . 19 4.2 波動方程式の分解と差分方程式化と数値計算 . . . . 19 4.3 波動方程式の数値計算 . . . . 20 4.4 互い違いグリッド（Staggard Grid） . . . . 20
海洋各論
2
見延 庄士郎、渡邊 祥史
扱う微分方程式
この授業で扱う微分方程式は，大きく分けると，初期値問題と境界値問題に分けることが出来る．まず は，初期値問題から学んでいこう．初期値問題は，もちろん時間微分を含んだ方程式で表わされ，そのよう な方程式は特に時間発展（偏）微分方程式とも呼ばれる．時間発展微分方程式にも，様々なものがあるが， まずは時間微分項が一回の微分となっているものを対象としよう．簡単のために，空間は一次元・従属変数 は一つとする． 簡単な空間一変数の一次元の時間発展線形偏微分方程式は，次のような形をしている． ut= a + bu + cux+ duxx (a, b, c, d は定数) (1) 今，変数 u は時間 t と空間上の位置 x の関数である．もちろん右辺第三あるいは第四項がなければ，常微 分方程式といわれる．まずはこの各項を順々にやっていこうと思う．ただし，あまりにもおもしろくない ut = a は，それだけを特に取り上げることはしない． 上記の，時間微分が一階である，空間一次元・一変数の微分方程式で数値計算の基本的な取り扱いを十分 に学んだ後に，波動方程式，あるいは空間二次元，複数変数の微分方程式の数値計算について学ぶ．そのの ち，境界値問題を扱おうと思う． 今後取り扱う方程式の一覧を，次に挙げておく． ut= bu b が実数ならば指数的増加・減衰の式 b が純虚数ならば振動の式 ut= cux 移流方程式 utt= c2uxx 波動方程式 uxx+ uyy= f (x, y) ポアソン方程式 (境界値問題)第
1
章 常微分方程式
u
t
= au
の差分解法
1.1
一階微分の差分
(finite difference)
近似
この授業では一貫して，差分法で微分方程式を解いていく．微分方程式で表される従属変数 (例えば u) は 一般に連続な独立変数 (例えば t) 上で定義されている．差分法は，離散的な格子点 (grid, mesh その他もろ もろの呼び方をされる場合がある) のみで定義された変数の値を，もとの微分方程式の近似方程式である差 分方程式によって求める方法である． 今，u をある離散的な時間 (t = ∆t× n, n = 1, 2, 3, · · ·) 上で定義されたものとして，微分の近似表現を与 えよう．微分の定義は， du dt(t) = δtlim→0 u(t + δt)− u(t) δt (1.1) であるから，(t = ∆t× n) 上のみで定義されている u については，次のような近似表現が可能である． du dt(t) un+1− un ∆t (1.2) ここで，unは，t = ∆t× n における u を表している．また，(1.1) の代わりに du dt(t) = δtlim→0 u(t)− u(t − δt) δt (1.3) あるいは， du dt(t) = δtlim→0 u(t + δt)− u(t − δt) δt (1.4) と表しても，δt → 0 の極限においては同じである．しかし，(1.3) あるいは (1.4) に対応する差分近似は (1.1) に対する差分近似と同じでは無い．従って，一階微分の差分近似は次の三つが考えられる． 前方差分 du/dt|t=n∆t = (un+1− un)/∆t 中央差分 du/dt|t=n∆t = (un+1− un−1)/2∆t (1.5) 後方差分 du/dt|t=n∆t = (un− un−1)/∆t (1.5) の右辺は ∆t を小さくして行けばもちろんすべて du/dt になるが，実際には ∆t は有限の大きさを持つ ので，下の図のように各々違いが出る． これらの差分近似の中でどれを実際の数値解に用いるべきかは，方程式に応じて定まる．この授業では 後に学 Neumann の安定性解析などで，それを論ずる．1.2
u
t= au
差分方程式とその解法
手始めに，まず常微分方程式 ut= au (1.6)の数値計算をやってみよう．この式は良く見る式である．定数 a が実数であれば，その正負に従って，指 数関数的に増減する式だし，虚数ならば，振動の式である．この式を解くためには，この式のほかに，初期 条件 u(t = 0) が必要なことに注意して欲しい．解は
u(t) = exp (a× t) × u(t = 0) (1.7)
となる． (1.6) 式の左辺を前方差分で表すこととして，(1.6) 式を離散化すると， (un+1− un)/∆t = aun (1.8) と書き表わせる．このように離散化された方程式をこれを差分方程式と呼ぶ．今，un は第 n ステップ (t = ∆t× n) まで既知であるとすると，未知変数 un+1は， un+1 = ∆taun+ un (1.9) で求めることが出来る．これを，n = 1, 2, 3, 4, と繰り返していくのが，時間発展常微分方程式を数値的に 求める (数値的に解く) ことになる．では，実際に数値計算をやってみよう．次に，例題プログラムを載せ るので，ASPEN で入力してもらう．チェックの便利のために，このプログラムに限らず，メンバ名の英字 部は EX でなくとも良いが，数字部は必ずこちらが提示する番号を付けること． これは，全部で 65 行のプログラムである．ここまで長くなったのは，絵を描くとか，解析解を求めるとかし ているからで，本当にこのプログラムの中で数値計算をしているのは，100 番の DO LOOP と SUBROUTINE のホンの数行である． このプログラムでは，数値計算の解 unに相当する変数は，U1 と U2 の二つだけである．配列は数値計 算のためではなく，絵を描くためにとられている．従って，たとい 100, 1000 の時間ステップにおける u の 値を求めるにも，計算に必要なメモリー量は U1 と U2 のためのわずか 2 変数分でしかない．n = 1, 3, 5, ... の，奇数ステップでは，U1 から U2 を求め，偶数ステップでは，U2 から U1 を求めている．念のために図 示すると，下のようになる． 時間ステップ (n) : 1 2 3 4 5 MAIN における変数 : U1 U2 U1 U2 U1 即ち，前方差分による数値計算を行なうには直前の結果のみ分かっていれば良い．また，MOD は剰余を与 える FORTRAN の関数で，奇数ステップか偶数ステップかを判断している．n=1 の値 (初期値) を U1 に代 入しているために，上図の順番になっている．もし，初期値を U2 に代入したなら，上図の U1 と U2 が逆 にならなくてはならない． 数値計算の結果は絵を描くために，DIMR(実数部) と DIMI(虚数部) の二つの配列に格納されている．こ こで，REAL と IMAG は，それぞれ複素数の実数・虚数部分を取り出す FORTRAN の関数である．また， 100 番の DO LOOP 中の GOTO 文は，初期値の 100 倍 (THRESH) よりも数値計算の値が大きくなってし まったら，そこで計算を打ち切るために使われている．
MAIN プログラム (SUBROUTINE ではないプログラム) での，CALL OPENPG 以下の CALL 文は全 て，こちらで用意した描画ルーチンを呼んでいる．資料４に，ここで使用されている描画用サブルーチンの 簡単な説明と，そのソース (プログラム・リストのこと) を載せた．
上の EX1 が入力できたのなら，READY が表示されている状態で， GO EX1
と入力してみよう．すると EX1 が起動され， ’Input, number of points’
’Input, time step’
’Input, real part of constant’ ’Input, imaginary part of constant’
という入力要請が順々に出る．最初のは計算するステップの数，t = n∆t としたときの n の最大値 (NMAX) である．次のは，∆t，最後の二つは，(1.6) 式の定数 a の実部と虚部である．とりあえず， (CASE あ) 1000 最大ステップ数 0.01 ∆t -0.1 定数実部 0 定数虚部 と打ち込んでみよう．これは，実数定数だから指数関数的減衰であるが，あまりおもしろくない．まあ試し にやってみただけである．ちなみに，数値解はペン番号１番 (画面上で白線，紙の上では細線)，解析解は ペン番号２番 (画面上で黄色線，紙では中太線) で表わされている． 次に， (CASE い) 1000 最大ステップ数 0.1 ∆t 0 定数実部 1 定数虚部 と打ち込んでみよう．これは，定数が純虚数なので，振動を表わしている．しかし，計算結果が爆発してい ることがよく分かると思う．このサンプル・プログラムではこのような爆発に耐え得るよう，ある一定値以 上に数値解が大きくなると計算をやめるようになっているが，そのような処置がとられていない場合は，” 指数オーバーフロー例外 ”のエラーが出るはずである． というわけで，差分で時間発展偏微分方程式を一応差分方程式で書くのは簡単だが，それだけでは安定に 計算されないことが分かってもらえたと思う．
1.3
振動方程式・前方差分近似における
Von Neumann
の安定性解析
数値計算の安定性を，求める方法はいろいろあるが，最も実用的で広く使われているのは Von Neumann の安定性解析である．この授業では最後まで愛用する．この方法では， u(n + 1) = λu(n) (1.10) と置いて λ を求め，λ は一般に複素数なのでその絶対値の大きさを評価する．つまり， |λ| < 1 安定 (減衰) |λ| = 1 安定 (中立) |λ| > 1 不安定 (増大) と判断する．外力がある場合には，外力項を落とした方程式で評価する．安定であるためには|λ| ≤ 1 であ る．(1.9) 式と上の (1.10) 式を一つにして，a が純虚数 iω の場合を考えよう．ここで ω としたのは，a の 虚部は角周波数になり，角周波数は慣例として ω・σ などで表記することになっているからである．そうす ると，λu(n) = i∆tωu(n) + u(n) (1.11)
λ = 1 + i∆tω (1.12)
|λ|2= 1 + (∆tω)2> 1 (1.13)
1.4
振動方程式・中央差分近似における
Von Neumann
の安定性解析
前節では振動方程式は，前方差分を用いるなら，必ず不安定になることを示した．それでは，他の差分近 似を用いたら安定に計算できるであろうか．一般に時間発展方程式は前方差分か中央差分で時間微分項を 近似することが多い．ここでは，中央差分近似で数値的に安定に計算できるかどうかを調べよう． 振動方程式を中央差分近似を用いて書くと， (un+1− un−1)/(2∆t) = iωun (1.14) 従って，un+1について解くと， un+1= 2i∆tωun+ un−1 (1.15)この式において，un+1= λ2un−1，un= λun−1として，共通因数 un−1を辺々割って消去すると， λ2− 2i∆tωλ − 1 = 0 (1.16) なる式が得られる．ここで上式を，解の公式で解くと， λ = i∆tω±√1− (∆tω)2 (1.17) となる．安定・不安定は|λ| の大きさできまる．これは，ルートの中が正なら， |λ| = (∆tω)2+ 1− (∆tω)2= 1 (1.18) であるから．ルートの中が正，つまり，∆tω≤ 1 ならば，Leap Frog は中立 (安定) である．1-2 節で調べた 前方差分の場合は|λ|2= 1 + (∆tω)2なので，∆t に応じて不安定の度合 (成長率) が決まるが，Leap Frog では ∆t≤ 1/ω でありさえすれば，∆t によらず中立である． 逆に ∆tω が 1 より大きいなら， |λ| = ∆tω +√(∆tω)2− 1 (1.19) (∆tω > 1, (∆tω)2− 1 > 0) (1.20) なので，|λ| > 1 であり，Leap Frog は必ず不安定になる．
1.5
振動方程式・中央差分近似での数値計算
前節では，時間微分を中央差分で近似するなら，安定に数値計算可能であることが分かった．早速，実際 に数値計算を行なって確かめてみよう．(1.6) 式に中央差分近似を用いて得られる un+1 = 2a∆tun+ un−1 (1.21) が，数値計算されるべき式である． ここで注意しなくてはならないのが，初期値である．中央差分では，n+1 ステップ目を求めるためには， 第 n ステップと n-1 ステップが必要なので，n=2 ステップ目の計算は，中央差分ではなし得ない．従って， n=2 ステップの値を求めるのは，EX1 で用いた前方差分によることにしよう． ではさっそくプログラムを書いてみよう．先週作ったプログラムを流用するのが手っとり早いので，ま ず，READY の出ている状態で，COPY NUMERCL2.FORT(EX1) NUMERCL2.FORT(EX2) として，EX1 を EX2 にコピーしよう．さて，ここで問題です．
EX.2 EX1 を元に，時間微分項を中央差分近似した数値計算のプログラムを作成せよ． 参考：この場合でも，U1 と U2 のみを用いて，数値計算を行なう．中央差分の数値計算を行なう SUBROUTINE は，例えば，
SUBROUTINE CENTER(UNXT, UNOW, CNST, DT)
のように変数を引き渡して，入力された UNXT を un−1，UNOW を un，出力時の UNXT を un+1とすれ
ば良い．納得が行かない人のために，変数がどのように使われるべきかを図示しておこう． 時間ステップ : 1 2 3 4 5 6 MAIN における変数 : U1 U2 U1 U2 U1 U2 図の矢印は，例えば n=3 ステップにおける u の値を求めるためには，n = 1 および 2 ステップ目の u の 値から求めることを示している． 参考の参考：実は EX1 では，U1,U2 の二つの変数は必要なく，一つだけでも数値計算は可能であった． ただ，EX2 や今後の問題を解くのに二つ変数を取り扱うことに慣れておいた方があとがスムーズであろう と思ったので，あのような例題とした．
さて EX2 が出来たのならば，EX1 でみごとに解が爆発した CASE(い) 最大ステップ数 100 ∆t 0.1 定数実部 0 定数虚部 1 を入力してみよう．うまくいったであろう．このパラメータのうち，最大ステップ数だけを十倍の 1000 にし て試してみよう (計算する時間範囲が CASE(う) の十倍になる)．やはりうまく行っただろう．しかも，殆ど 解析解と数値解が重なっていることが，つまり，正しく数値計算がなされていることが分かるだろう．次に， (CASE う) 最大ステップ数 1000 ∆t 0.3 定数実部 0 定数虚部 1 を入力してみよう．数値解の振幅は変化しないが，解析解よりも周期がやや短いのが，注意して見るならば 分かるだろう．従って，安定に計算できても数値解が常に正しいとは限らない．さて，もうひとつ， (CASE え) 最大ステップ数 100 ∆t 1.01 定数実部 0 定数虚部 1 を，入力してみよう．この結果から，Neumann の安定性解析で示されたように，∆tω が１よりも大きいな ら，数値的計算は不安定であることが確認できる．
1.6
時間に中央差分を用いる解法
(Leap Frog)
一般の偏微分時間発展方程式 ∂u/∂t = f (x, t) を，時間に中央差分を用いて書くと， (un+1− un−1)/(2∆t) = fn (1.22) un+1 = (2∆t)fn+ un−1 (1.23)となる．この数値解の方法は Leap Frog (蛙飛び) とも呼ばれる．これは，下図のように n + 1 の値を求め るには，一つ前の unは使わずに，二つ前の un−1を用いるからである．蛙飛びの雰囲気が分かってもらえ るだろうか？ さて，振動方程式の数値計算においては，leap frog は前方差分よりもよい結果をもたらした．では，減 衰方程式 (ut=−au, a > 0) についてはどうであろうか？ 問題 Neumann の安定性解析を用いて，減衰方程式を中央差分近似した場合 および前方差分近似した場合に，安定に数値計算し得るかどうかを調 べよ．また，安定に数値計算し得るなら，安定条件を求めよ．得られた 結果は，Ex1 と Ex2 で確認せよ．
1.7
Leap Frog
における，物理モードと計算機モード
Leap Frog のような，3 レベル (n + 1 ステップを求めるのに，n と n− 1 を用いる．合わせて 3 レベル．) の解法では，かならず 2 つのモードが存在する．一つは物理的 (数学的に) 真の解の近似である物理モード で，もう一つは数値計算で出てくるだけの物理的には意味のない計算機モードである．実際に得られる数 値解はこの両モードの混合で与えられる． ここで，この世でいちばん簡単な常微分方程式 ∂u/∂t = 0 で計算機モードを考えて見よう．この数値解 は，1.6 節の図の特殊な場合 f (t) = 0 で，次の図のように求められる．ここでは，偶数ステップと奇数ス テップが何の関係も持たないことが分かるだろう．すなわち，解は U (奇数ステップ)=定数 1，U (偶数ステッ プ)=定数 2 となる．次に Neumann の安定性解析の視点から，計算機モードを考えよう．Leap Frog についての Neumann の 安定性解析において，安定・不安定の他に重要なのは，λ が二つの解を持つということである．結局のとこ ろ，安定性解析では unは Λn−1と表現されるのだから，一般に n レベルの解法では，n− 1 次方程式が得 られ，λ は n− 1 個の解を持つ．もちろんこの中の一つだけが物理モードで，他は計算機モードである． 中央差分近似された振動方程式の安定性解析では，λ は次のように求められた． λ = i∆tω±√1− (∆tω)2 ∆t→ 0 の極限では，λ は，ルートの正負に応じて， λ→ 1, −1 (1.24) の，二つの極限に近付く．もちろん ∆t→ 0 なのだから，λ → 1 が正しい解であり，λ → −1 は，1 ステッ プ毎に符号の変わる計算機モードである．なお，前方差分のような 2 レベルの解法では，計算機モードは 存在しない． EX3 下の説明に従ってプログラム Ex.3 を作成・実行せよ．
EX1 を EX2 にコピーしたのと同様に，まず EX2 を新しいメンバ EX3 にコピーする．EX3 では，1・2STEP 目の値の平均値を 1 とし，両 STEP の値の差を入力することにしよう．解析解の初期値は，数値解用の２ つの初期値の平均値である 1 とする． 話は少々ややこやしいが，要は下図の通りである． つまり， u1 = 真の初期値− 初期値の差/2 (1.25) u2 = 真の初期値 + 初期値の差/2 (1.26) である．従って，U1 と U2 は次のように与えればよい．
WRITE(6,*) ’Input, difference of initial value’ READ(5,*) DIFFER
U1 = DINIT - DIFFER/2 U2 = DINIT + DIFFER/2
また，DO 100 の DO LOOP は n=2 からはじめることと，DIMR(2) と DIMI(2) を U2 の値から求めてお くのを忘れないように． DEFFER には，0, 0.2, 0.5, 1, 2 あるいはこれらの負の値等を入力してみよう．すると，ギザギザの数値 解が得られるだろう．この数値解は，物理モードの解の上に，1 ステップ毎のギザギザの計算機モードが重 ね合わされたものである． 2 つのモードがどの程度の割合で発生するかは，線形方程式の場合完全に初期値によっているので，Ex2 のように第 2 ステップの値をある程度正しく求めるならば，計算機モードはその後の数値計算には悪影響 を及ぼさない．しかし，一般に非線形方程式では，初期値が計算機モードを発生させないように与えられて いた場合でも，計算を進めるにつれて計算機モードが発生してしまう．したがって，非線形方程式の数値計 算においては，20 ステップ程度に一回計算機モードを殺し物理モードを生かす 2 レベルの解法 (前方差分な ど) を用い，他の 19 回の計算を Leap Frog で行なうべきである．ただし当分の間，線形問題に限って授業 を進めるので，このような操作は不要である．しかし，研究の場面で数値計算をする方程式の多くは非線形 方程式であるので，覚えておくことは必要である．
1.8
離散化された微分の打ち切り誤差
1.1 節の図では，微分を差分で近似するには，真ん中で取る (中央差分) のがいちばん精度がよさそうに見 えるこの節では，きちんと定量的に，離散化に伴う微分の誤差を考えよう． Un+1および Un−1を Unのまわりで，テーラー展開すると，Un+1 = Un+ ∆t(dU/dt)n+ (1/2!)∆t2(d2U/dt2)n+ (1/3!)∆t3(d3U/dt3)n+· · · (1.27) Un−1 = Un− ∆t(dU/dt)n+ (1/2!)∆t2(d2U/dt2)n− (1/3!)∆t3(d3U/dt3)n+· · · (1.28)
ただし，ここで (dU/dt)nは第 n ステップ (t = ∆t× n) における dU/dt を表わす．従って，前方差分，お
よび中央差分は，
(Un+1− Un)/∆t = (dU/dt)n+ (1/2!)∆t(d2U/dt2)n+ (1/3!)∆t2(d3U/dt3)n (1.29)
(Un+1− Un−1)/(2∆t) = (dU/dt)n+ (1/3!)∆t2(d3U/dt3)n (1.30)
となる．いま本当に表わしたいのは右辺第一項 (dU/dt)nだけなのだから，第二項以下は誤差である．数値
計算では，安定性などの要請で，ステップ幅 ∆t は 1 より小さいので，∆t > ∆t2である．結局前方差分は
オーダー ∆t 程度，中央差分は ∆t2程度の誤差を持ち，前方差分の方が 1/∆t だけ誤差が大きいことがわか
第
2
章 拡散方程式
U
t
= aU
xx
の差分解法
2.1
拡散方程式の差分方程式
拡散方程式は， Ut= aUxx (2.1) の形をしている．式の性質については，第一回目の授業ノートを参照されたい． 数値計算を行なうために，拡散方程式 (2.1) を差分化しよう．空間については中央差分とし，時間につい ては前方差分と中央差分（Leap Frog) の 2 通りを考えよう．また以下では常に，Un(j) のように時間ステップを n，空間ステップを j で表す事としよう．空間ステップを i ではなく，j にしたのは虚数単位と区別しや すくするためである．更に，時間ステップ幅を ∆t，空間ステップ幅を ∆x とする． 問題 Ut= aUxxを，空間について中央，時間について前方の差分方程式で書け． 解は，U (n + 1, j) = の形で下に書くこと． (2.2) ヒント:下図のように考えれば，離散化された２階微分表現を求められるだろう． 問題 Ut= aUxxを，空間について中央，時間について中央の差分方程式に書き直せ． (2.3)
2.2
拡散方程式についての
Neumann
の安定性解析
ここでは，拡散方程式について Neumann の安定性解析を行ない．(2.2) および (2.3) が安定に数値計算 可能であるかを調べよう． 拡散方程式でも，後に行なう移流方程式でも，ある方程式を解くときにこちらが任意に決めれるのは，グ リッド幅 ∆t, ∆x である．∆t は小さいほど，∆x は大きいほど安定になる．大抵の場合空間グリッド幅 ∆x は，見たい現象によって定まる．つまり，10km の現象をある程度の精度で計算しようと思えば，少なくと も 2∼3km 程度の ∆x でなくてはならない．従って安定にするためには，ある程度 ∆t を小さくする必要が あるが，むやみに ∆t を小さくすると，計算機代や計算時間がかかりすぎる．よって，∆t を最大どれ位ま でおおきくしていいか？が問題になる． 前方差分についての安定性解析 移流方程式でも出て来るが，空間一次元の時間発展問題に於いては， Un+1(j) = U0,nλ exp (ik∆xj) (2.4) Un(j) = U0,nexp (ik∆xj) (2.5) あるいは，このほうが一般的だが Un+1(j) = U0,nλn+1exp (ik∆xj) (2.6) Un(j) = U0,nλnexp (ik∆xj) (2.7)と置く．どっちでも無論結果は同じである．
まず，(2.6)，(2.7) を (2.2) に代入した結果を下に書け．
(2.8) これを共通因数，U0λneik∆xjで割って得られる λ の一次式を，λ = の形で下に書け．
(2.9) 上式に e±ik∆x= cos (k∆x)± i sin (k∆x) を代入して得られる式を次に書け．
(2.10) ここで，安定であるための ∆t の条件を考えるのだが，その前に k∆x の取り得る範囲を調べよう．まず，
k > 0, ∆x > 0 なので，k∆x > 0 である．k∆x の最大値は，グリッド上で分解できる最小の波長で決まる．
これは，下図のような二波長波 (two grid interval wave）である．この波長は 2∆x，すなわちこの波数は
π/∆x，従って k∆x の最大値は π である．もちろん最小値はゼロ． 安定（λ < 1) のための，∆t の条件を求めよ． (2.11) では Leap Frog ではどうか 以下で空白となっている (2.12)∼(2.15) 式を書け． (2.3) 式を，U0λneik∆xjの表記で書けば， (2.12) 共通因数，U0λn−1eik∆xjで割って，λ で整理すると， (2.13) ここで，前と同じように e±ik∆x= cos (k∆x)± i sin (k∆x) を入れて整理すると，
(2.14) これを d = ∆ta/∆x2と置いて，解の公式に入れると (2.15) この式から，Leap Frog が数値的に安定であるかどうかを論じよ．
2.3
前方差分による，拡散方程式の数値計算
Neumann の安定性解析の結果，安定に計算できるのは，時間に前方・空間に中央で差分化した場合であ ることが理解されたと思う．ここでは，その数値計算を行なおう（EX.4)．次の頁にプログラムの大体の部 分を載せたが，肝心の前方差分を行なう SUBROUTINE である FORWRD の中身は示していない．これを， みんなに作ってもらおう．FORWRD の仕様はであり，サイズ JMAX の実数型配列 DIMNOW（U (j) at t = n× ∆t に相当を）与えられて，次の時
間ステップの値を時間に前方差分を用いて求め，やはりサイズが JMAX である実数型配列 DIMNXT に格 納する． ACOEF は拡散係数 (a)，DX は空間ステップ幅，DT は時間ステップ幅である．また，境界 j=1, JMAX での境界条件は，ux = 0 としよう．これは no-flux condition（境界を通じてやり取りが無い）で
ある．
なお EX.4 において，時間ステップの進行と，MAIN PROGRAM における配列 DIM1・DIM2 との対応 は下図のようになっている． 図の矢印は，n+1 ステップにおける変動の空間分布 DIM1 を求めるためには，n ステップでの配列 DIM2 を用いていることを示している． EX.4 プログラムを完成させ，動かしてみよ． ただし，リスト中の日本語のコメントは打ち込まないこと． プログラムを走らせると，入力要請が出て来るので，まず次のように代入してみよう． 時間ステップ数 1000 初期値の e-folding scale 20 時間ステップ幅 0.01 拡散していく様子が，よく分かるであろう．では時間ステップ幅を 0.1 にしてみよう．今度は爆発してしま うだろう．後で詳しく触れるが，安定性を決めるのは，拡散係数，時間ステップ幅，そして空間ステップ幅 である．今，空間ステップ幅は ∆x = 0.1，拡散係数 = 0.1 に固定しているので，時間ステップ幅のみで安 定・不安定が決まる．いろいろやってみれば，時間ステップ幅 0.05 以内で安定であるという，安定性解析 の結果が正しいことが分かるだろう．なお，空間ステップ幅を固定にしているのは，大抵の場合見たい現象 のスケールを充分解像できるようにという要請で空間ステップ幅は定まり，それに対して安定な時間ステッ プ幅を選ぶ，というステップ幅の決め方になるからである．
第
3
章 移流方程式
u
t
=
−cu
x
3.1
移流方程式の差分式
ここから、移流方程式 ut=−cux (c > 0) (3.1) である。数値計算を始める前に、まずこれを差分方程式で書かなくてはならない。この方程式の差分は、時 間については中央か前方差分、空間については中央か後方差分 (c > 0 の場合、c < 0 なら前方差分) が取り 得る。c > 0 の場合に空間について後方差分なら取り得るが、前方差分は取り得ないのは上記の方程式が一 回目の授業で述べた通り、時間が進むにつれて後方から前方へ情報が伝わっていくという式だからである。 つまり、t = (n + 1)∆t の時点で x = j∆x に存在する情報は、t = n∆t あるいは t = (n− 1)∆t の時点に は、x = j∆x よりも後方にあるのである。この点に関しては、この説明では分かりにくいかも知れないが、 3.5 節での影響領域の議論を読めば、もっとすっきりするだろう。また、この式を空間について後方差分を 取る場合、この差分形を風上差分あるいは上流差分とも呼ぶ。これらの呼称は、情報が風上から風下へ、あ るいは上流から下流へ伝わってく際に、x = j∆x とその風上（上流）側の x = (j− 1)∆x とで差分を求め るという意味を、強く表している。結局、時間・空間に各々2 通りの差分の取り方が可能であるのだから、 計 4 通りの差分方程式が成立する。 問題 (3.1) 式を、Un+1(j) = の形の差分方程式で書き直せ。 時間に前方、空間に後方を (3.2) 式 時間に前方、空間に中央を (3.3) 式 時間に中央、空間に後方を (3.4) 式 時間に中央、空間に中央を (3.5) 式 とし、それぞれ下に書け。 (3.2) (3.3) (3.4) (3.5)3.2
移流方程式の
Neumann
の安定性解析
(3.2)∼(3.5) 式の数値計算を行なう前に、Neumann の安定性解析を行ない、安定な差分式を選ぼう。 問題 (3.2)∼(3.5) 式に対し、それぞれ Neumann の安定性解析を行なえ。 Neumann の安定性を行なうと、c∆t/x というのがよくでてくる。これはクーラント数と呼ばれ、慣例的に µ で書き表わされることが多い。 ヒント 1．(3.2)∼(3.5) の 4 式の内、安定なのは 2 つである。Neumann の安定性解析の結果、安定に数値計算が可能であることが示された差分方程式には、アンダー ラインを引い等の印を付けておこう。 ヒント 2．Neumann の安定性解析は、増幅率 λ の一次式あるいは二次式になる。このうち、二次式でまっ とうには解けない場合があるだろう。 λ2+ aλ + b = 0 ここで、a または b が複素数ならば、 λ = (−b ±√b2− 4ac)/2a のルート内が複素数になってしまい，そう簡単には解けない。こうなってしまった場合、最も不安定 なりやすいのは 2grid interval wave (k∆x = π) であったという拡散方程式で得た経験的知識をもと に、この波長が安定かどうかを解析的に調べ、安定ならばならば数値計算をして、実際に（他の全て の波長についても）安定であるかどうかを確認することとしよう。 ヒント 3．Neumann の安定性解析から、移流方程式の安定な数値計算のためには、µ≤ 1 が条件であるこ とが分かる。これは CFL condition という安定性のための非常に有名な条件である。この条件から時 間ステップ幅は次のように制約される。 µ = c∆t ∆x ≤ 1 → ∆t ≤ ∆x c
3.3
移流方程式の境界条件（radiative boundary condition
と
cyclic
boundary condition)
数値計算を行なう前に、j=1 と jmax での境界条件について考えよう。空間差分は風上差分と中央差分が 可能であった。この二つの差分法で un+1(JMAX) あるいは un+1(1) を求めようとすると何が起るだろうか？
まず、風上差分について考えよう。この差分方程式では、un+1(JMAX) を求めるのに必要なのは、un(JMAX)
と un(JMAX-1) であり、これらの J は共に計算領域内 (1∼JMAX) に入っているので、un+1(JMAX) は問
題無く求めることができる。一方 un+1(1) を求める際には、un(0) という定義されていない値が必要にな る。しかし、計算領域に有意な影響をもたらすエネルギーの供給源が計算領域外に存在することはないと 仮定すると、計算領域の最も上流点では u = 0 とすることができる。すなわち、un+1(1) = 0 であるから、 風上差分の場合には、境界での問題を生じずに、移流方程式の数値計算が可能である。 上記のように風上差分を用いて下の問題で計算を行なうと、初期に与えたでっぱりが、計算領域から抜け ていくのが分かるだろう。このように、境界から波あるいは移流など伝播する性質のものを逃がしてやるの が radiative boundary condition （放射境界条件）である。
空間に中央差分を使用した差分方程式では、un+1(JMAX) を求めるには、un(JMAX-1) と un(JMAX+1)
が必要である。しかし、un(JMAX+1) は計算領域外であり定義されていないので、普通の差分方程式では
計算できない。従って、un(JMAX+1) に相当する値を何らかの方法で定めてやらなくてはならない。ここ
では、cyclic boundary condition を用いよう。これは両端の境界をつないでしまう境界条件である。つま り、un+1(JMAX) を計算するには、un(JMAX+1) の値が必要だが、それを un(1) の値と等しいとするので
ある。同様に un+1(1) を求めるのに必要となる un(0) の値は、un(JMAX) の値と等しいとする。これは、
下図のように一次元空間をくるっと曲げて、両端を 1 グリッド分だけ離してくっつけたものと理解するこ とができる。
3.4
移流方程式の数値計算
EX.5 & 6
(3.2)∼(3.5) のうち、安定に数値計算可能な二つの差分方程式の数値計算を行なえ。EX.4 をコピーし 変更するのが近道であろう。 時間ステップ数は 100∼300、移流速度は 2.0 とする。境界条件は、空間に 風上差分を用いる場合は放射境界条件を、中央差分を用いる場合は cyclic boundary condition を用いる ものとする。
3.5
CFL condition
の物理的解釈
移流方程式は、ある速度での移動を表わしているが、差分方程式で表現し得る最も速い情報の伝搬スピー ドはどのようなものであるかを、考えてみよう。空間・時間に中央の差分方程式をもう一度書く。 U (n + 1, j) = U (n− 1, j) − ∆tc(U(n, j + 1) − U(n, j − 1))/∆x (3.6) これを見ると、(n+1,j) 点は、(n-1,j),(n,j+1),(n,j-1) の点から計算される事が分かる。つまり、下図で○は ◇から、計算されるまた、◇は×から求められる。 同様の操作を繰り返して行くと、下図の斜線部を点○を求めるのに使用された領域が下図斜線部で示さ れる。この領域のことを点○に対する影響領域という。 重要なことは、この傾きの外に何があろうと点○は影響されないということである。斜線部の境界の斜め 線の傾きは、∆t/∆x でこれは遅さ（速度の逆数）の次元をもつ。逆に ∆x/∆t は、この数値計算で表わし 得る最大の速度を示している。物理的な伝播解が ∆x/∆t よりも大きい場合、物理的な解は影響領域の外に あり、数値的に表現できない。従って、伝播が数値的に表現され得るためには、 c ≤ ∆x/∆t (3.7) が必要である。これは、CFL condition である。つまり、CFL condition の物理的意味は、物理解の伝播ス ピードよりも、数値的に表現し得る最大速度のほうが大きくなくてはならないということである。一般に 系に多種の情報の伝播スピードが存在する場合には、最も速い伝播スピードに対し、CFL condition が満た されていなくてはならない。3.6
時間に前方・風上（空間後方）差分の影響領域
時間に前方差分・空間に後方差分を取った場合の影響領域は、下図のようになり、求める点の後方のみか ら情報は伝搬してくる。このことから、移流速度が正である移流方程式についてへ、空間に前方差分を取り 得ないことがわかるであろう。 EX.5 と 6 を動かしたときに気付いた人がいるかもしれないが、空間に中央差分を用いると、上流にも擾 乱が伝搬してゆく。一方、風上差分ではそのような現象は生じない。空間の中央差分における情報の上流へ の伝搬は、再現する波が grid のサイズに近ければ近いほど顕著に現われる。その場合でも、擾乱の重心は 真の伝搬速度で下流に移動する。 これを確認するために、EX.5 と 6 の初期値幅を 1 にして、計算を行なえ。 また、EX5A と 6A に、EX5 と 6 の初期値の計算を下のように変更して矩形を 与えてみよ。 DO 10 J = 1, JMAXを次のように変更する。 DO 10 J = 1, JMAX IF (ABS(J-FJMID).LT.WIDTH/2) THEN DIM1(J) = AMP ELSE DIM1(J) = 0 ENDIF 10 CONTINUE 中央差分による上流への情報の漏れは、Grid を細かく取ることで、たいていの場合その影響を無視する ことができる。しかし、非線形性の強い現象では、それでかならずしもうまくいくとは限らない。場合に よっては、風上差分を使用した方がよいだろう。
3.7
風上差分における人工粘性
EX.5A において風上差分で移流方程式の数値計算を行なって、初期値に矩形波を入れると、伝播にとも なってその形がなまるのが観察されたであろう。このなまりの原因を考えよう。実は EX.5A と 6A のなま りの効果の違いは、空間微分の取り方に起因する。そこで、時間に微分・空間に差分を用いた微分差分方程 式で議論しよう。風上差分における微分・差分方程式は次のように書ける。 ∂Uj/∂t =−c(Uj− Uj−1)/∆x (3.8) ここでは Uj は時間の関数である。Uj−1を Ujのまわりでテーラー展開すると、 ∂U ∂t = − c ∆x{Uj− (UJ− ∆x ∂Uj ∂x + ∆x2 2 ∂2U j ∂x2 − ∆x3 6 ∂3U j ∂x3 )} (3.9) となる。これを整理すると、 ∂Uj ∂t + c ∂Uj ∂x = c∆x 2 ∂2U j ∂x2 − c∆x2 6 ∂3U j ∂x3 (3.10) である。右辺第一項は空間二階微分である。Section 2 を思い出せば、空間二階微分は拡散項であることが分 かる。つまり、本来の移流方程式にはない拡散の効果が、風上差分には入ってしまうのである。これを、人 工粘性あるいは人工拡散という。右辺第二項は第一項よりもはるかに小さいので、まあ気にしないでよい。 問題 空間に中央差分を取るなら、人工粘性は存在しないことを示せ。 この問題から、空間に中央差分を取るなら、人工粘性が効かないことが分かる。また (3.10) 式から人工 粘性は c∆x に比例する、すなわち、グリッド間隔が広いほど大きな人工粘性が生じることが分かる。3.8
分散性 −空間差分の影響−
さて、微分差分方程式を用いて更に空間微分の影響を考えよう。Neumann の安定性解析でも、ある波数 成分に着目する方法（eikxの考え方）が有効であった。つまり、線型方程式については、重ね合わせが成り 立つので、任意の関数はそのフーリエ成分の和で表わされる。各フーリエ成分、すなわち波数成分について 性質が分かれば、その和で表わされる任意関数についても性質が分かったことになる。 空間に中央差分を取った微分差分方程式 ∂Uj ∂t + c Uj+1− Uj−1 2∆x = 0 (3.11)において、ある波数成分の振幅を Ukとし、U
j(t) = Re{Uk(t) exp (ik∆xj)} とすると、 dUk
dt e
ik∆xj+ cUkeik∆x(j+1)− Ukeik∆x(j−1)
2∆x = 0 (3.12)
上式を共通因数 exp (ik∆xj) で割って、exp (iθ) = cos (θ) + i sin (θ) を用いて整理すると、
dUk dt + ikc sin (k∆x) kδx U k= 0 (3.13) となる。一方微分方程式 Ut+ cUx= 0 に U (t) = Re{Uk(t) exp (ikx)} を代入すると、 dUk dt + ikcU k = 0 (3.14) これと (3.13) を比較すると、(3.13) では sin (k∆x)/(k∆) だけ伝播速度が遅くなっていることが分かる。既 に何回か出てきたように、k∆x は 0∼π の範囲を取るので、sin (k∆x)/∆x は 0∼1 の範囲を取る。すなわ ち、二波長波 (k∆x = π) は伝播速度ゼロになり、長波 k∆x→ 0 の極限では sin (k∆x)/(k∆x) → 1 なので 伝播速度は正しく c である。伝播速度が波数によることを、分散性 (dispersion) という。数値計算では、本 来非分散の波も離散化に伴う分散性を持ってしまう。 問題 空間に後方に差分を取った場合の、離散化にともなう分散性を求めよ． 以上の様に、安定に計算できても、正しい結果が得られるとは限らないのが数値計算の恐いところであ る。人工粘性や離散化に伴う分散性などがどうしても入ってしまう。一般にこれらの悪い性質は、グリッド を見たい現象に対して十分細かく取ることで、ある程度抑えられる。しかし、計算機の容量や予算の関係で そう細かく取れない場合は、十分注意する必要がある。 風上差分を用いると情報の逆流を防ぐことができるが、人工粘性がかかってしまう。このように長所があ れば短所もあって、風上差分を用いるか用いないかはどちらの長所・短所に重きを置くかで決まる。
3.9
移流拡散方程式
EX.7 移流拡散方程式 Ut=−V Ux+ ϵUxxを差分化して下に書け。 次に、その数値計算プログラムを作成・実行せよ。 拡散項の計算は、当然時間に前方・空間に中央差分で行なう。また移流項の計算には、 空間に後方差分をとると人工粘性が効き、本来の拡散の効果をみづらくなるので、 空間に中央・時間に中央の Leap Frog とする。また、境界条件は cyclic boundary condition を適用する。パラメータは、DX=1.0、C(移流速度)=0.1、ϵ(拡散係数)=0.01、 与えた初期値が一回りしてはじめの場所に戻ってくるまで計算を行なうこと。作業用配列、WORK(JMAX) を MAIN PROGRAM で宣言すること。そして、 移流項を 時間・空間に中央差分で計算するのと、拡散項を時間に前方・空間に中央で 計算するのは一つの SUBROUTINE で行なう。仮にこの SUBROUTINE を ONESTP と名付けると、その呼び出しは次の様にすること。
CALL ONESTP(DIM2, DIM1, WORK, C, EPS, DX, DT, JMAX)
DIM1,DIM2 はサイズ JMAX の配列である。つまり、WORK を SUBROUTINE の中で 有効に使用することを予想している。なにに使わなくてはならないのか、よく考えて みよう。（ヒントを次の頁に載せた）
時間ステップ幅は、移流方程式・拡散方程式の安定性条件からまずそれぞれ求め、その小さい方の 8 割 り程度の値を用いる。移流拡散方程式では、移流のみ・拡散のみの数値計算よりも、時間ステップ幅をやや 小さく取らなくては、安定性に計算できない。
一般に、いつも Neumann の安定性解析などの手法で、安定な時間ステップ幅を決められるわけではな い。その場合、ある程度の目安を Neumann の安定性解析などでもとめ、それから多少時間ステップ幅を変 えて数値計算を行ない、どの程度の値ならば安定に計算できるかを調べる。ちなみに、移流拡散方程式は、 Neumann の安定性解析が可能なのだが、ここでは割愛する。 ・WORK の使い方のヒント まず、差分方程式を ｕ n+1(j) = の形で表して下に書け。 (3.15) EX6 で行なったように、時間に中央差分を用いて、移流方程式を差分法で計算する場合は、次のように 配列を使用する。すなわち、unは UNOW の配列に格納され、un+1と un−1は UNXT の配列に格納されて
いる。EX6 で un+1と un−1を同じ配列に取って何も問題が無かったのは、un+1(j) を求めるのに、n-1 の時 間ステップにおいては、un−1(j) のみを用いているためである。すなわち、UNXT(j) は un−1(j) として参照 されるやいなや、un+1(j) の値を代入され、その後は、un−1(j) として参照されることはなかったのである。 しかし、(3.15) は un+1(j) を求めるのに、n-1 の時間ステップにおいて un−1(j-1), un−1(j),un−1(j+1) を用 いるようになっているはずである。念のために、移流方程式と移流拡散方程式で un+1(j) を求めるのに必要 な時間・空間グリッドを図示しておこう。 逆に言うと、un−1(j) は un+1(j) と un+1(j+1) の計算に用いられる（un+1(j-1) の計算にも用いられる が、これはここでの議論に本質的ではない）。従って、の差分方程式を計算する時点で（SUBROUTINE の引き渡し時ではなく）un+1と un−1に同じ配列を用いるなら、un+1(j) の計算をする瞬間に UNXT(j) は un−1の値から un+1の値に切り変わってしまう。従って、次のｊの un+1(j+1) の計算に用いられるべき、 un−1(j) を格納しているはずの UNXT(j) の値は un+1(j+1) に置き変わっているので、un+1(j+1) の計算が
正しく行ない得ない。
これは、結局は un+1, un, un−1のそれぞれに別な配列を用意しなくてはならないことを示している。し
かし、MAIN PROGRAM で 3 時間ステップ分の配列と明示的に用意し、それを MAIN の DO 文において 3 通りの場合分け（IF 文）をして用いるのは効率が悪い。SUBROUTINE の中では 3 時間ステップ分の値 を一時的に記憶しなくては計算ができないが、MAIN で記憶するべきなのはあくまで 2 時間ステップ分な のである。すなわち (3.15) の右辺には unと un−1しかない。
そこで、WORK の用い方は例えば次のようになる。SUBROUTINE の中で un−1を WORK に代入し、そ
れ以後では un−1としては配列 WORK を用いる。この方法では、SUBROUTINE 以外では WORK の値を
保持する必要が無い。すなわち、より複雑なプログラムを用い、WORK を必要とする複数の SUBROUTINE があったとしても、MAIN プログラムで宣言する配列は WORK 一つで済むのである。
第
4
章 波動方程式
u
tt
= c
2
u
xx
4.1
波動方程式の解釈
波動方程式は utt= c2uxx (4.1) なる形をした式である。これはスピード c での、x の正負の方向への情報の伝搬をしめす。つまり、(4.1) 式 の右辺を左辺に移行し、次のように変形すると、 (∂ ∂t + c ∂ ∂x)( ∂ ∂t − c ∂ ∂x)u = 0 (4.2) すなわち、この式を満たすためには、+c あるいは−c の速度で伝播する移流方程式のどちらかが、満足さ れなくてはならない。4.2
波動方程式の分解と差分方程式化と数値計算
波動方程式は、2 階の時間微分を含んでいる。この数値計算を行なうには、2 階の時間微分をあらわに扱 う方法と、変数を一個ふやし一階の時間微分の連立方程式を解く方法がある。この授業では後者を選ぼう。 また地球物理関係では、問題とする系が連立微分方程式で与えられている場合も多い。連立方程式をここ で扱っておけば、そのような問題にもなじみやすいであろう。 ある変数 η に対し、波動方程式 (4.1) が次の様に分解できると、仮定しよう。 ut = f (η) (4.3) ηt = g(u) (4.4) ここで、f・g は関数である。これらを、(4.1) utt= c2uxxに代入すると、 utt= (f (η))t= f (ηt) = f (g(u)) (4.5) が得られる。いま f と g はその関数の積が c2∂2/∂x2になればどのようなものでもよいのであるから。簡単 のために、f = g = c∂/∂x としよう。すると、 ut = cηx (4.6) ηt = cux (4.7) が得られる。今 u を速度とすると、η は変位に相当するものである。 (4.6，4.7) の数値計算を行なうために、この式を差分化しなくてはいけない。空間差分を考えると、情報 の伝搬方向が前もって分からないため、風上差分はつかえない。ここで、波動方程式の性質は移流方程式の それとよく似ていることを考えると、安定に数値計算可能な差分の取り方をどうしたらよさそうか？ここま でやってきた皆さんにはピンとくるものがあるだろう。差分の取り方を決定する前に、grid をどのように空 間・時間上に配置するかを見ておこう。ここでは、単に u と η を同じ空間・時間上の点に取ることにしよう。 問 前図のグリッド配置に応じて、(4.6，4.7) の差分方程式を下に書け．(4.8) (4.9) 問 上で出した差分方程式に、Neumann の安定性解析を適用し、やはり CFL condition が 安定のための条件であることを証明せよ。
