土 木・建 築 関 連 技 術 者 の た め の
水理計算における差分法入門講座
目 次 １．はじめに ２．計算機の発展史 ３．微分と差分法 ４．現象のモデル化と微分方程式 ５．偏微分方程式の分類と特性 ６．陽解法と陰解法 ７．連立一次方程式の数値解法 ８．楕円型方程式の差分解法 ９．放物型方程式の差分解法 １０．双曲型方程式の差分解法 演習問題（１）
１．はじめに 水理学は、土木技術者にとって必須の知識であることは間違いないが、他の専門 科目、例えば土質力学やコンクリート工学などと比較すると敬遠される嫌いがある。 水理学がこのような扱いを受ける理由の一つに、かつての流体力学的な高等数学を 駆使した説明が多かったことが挙げられる。もちろん、理論の裏付けとしての基礎 知識は必要であるが、そこから引き出されるさまざまな定理を如何に用いるかとい うことにもう少し時間を割いてもいいのではないかと思うこともある。一方、全く 実用的な公式のみを教え、流れの本質を教えないのも片手落ちというものである。 高専や大学を卒業して技術者として活躍されている方々にとっては、実際の設計や 計算を行う際に、基礎的な知識と応用の両者に精通しなければならないのは当然で あるが、時間の関係もあり、どうしても応用的な事柄に重点を置かざるを得ないの が現状ではなかろうか。 本講座では、このように現場で水理計算を行う必要のある、あるいは今後必要に なる可能性のある技術者のために、数値計算、特に差分法の基礎的な知識を紹介し、 今後の参考にしていただきたいと考えている。
２．計算機の発展史 本論に入る前に、数値計算を行うための道具としての計算機の歩みについて振り 返って見たい。そもそも計算機というものがいつ頃出現したのか調べてみると、１ ７世紀のヨーロッパに遡る。イギリスのEdmund Gunter による計算尺やドイツの W.Schickard の四則計算機が始まりだったかも知れないが、よく知られているのは フランスのパスカル（Blaise Pascal、1623-1662）によって発明された世界初のデ ジタルコンピュータであろう。デジタルといっても軸と歯車を使った機械的なもの であった。彼とほぼ同時代には、あのニュートン（Isaac Newton、1642-1727）も 活躍していた。当時ニュートンと大論争したドイツのライプニッツ（Gottfried W. Leibniz、1646-1716）も手回しの計算機を作り天文学の研究に役立てたようである。 １８・１９世紀になるとイギリスのバベッジ（Charles Babbage、1791-1871） が微分解析のための計算機を発明した。記憶装置（メモリ）として1000 ワードを 有しており、このアイデアは IBM 社により後に採用されている。また Jacguard は、この頃レジスターを開発している。 今世紀に入ってはじめの半世紀には、アメリカのハーバード大学 Howard H. Aiken（IBM の創業者）が 1944 年、Mark 1という計算機を開発している。これ は、132 ワードのメモリ内蔵で、掛け算に約 6 秒、足し算に約 3 秒、割り算に約 12 秒を要する計算機であった。また、H. Hollerith が開発した初めてのパンチカー
ドマシンは IBMに採用されている。1943 年、J. Manchly & J. P. Eckert は、アメ
リ カ Aberdeen 研 究 所 に て ENIAC (Electronic Numerical Integrator and
Calculator、電子数値微積分機) を開発した。この計算機には 18000 個の真空管が
使われ、総重量が 30 トンという化け物であった。計算速度としては、加減算に約
200 マイクロ秒を要した。
1950 年代に入り、計算機は急速な発展を遂げる。IBM は 601, 602, 650 などの 機種に機械語を用い、UNIVAC、EDVAC などは FORTRAN を用いた。IBM 709
は32768ワードのメモリで足し算に 24 マイクロ秒の時間を要する計算機であった。
クロ秒の記録を出した。マイクロコンピュータとしてApple II は 128kワードのメ モリを装備していた。1980 年代に入ると集積回路が用いられ、コンピュータ言語 としてALGOL（Algorithmic Language）、BASIC、COBOL、PL/I などが生まれ た。 この計算機の発展によって、計算にかかるコストや時間も変化を遂げてきた。図 ２−１は時間とともに計算コストがどれだけ減少しているかを示した片対数グラ フである。このグラフから言えることは、計算コストを１／１０にするためには約 ８年がかかっているということである。例えば、1950 年頃に$10,000 かかってい た計算は、1990 年頃には$0.10 で済む計算になる。 一方、計算速度については図２−２を参照すれば、1960 年頃に 0.1MFLOPS （FLOPS は１秒間当たりの平均的な実数の四則演算回数、M（メガ）は 106）だ ったのが、1990 年頃には 10GFLOPS（G（ギガ）は 109）に達している。30 年間 に約100,000倍の速さになっている計算となる。おそらく 2000年頃には1TFLOPS （T（テラ）は 1012）の計算速度を実現しているのではないだろうか。例を挙げる と、翼まわりの層流境界層の乱流境界層への遷移計算を行うのに、1950 年頃は人 間が100 年かかる計算だったのが、1950 年頃には IBM650 によって 100 時間で済 むようになり、1970 年頃には IBM370 で計算が 10 分で終わり、1990 年頃には CRAY-YMP を用いればわずか数秒で終わる計算になってしまったということであ る。 さらにこの大型計算機の時代も少しずつ変わり始めている。すなわちパーソナル コンピュータの普及とインターネットが隆盛を極める今日、パソコンでかなりの計 算までやり遂げることが可能になり、インターネットで瞬時に地球の裏側まで情報 交換が行えるようになった。 しかしながら、このような恵まれた状況下でも一つの問題が残っている。それは、 計算を命令する側、つまり人間がどのような指令を下すかということである。単純 な数値シミュレーションでも使う側がそれを理解していなければ唯の時間つぶし と電気の無駄使いにしか過ぎないであろう。したがって、何をどのように計算させ
るかということを認識しておかなければならない。
図２−１ 計算コストの時間的変化
図２−２ 演算速度の高速化
３．微分と差分法 ある物理量が空間内に連続的に分布していて、数学的にいうと微分可能な連続関 数f x( )で与えられているものとする。例を挙げれば、大気の温度が地表から上空に 向けて高さ x の関数として分布している状況を想像することは難しくはないであ ろう。このとき、x のある値x0における温度の微係数f '( )x を求める問題を考える。 微分の定義に従えばこの微係数は 0 0 0 0 ( ) ( ) '( ) lim h f x h f x f x h → + − = (3.1) と表現することができる。計算機においては無限小を扱うことはできないので、あ る有限な大きさの h で代用する。すなわち、 % 0 0 0 ( ) ( ) '( ) f x h f x f x h + − = (3.2) このときの近似の様子は図３−１に示されている。 h が小さければ小さいほど直線 B は直線 A に近づいていく。すなわち、(3.2)式の近似式は、限りなく微係数(3.1) に近づいていくことが分かる。式(3.2)の分子を前進差分、hで割ったものを前進 差分商と呼ぶ。見方を変えると、式(3.2)は、x =x0における微係数をその近傍の関 数値で近似したものと考えることもできる。すると、x=x0における微係数の近似 x y h 0 x x0+h ( ) y= f x 図３−１ 微係数と差分商 A B
方法には他にいろいろなものがあることに気づくであろう。例えば、x= x0−hと 0 x=x における関数値を用いた差分商 % 0 0 0 ( ) ( ) '( ) f x f x h f x h − − = (3.3) や式(3.2)と(3.3)の算術平均から得られる別の差分商 % 0 0 0 ( ) ( ) '( ) 2 f x h f x h f x h + − − = (3.4) などが考えられる。式(3.3)は後退差分商、式(3.4)は中心差分商（中央差分商）と それぞれ名前がつけられている。 以上の３つの差分商は、どれも厳密には両辺が等しくない。ある近似度（誤差、 近似のオーダー）で近似する手段なのである。この誤差を見積もってみよう。例え ば、式(3.2)の誤差を求めるために、関数値f x( 0+h)をx=x0のまわりに Taylor 展 開すると 2 3 0 0 0 0 0 ( ) ( ) '( ) ''( ) '''( ) ... 2! 3! h h f x + =h f x + f x h+f x +f x + (3.5) となる。この式よりf '(x0)を求めると 2 0 0 0 0 0 0 0 ( ) ( ) '( ) ''( ) '''( ) .... 2! 3! ( ) ( ) ( ) f x h f x h h f x f x f x h f x h f x O h h + − = − − + + − = + (3.6) が得られる。ここで、 ( )O h は、それ以降の項が大きくてもせいぜいh の１次のオー ダーであることを表している。したがって式(3.2)の差分商の誤差は ( )O h であるこ とが分かる。同様にして、 f x( 0−h)をx=x0のまわりに Taylor 展開すると 2 3 0 0 0 0 0 ( ) ( ) '( ) ''( ) '''( ) ... 2! 3! h h f x − =h f x −f x h+f x −f x + (3.7) であるから(3.3)の後退差分商の誤差もhの１次のオーダー ( )O h であることが求め られる。一方、式(3.4)の中心差分商については、(3.5)と(3.7)の差を取ることで 求められる。すなわち、
3 0 0 0 0 2 '''( ) ( ) ( ) 2 '( ) ... 3! h f x f x + −h f x −h = hf x + + (3.8) したがって、 2 0 0 0 ( ) ( ) '( ) ( ) 2 f x h f x h f x O h h + − − = + (3.9) が得られ、中心差分商はh の２次のオーダーの誤差をもつことが分かる。 さらに２次の微係数も差分商の差分商で求めることができる。すなわち、 0 0 0 0 0 0 0 0 0 0 2 '( /2) '( /2) ''( ) ( ) ( ) ( ) ( ) ( ) 2 ( ) ( ) f x h f x h f x h f x h f x f x f x h h h h f x h f x f x h h + − − = + − − − − = + − + − = % % (3.10) 式(3.10)の誤差は、簡単な計算によって h の２次のオーダーであることが分かる。 同様にして高次の微係数は差分によって近似することができるが、次数が増えれば 増えるほど、用いる f の数が多くなることも理解されるであろう。
４．現象のモデル化と微分方程式 自然界に起こるいろいろな事象をできるだけ数学的に取り扱うためには何ら かのモデル化が必要である。もし、その事象が何らかの物理的な法則に基づく ものであれば、一般的には独立変数として３次元空間座標と時間座標を考えた 系において、ある物理的な考察をもとにモデル式を組み立てることになる。 ここでは、一つの例として熱伝導の問題を取り上げる。熱が伝わる物質中の 熱の流れをモデル化する際、次のような実験結果を用いることとする。すなわ ち、 (1) 熱は温度の高い方から低い方へ流れる (2) 熱がある面を通して伝わる量はその面積に比例し、また、法線方向の温度 勾配に比例する（比例定数：熱伝導率 k ） (3) 物体の温度が変化するとき、物体内部の熱量の増減は、物体の質量や温度 変化に比例する（比例定数：定積比熱 c ） 以上のような仮定を用いて、図４−１のような３次元空間内の微小な直六面 体要素 V∆ = ∆ × ∆ × ∆x y zの熱の出入りを考えることにする。この要素の中心点の 座標を ( , , )x y z とおく。
z
x
y
Δ
y
Δ
x
Δ
z
(x, y, z)
2 x x f −∆ 2 y y f −∆ 2 x x f +∆ 2 y y f +∆ 図４−１ 微小要素の界面における熱の流れ単位体積当たりの物質の質量（密度）をρとすれば、この要素の質量は m ρ V ρ x y z ∆ = ∆ = ∆ ∆ ∆ (4.1) となる。上記の仮定(3)により、微小要素内で微小時間 t∆ に温度変化 u∆ が生じ たとき、要素内に蓄えられる熱量は H c m u ρc x y z u ∆ = ∆ ∆ = ∆ ∆ ∆ ∆ (4.2) であり、熱が蓄えられる割合は、 H u c x y z t ρ t ∆ = ∆ ∆ ∆ ∆ ∆ ∆ (4.3) で表される。また、もし要素内で何らかの理由で単位時間、単位体積当たり ( , , , ) q x y z t の熱量が発生する場合は、要素内全体で ( , , , ) q x y z t ∆ ∆ ∆x y z (4.4) の熱量が追加されることになる。 さて、この微小要素の各面より出入りする熱量を考えることにする。上記の 仮定(1),(2)より、x 軸に垂直な面 , , 2 x x y z y z −∆ ∆ ∆ を通して単位時間に伝わる熱量は、 , , 2 , , 2 x x y z x x y z u f k y z x ∆ − ∆ − ∂ = − ∆ ∆ ∂ (4.5) 同様にして、面 , , 2 x x y z y z +∆ ∆ ∆ を通して伝わる熱量は、 , , 2 , , 2 x x y z x x y z u f k y z x ∆ + ∆ + ∂ = − ∆ ∆ ∂ (4.6) となる。同様にしてｘ方向、ｙ方向、ｚ方向の全ての熱の流れを列挙し、界面 より要素内に入っていく成分を考えると以下のように記述することができる。 2 2 2 2 2 2 x x y y x x y y z z z z u u u u k y z k y z k x z k x z x x y y u u k x y k x y z z ∆ ∆ ∆ ∆ − + − + ∆ ∆ − + ∂ ∂ ∂ ∂ − ∆ ∆ + ∆ ∆ − ∆ ∆ + ∆ ∆ ∂ ∂ ∂ ∂ ∂ ∂ − ∆ ∆ + ∆ ∆ ∂ ∂ (4.7)
さて、式(4.3)の要素内に熱量が蓄えられる割合は、(4.4)の内部熱発生量と(4.7) の要素境界面より入ってくる熱量の成分の和に等しい。したがって 2 2 2 2 2 2 ( , , , ) x x y y x x y y z z z z u c x y z q x y z t x y z t u u u u k y z k x z x x y y u u k x y z z ρ ∆ ∆ ∆ ∆ + − + − ∆ ∆ + − ∆ ∆ ∆ ∆ = ∆ ∆ ∆ ∆ ∂ ∂ ∂ ∂ + ∆ ∆ − + ∆ ∆ − ∂ ∂ ∂ ∂ ∂ ∂ ∆ ∆ − ∂ ∂ + (4.8) が得られる。両辺を c x y zρ ∆ ∆ ∆ で割り、∆ ∆ ∆ →x, y, z 0ならびに∆ →t 0とすれば、 次の熱伝導方程式が得られる。 2 2 2 2 2 2 ( , , , ) u u u u q x y z t t x y z k κ κ ∂ ∂ ∂ ∂ = + + + ∂ ∂ ∂ ∂ (4.9) ここに、 k c κ ρ = は熱伝導係数と呼ばれる定数である。 以上は、熱伝導方程式を導く一般的な手順であるが、初めに示したある種の 仮定（モデル化）を用いることによって、対象としている物理量を数学的に取 り扱うことが可能になっていることが分かる。水理学を含むその他の力学にお いても、このように我々の身の回りの現象をある連続で微分可能な量の関係と して取り扱い、偏微分方程式などの基礎式を導いてさまざまな解析が行われて いる。
５．偏微分方程式の分類と特性 工学で取り扱う偏微分方程式の多くは二階線型偏微分方程式である。 二階線型偏微分方程式の一般形は、 2 2 2 2 2 ( , ) u ( , ) u ( , ) u ( , ) u ( , ) u ( , ) ( , ) A x y B x y C x y D x y E x y F x y u G x y x x y y x y ∂ + ∂ + ∂ + ∂ + ∂ + = ∂ ∂ ∂ ∂ ∂ ∂ (5.1) であるが、次の基準により３つの型に分類されている。 2 4 0 B − AC> ならば 双曲型 2 4 0 B − AC= ならば 放物型 2 4 0 B − AC< ならば 楕円型 実際に、簡単な二階線型偏微分方程式の例をあげると次のLaplace 方程式（ラプラ ス方程式） 2 2 2 2 0 x y ∂ Φ ∂ Φ+ = ∂ ∂ (5.2) や、Poisson 方程式（ポアソン方程式） 2 2 2 2 q 0 x y ∂ Φ ∂ Φ+ + = ∂ ∂ (5.3) は、B2−4AC=02− × × = − <4 1 1 4 0であるから楕円型の方程式に分類される。 また、次の１次元熱伝導方程式 2 2 u u t κ x ∂ = ∂ ∂ ∂ (5.4) は、 2 2 4 0 4 0 0 B − AC= − × × =κ であるから放物型方程式となる。 一方、１次元波動方程式
2 2 2 2 2 u u c t x ∂ = ∂ ∂ ∂ (5.5) は、判別式より双曲型の微分方程式として分類することができる。 さて、以上３つの型の方程式の特性について考えて見よう。最初の楕円型偏微分 方程式は、従属変数の独立変数に関する２階微分の和が一定、またはゼロとなって いる。物理的には、ある点の従属変数の値が周辺の値より大きくなっているか、小 さくなっているか、等しいかを表している。例えば、ラプラス方程式(5.2)は、ある 点における従属変数の値の分布は周囲と調和して、偏差がなく滑らかにつながって いることを表している。これは主に現象が定常状態に移行したときにある種の平衡 状態として形成されるものである。また、ポアソン方程式においては、その平衡状 態のずれは定数項が原因で起こっていることを示唆している。このように楕円型方 程式の性質はある平衡状態で時間的な変動がない定常状態を表すことが多い。 次に放物型方程式については、独立変数として時間と空間があり、時間に関して は１階微分、空間に関しては２階微分となっている。空間的には周囲との平均化に より徐々にならされていくが、時間に関しては１方向しか自由度がなく、時々刻々 空間的平均化が時間進行で進んでいく現象を表す。熱伝導方程式はこの中に含まれ るが、温度の平均化によって熱は広がっていくことを表している。同様に物質拡散 を表す拡散方程式もこの部類に属している。 最後に双曲型方程式についてはついては、時空間とも２階の微分となっている。 従属変数の空間に関する平均化は時間に関する平均化とある関係で結ばれており、 時間と空間が密接に関係しあっている。これも非定常現象を表す方程式ではあるが、 時間的に平均化が行われるのではなく、時間−空間の間のある決まった関係で解は 波のように振る舞う。この式が波動方程式を表す理由もここにある。
６．陽解法と陰解法 差分解法を大きく分けると陽解法と陰解法に分けることができる。違いを一言で 言えば、離散化した式に１つの未知数しか含まれないのが陽解法であり、一方、式 中に２個以上の未知数が含まれ、式を連立して解かなければならないものが陰解法 となる。例えば、次の偏微分方程式 2 2 u u t κ x ∂ = ∂ ∂ ∂ (6.1) を差分化することを考える。 一つの差分式として次式が考えられえる。 1 1 1 2 2 ( ) n n n n n i i i i i u u u u u t κ x + + − − = − + ∆ ∆ (6.2) 一方、次のような差分式も考えられる。 1 1 1 1 1 1 2 2 ( ) n n n n n i i i i i u u u u u t κ x + + + + + − − = − + ∆ ∆ (6.3) ここで、変数uの上添字n+1 は時刻 t=(n+1)Δt を表し、下添字 i は場所 x=iΔx を 表すものとする。Δt、Δx は時間軸方向及び空間軸方向の刻み幅であり、場合によ っては変化させることも可能である（非一様格子）。簡単のためにここでは刻み幅 は一定の値をもつものと仮定する（一様格子）。 図６−１は、差分式(6.2)の差分形式を表す計算格子の模式図である。黒丸は値が 既知の格子点を表し、白丸は値が未知である格子点を表している。ある１点の未知 の値を求めるのに数点の計算格子における既知の値が使われていることが分かる。 一方、図６−２は、差分式(6.3)に関する同様の図を表している。この場合、式中の 未知数の数が複数個存在するのでこの式だけでは解くことはできない。したがって、 図６−２の場合は、この差分の計算格子を横方向にずらしてそれらを全体として解 くことになる。すなわち、複数の未知数に関する連立方程式を解かなければならな
い。横にずらしていくだけだから、やはり未知数の数が式の数より多くなるのでは ないか、と危惧する人もいるかも知れないが、実は計算領域の端では何らかの境界 条件を与えることになっているので、この連立方程式は解くことができる。 以上の説明の中に出てきた式(6.2)が陽解法の一つで、式(6.3)が陰解法の一つであ るということは容易に理解できるであろう。 ｘ t n n+1 i-1 i i+1 図６−１ 陽解法 ｘ t n n+1 i-1 i i+1 図６−２ 陰解法
７．連立一次方程式の数値解法 差分法として陰解法を採用すれば、連立一次方程式を解く必要があるので、ここ では簡単にその数値解法について述べる。 7.1 直接法（消去法、初等解法） 通常手計算で連立１次方程式を解く際の消去法について手順をまとめたものと 思えば良い。直接法の中でもっとも基本的でポピュラーなのはガウスの消去法であ る。例えば、次の３元連立１次方程式、 11 1 12 2 13 3 1 14 21 1 22 2 23 3 2 24 31 1 32 2 33 3 3 34 ( ) ( ) ( ) a x a x a x b a a x a x a x b a a x a x a x b a + + = = + + = = + + = = (7.1) を解くことを考えてみよう（ただし、a はゼロでないものとする。ゼロの場合は、11 第２式もしくは第３式と入れ替えれば良い）。 まず、第２式と第３式よりx を消去することとする。第１式を1 a で割ると次式 11 が得られる。 13 12 14 1 2 3 11 11 11 a a a x x x a a a + + = (7.2) 第２式、第３式のx を消去するため、式(7.2)を1 a 倍した式と21 a31倍した式、 13 12 14 21 1 21 2 21 3 21 11 11 11 13 12 14 31 1 31 2 31 3 31 11 11 11 a a a a x a x a x a a a a a a a a x a x a x a a a a + + = + + = (7.3) を(7.1)の第２式、第３式よりそれぞれ差し引くこととする。結果として(7.1)式は以 下の連立１次方程式に書き直される。 11 1 12 2 13 3 14 13 12 14 22 21 2 23 21 3 24 21 11 11 11 13 12 14 32 31 2 33 31 3 34 31 11 11 11 ( ) ( ) ( ) ( ) a x a x a x a a a a a a x a a x a a a a a a a a a a x a a x a a a a a + = − + − = − − + − = − (7.4)
次に、式(7.4)の新しい係数を以下のようにおけば、 (1) 12 (1) 13 (1) 14 22 22 21 23 23 21 24 24 21 11 11 11 (1) 12 (1) 13 (1) 14 32 32 31 33 33 31 34 34 31 11 11 11 , , , , a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a = − = − = − = − = − = − (7.5) 連立１次方程式は次式のようになる。すなわち、 11 1 12 2 13 3 14 (1) (1) (1) 22 2 23 3 24 (1) (1) (1) 32 2 33 3 34 a x a x a x a a x a x a a x a x a + + = + = + = (7.6) ただし、 (1) 22 0 a ≠ とする。 次に、(7.6)の第２式、第３式で上と同様の方法を用いて式を変形することにする。 (7.6)の第２式の両辺を (1) 22 a で割れば次式が得られる。 (1) (1) 23 24 2 (1) 3 (1) 22 22 a a x x a a + = (7.7) (7.6)の第３式よりx3を消去するため、(7.7)式に (1) 32 a を乗じて (1) (1) (1) (1) 23 (1) 24 32 2 32 (1) 3 32 (1) 22 22 a a a x a x a a a + = (7.8) (7.6)の第３式より差し引けば、 (1) (1) (1) (1) 23 (1) (1) 24 33 32 (1) 3 34 32 (1) 22 22 (a a a )x a a a a a − = − (7.9) が得られる。ここで、 (1) (1) (2) (1) (1) 23 (2) (1) (1) 24 33 33 32 (1) 34 34 32 (1) 22 22 , a a a a a a a a a a = − = − (7.10) とおけば連立１次方程式は次式のように書き直される。
11 1 12 2 13 3 14 (1) (1) (1) 22 2 23 3 24 (2) (2) 33 3 34 a x a x a x a a x a x a a x a + + = + = = (7.11) 解x x x1, 2, 3は下から次々に代入（後退代入）することで求めることができる。すな わち、 (2) 34 3 (2) 33 (1) (1) 2 (1) 24 23 3 22 1 14 12 2 13 3 11 1 ( ) 1 ( ) a x a x a a x a x a a x a x a = = − = − − (7.12) によって解が得られる。 以上の演算は行列を用いて行うことができる。連立１次方程式(7.1)の係数と定数 項を含む行列 11 12 13 14 21 22 23 24 31 32 33 34 a a a a a a a a a a a a (7.13) を考える。まずa で第１行を割ったものを用いて、第２行、第３行の第１列を消 11 去すると、 11 12 13 14 11 12 13 14 12 13 14 22 21 23 21 24 21 (1) (1) (1) 11 11 11 22 23 24 (1) (1) (1) 32 33 34 12 13 14 32 31 33 31 34 31 11 11 11 0 0 0 0 a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a − − − = − − − (7.14) となる。さらに、得られた行列の２行２列の (1) 22 a で第２行を割ったものを用いて３ 行２列の要素を消去すれば、
11 12 13 14 11 12 13 14 12 13 14 22 21 23 21 24 21 (1) (1) (1) 11 11 11 22 23 24 (2) (2) (1) (1) 33 34 (1) (1) 23 (1) (1) 24 33 32 (1) 34 32 (1) 22 22 0 0 0 0 0 0 a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a − − − = − − (7.15) が得られる。最終的には、式(7.12)のような後退代入によってすべてのx について 解くことができる。 以上の方法は一般的なｎ元連立１次方程式にも適用することができる。ｎ元連立 １次方程式の各係数と定数項で構成される行列 11 12 1 1, 1 21 22 2 2, 1 1 2 , 1 n n n n n n nn n n a a a a a a a a a a a a + + + L L M M O M M L (7.16) を変形して最終的に以下の行列の形になるまで上記の演算を繰り返せば良い。 11 12 13 1 1, 1 (1) (1) (1) (1) 22 23 2 2, 1 (2) (2) (2) 33 3 3, 1 ( 1) ( 1) , 1 0 0 0 0 0 0 n n n n n n n n nn n n a a a a a a a a a a a a a a + + + − − + L L L M M O M M L (7.17) 以上がガウスの消去法の概略である。他にガウス・ジョルダンの消去法、LU 分解 法、コレスキー分解法などがあるが、紙面の関係で割愛することとする。
7.2 反復法 7.1 節の直接法と比較して反復法とはどのようなものか一言でいうならば、連立 １次方程式Ax=bの適当な解の初期値（第０次近似解）x(0)を与え、それを元に第 １近似解、２次、３次….というように近似していき、ある収束半径（計算が正解値 に近づいたかどうかの判定条件）に達するまでこの計算を行い、解を反復して求め る方法である。 具体例として以下に反復法の１解法であるヤコビ法を例示する。ｎ元連立１次方 程式Ax=bの係数行列Aは以下のように分解することができる。すなわち、 1 2 12 13 1 11 21 23 2 22 31 32 1 33 2 1, 1 2 , 1 0 0 0 0 0 0 0 , 0 0 0 0 0 0 0 0 n n n n n n n n nn a a a a a a a a a a a a a a a a − − = + = = A A A A A L L M O O M O M M M M M O L L (7.18) 行列A は行列1 Αの対角成分のみを有する対角行列で、行列Α2は非対角成分のみを 有する行列である。まず、 (0) 1 = A x bを満たす解x(0)を第０次近似解として求めてみ よう。A は対角成分しか持たないので逆行列は 1 11 22 1 1 33 1 0 0 0 1 0 0 1 0 0 0 1 0 0 0 nn a a a a − = A L M M M M O L (7.19) と簡単に求まり、第０次近似解 (0) x は以下のように容易に得られる。 (0) 1 1 2 1 11 22 , , , T n nn b b b a a a − = = x A b L (7.20)
次に第１次近似解x(1)は、上記のx(0)を用いて (1) + (0) = 1 2 A x A x b (7.21) によって求まるものとすれば (1) 1 1 (0) 1 1 − − = − 2 x A b A A x (7.22) 同様にして第ｋ次近似解は、 ( ) 1 1 ( 1) 1 1 k = − − − k− 2 x A b A A x (7.23) となる。このときのx( )k のi番目の成分は一般的に書けば ( ) ( 1) 1 ( ) n k i k i ij ii j ii b x a x a j i a − = = − ≠
∑
(7.24) となる。 以上がヤコビ法の概略である。この他にガウス・ザイデル法、SOR 法などが反 復法として挙げられるが紙面の関係で割愛する。 反復法と前節の直接法はどのような場合に使い分けるかであるが、基本的には元 数が多くなればなるほど反復法による計算の方が計算が少なくて済む。また、疎な 行列（スパース・マトリックス）においても反復法が威力を発揮する。逆に各要素 が０でない密な小さい行列ならば直接法による計算がいい場合もある。 以上の直接法、反復法以外に、特別な行列の連立１次方程式を解く解法がいくつ かある。例えば、１次元の流体解析でしばしば現れる３重対角行列の場合がそうで ある。３重対角行列とは、行列の対角成分とその両側の成分のみ非ゼロ要素でそれ 以外はゼロという行列で、この場合、非常に簡単な漸化式で解が表せる。これも紙 面の関係で割愛するが、以上のいろいろな連立１次方程式の数値解法については既 にいろいろなサブルーチン集がいくつか販売されているので、あらためて自分で作 る必要はないのかも知れない。水理技術者は、これらを使う際にどのような計算を するのか概略を掴めていれば良いと思われる。８．楕円型方程式の差分解法 前述のように楕円型方程式の代表的なものは２次元のLaplace 方程式 2 2 2 2 0 x y ∂ Φ ∂ Φ + = ∂ ∂ (8.1) である。この方程式は、ポテンシャル問題で頻繁に出てくる式で、流れや熱伝導が 定常状態になっているときに現れる。ポテンシャルであるから浸透流問題や地下水 流、物体から離れたところで粘性を無視しうる領域における流れ解析などに適用さ れる。 このLaplace 方程式の境界条件は、領域を囲む外側境界で従属変数の値（ディリ クレ問題）かその境界に垂直な微係数として（ノイマン問題）、あるいはその組み 合わせの条件として与えられる。 代表的な差分形式は次の５点差分公式である。
( )
( )
1, , 1, , 1 , , 1 2 2 2 2 0 i j i j i j i j i j i j x y + − + − Φ − Φ + Φ Φ − Φ + Φ + = ∆ ∆ (8.2) ここで、添字は計算領域を差分格子で切ったときのｘ、ｙ座標を表す。例えば、 , i j Φ は、Φ(
x0+ ∆i x y, 0+ ∆j y)
を表している（図８−１）。 ｘ ｙ i i+1 i+2 1 i− 2 i− j 1 j− 1 j+ 2 j− , i j Φ Φi+1,j 1, i− j Φ , 1 i j− Φ , 1 i j+ Φ 図８−１ 楕円型方程式の５点差分格子また、より精度の高い結果を得るために次の９点差分公式もよく用いられる。
( )
( )
( ) ( )
(
)
( ) ( )
( ) ( )
(
)
1, 1 1, 1 1, 1 1, 1 2 2 2 2 1, 1, , 1 , 1 2 2 2 2 , 5 5 2 2 20 0 i j i j i j i j i j i j i j i j i j x y x y x y x y + + − + + − − − + − + − Φ + Φ + Φ +Φ ∆ − ∆ ∆ − ∆ − Φ + Φ + Φ +Φ ∆ + ∆ ∆ + ∆ − Φ = (8.3) 打ち切り誤差（精度）は、５点公式が 2 2 ( , ) O ∆x ∆y 程度で、９点公式はO(∆x4,∆y4) 程度である。 実際の計算では、計算領域を∆ ∆x, yの等間隔に区切り、領域内部に従属変数 , ( 0,1, ; 0,1, , ) i j i ie j je Φ = L = L を設定し、個々の格子点について上記の差分式を用い る。結果として出てくるのは、Φi j, (i=0,1,Lie j; =0,1,L,je)に関する連立１次方程 式である。計算境界においては、Φもしくは n ∂Φ ∂ （n は境界に垂直方向を意味する） が与えられるので、連立方程式の一番外側の値については何らかの情報が与えられ る。したがって、未知数の数と方程式の数が一致し、この連立１次方程式を直接法 あるいは反復法で解けば良い。９．放物型方程式の差分解法 放物型の代表的な偏微分方程式は、１次元の熱伝導方程式（もしくは拡散方程式） である。 2 2 u u t κ x ∂ ∂ = ∂ ∂ (9.1) この式は、１次元等方性媒質中の熱伝導や拡散を支配する方程式であり、また、放 物型の境界層方程式の簡単なモデル式としても用いられている。 以下に幾つかの重要な放物型方程式の差分スキームについて概説する。 9.1 簡単な陽的解法 次式は陽的な１ステップ形式の差分解法である。
( )
1 1 1 2 2 n n n n n i i i i i u u u u u t κ x + + − − = − + ∆ ∆ (9.2) 打ち切り誤差は 2 ( , ) O ∆ ∆t x である。つまり時間に関して１次精度、空間に関して２ 次精度の差分スキームである。 9.2 簡単な陰的解法（Laasonen Method） 熱伝導方程式の簡単な陰的解法は、Laasonen(1949)によって提案された。その 差分形式は、次式のようである。( )
1 1 1 1 1 1 2 2 n n n n n i i i i i u u u u u t κ x + + + + + − − = − + ∆ ∆ (9.3) 打ち切り誤差は 2 ( , ) O ∆ ∆t x である。後の章で差分解法の数値的安定性についても論 じるが、一般的には陰解法は無条件安定（ t∆ を大きくとれる）であることが多く、 このLaasonen スキームもその１つである。9.3 Richardson Method 放物型方程式に対して時間と空間に２次精度を期待して提案された次式の Richardson Method(1910)は、全ての∆ ∆t, xに対して不安定なスキームであること が知られている。 すなわち、
( )
1 1 1 1 2 2 2 n n n n n i i i i i u u u u u t κ x + − + − − = − + ∆ ∆ (9.4)Forsythe & Wasow(1960)は、「この式は、いまでは形式的に考えられる計算過程に 関する数学的解析の必要性を想起させるものとしてのただ歴史的関心をひくもの となっている。」と述べている。つまり、この差分解法は不安定で使い物にはなら ないが、数学的な安定解析の重要性を考えさせるものだということである。 9.4 DuFort-Frankel Method（蛙跳び法、Leapfrog 法） 上述の不安定なRichardson Method における n i u を時間平均表示 1 1 ( n n ) / 2 i i u + +u − で 置き換えることでスキームを無条件安定とすることができる。最終的な式形は、
( )
1 1 1 1 1 1 2 2 n n n n n n i i i i i i u u u u u u t κ x + − + − + − − = − − + ∆ ∆ (9.5)となる。この差分スキームを初めて提案したのはDuFort & Frankel(1953)である。
彼らはこの差分解法を蛙跳び法（Leapfrog 法）およびピラミッド法と呼んだ（図 ９−１参照）。打ち切り誤差は 2 2 ( , ) O ∆ ∆t x であり、この差分は陽解法であるにも関 わらず無条件安定な差分スキームである。 9.5 Crank-Nicolson Method
Crank & Nicolson(1947)は、熱伝導方程式を解くために次式のような陰的差分解 法を用いた。
( )
1 1 1 1 1 1 1 1 2 ( 2 ) ( 2 ) 2 n n n n n n n n i i i i i i i i u u u u u u u u t κ x + + + + + − + − − − + + − + = ∆ ∆ (9.6)このアルゴリズムは無条件安定な陰的解法であり、クランク・ニコルソンスキー ム（Crank-Nicolson Scheme）として非常に良く知られている。このスキームの打 ち切り誤差は 2 2 ( , ) O ∆ ∆t x である。差分形式の構造を図９−２に示す。 以上が１次元放物型方程式に関する差分解法の概略である。２次元の同方程式に ついてはまた特殊な手法が存在するが、それについては後述することとする。 ｘ ｙ i 1 i− i+1 1 n− n 1 n+ 1 n i u − 1 n i u+ 1 n i u− 1 n i u + 図９−１ DuFort-Frankel Method ｘ ｙ i 1 i− i+1 n 1 n+ n i u 1 n i u+ 1 n i u− 1 n i u + 図９−２ Crank-Nicolson Method 1 1 n i u++ 1 1 n i u−+
１０．双曲型方程式の差分解法 双曲型方程式の典型的な例は、次式の１次元波動方程式である。 2 2 2 2 2 u u c t x ∂ = ∂ ∂ ∂ (10.1) この方程式は、一様な媒質中を波速 c で進む音波の伝播を支配する方程式である。 この式と性質が似ている１階の偏微分方程式は、 0 ( 0) u u c c t x ∂ + ∂ = > ∂ ∂ (10.2) 式(10.1)は(10.2)より求めることができる。式(10.2)は１階の１次元波動方程式であ るが、ここではこの式を双曲型方程式のモデル方程式として採用することとする。 この式は、ｘ方向に c のスピードで伝わる波を記述する線型双曲型方程式であり、 また、非粘性流れを支配する非線型方程式のモデル方程式としても利用されている。 注意すべきことは、ここでは式(10.2)を波動方程式として今後その差分解法を概説 していくが、本来の古典的な波動方程式は式(10.1)を指すものと考えてほしい。式 (10.2)は、しばしば、対流方程式（移流方程式）とも呼ばれる。 10.1 Euler 陽的差分法 次に示す簡単で陽的な１ステップの計算法 1 1 0 n n n n i i i i u u u u c t x + + − + − = ∆ ∆ (10.3) 1 1 1 0 2 n n n n i i i i u u u u c t x + + − − + − = ∆ ∆ (10.4) は、それぞれ打ち切り誤差がO(∆ ∆t, x)と 2 ( , ) O ∆ ∆t x であるが、安定解析を行うと残 念なことにどちらのスキームも無条件不安定となる。(10.3)式は風下差分、(10.4) 式は中央差分という名前がついているが、どちらも波動方程式の計算に用いること はできない。
10.2 風上差分（Upwind Difference Method） 上述のEuler 陽的差分において、もし波速c が正ならば、空間方向の差分に前進 差分の代わりに後退差分を用いることで安定な計算スキームが得られる。また、も し、波速cが負であれば、安定性を確保するために空間方向に前進差分を用いなけ ればならない。波速が正の場合の差分式は、 1 1 0 n n n n i i i i u u u u c t x + − − + − = ∆ ∆ (10.5) この差分スキームは打ち切り誤差が ( , )O ∆ ∆t x の差分解法であり、条件付き安定であ る。 10.3 Lax Method Euler 陽的差分法で空間方向に中心差分を用いた不安定な差分式(10.4)式におい て、 n i u を平均値( n1 n1) / 2 i i u+ +u− で置き換えるとスキームが条件付き安定となる。この 差分解法はLax Method(1954)と呼ばれ、以下のような差分形式となっている。 1 1 1 1 1 1 ( ) 2 0 2 n n n n n i i i i i u u u u u c t x + + − + − − + − + = ∆ ∆ (10.6) この陽的な１ステップの差分解法は１次精度を有している。 10.4 Euler 陰的差分法 いままでの差分解法は陽的解法のみであった。以下に陰的解法について述べるこ ととする。次式の陰的スキーム 1 1 1 1 1 0 2 n n n n i i i i u u u u c t x + + + + − − − + = ∆ ∆ (10.7) は、１次精度で打ち切り誤差は 2 ( , ) O ∆ ∆t x である。安定解析によれば無条件安定で あるが、陰的差分法のため毎ステップごとに連立１次方程式を解く必要がある。
10.5 Leap Frog Method
Leap Frog Method は、放物型方程式の差分解法においては無条件安定な陽的差 分解法として紹介されたが、双曲型方程式においては条件付き安定となる。また、
これまで波動方程式の１次精度差分法のみ扱ってきたが、このLeap Frog Method
は簡単な２次精度差分解法である。差分式形は、 1 1 1 1 0 2 2 n n n n i i i i u u u u c t x + − + − − + − = ∆ ∆ (10.8) となる。時間方向に３つのレベルの格子点を扱い、そのうち過去の２時間ステップ の値が既知である必要がある。プログラミングを行う際にはこの点に注意が必要で ある。 10.6 Law-Wendroff Method
Lax-Wendroff Method は、次のような Taylor級数展開を行うことで得られる。
( )
2 2( )
3 1 2 1 2 n n i i u u u u t t O t t t + = + ∆ ∂ + ∆ ∂ + ∆ ∂ ∂ (10.9) 波動方程式 2 2 2 2 2 , u u u u c c t x t x ∂ = − ∂ ∂ = ∂ ∂ ∂ ∂ ∂ (10.10) を用いて(10.9)式は次式のように書き換えることができる。( )
2 2( )
3 1 2 2 1 2 n n i i u u u u c t c t O t x x + = − ∆ ∂ + ∆ ∂ + ∆ ∂ ∂ (10.11) 最後に u x ∂ ∂ と 2 2 u x ∂ ∂ を２次精度の中央差分で置き換えれば良く知られたLax-Wendroff スキームが得られる。(
)
2( )
( )
2(
)
1 1 1 2 1 2 1 2 2 n n n n n n n i i i i i i i c t c t u u u u u u u x x + + − + − ∆ ∆ = − − + − + ∆ ∆ (10.12) この陽的な１ステップの差分スキームは２次精度の差分解法である。10.7 ２段階 Lax-Wendroff Method 非粘性流れの方程式などの非線型方程式に対しては元のLax-Wendroff スキーム を変形した２段階解法が用いられる。波動方程式に適用すると、 第１ステップ 1 / 2 1 / 2 1 1 1 ( ) 2 0 / 2 n n n n n i i i i i u u u u u c t x + + + + − + − + = ∆ ∆ (10.13) 第２ステップ 1 1 / 2 1 / 2 1 / 2 1 / 2 0 n n n n i i i i u u u u c t x + + + + − − − + = ∆ ∆ (10.14) となる。 このスキームは２次精度の差分解法である。第１ステップにおいて中間の時間ス テップの格子中間点i+1 / 2における値を求め、第２ステップにおいて残りの半分の
時間ステップに対して Leap Frog Method を用いている。もし、この２段階
Lax-Wendroff Method を線型方程式に適用すると、元の Lax-Wendroff Method と 一致する。したがって、この２段階解法は、特に非線形方程式を陽的に解くための 計算法と思えば良い。
10.8 MacCormack Method
MacCormack Method(1969)は、流れの方程式を解くために広く用いられている
手法である。このMacCormack Method は、２段階 Lax-Wendroff Method を修正
して、格子中間点i+1 / 2とi−1 / 2における未知数を計算しなくてよくしたものであ る。この特徴のためにこのスキームは、非線型偏微分方程式を解く際に特に有用で ある。線型の波動方程式へ適用すると、この陽的な予測子、修正子法は次式のよう になる。 予測子 1 1 ( ) n n n n i i i i t u u c u u x + + ∆ = − − ∆ (10.15)
修正子
(
)
1 1 1 1 1 1 2 n n n n n i i i i i t u u u c u u x + + + + − ∆ = + − − ∆ (10.16) 1 n i u + は時間ステップ n+1 における仮の予測値である。修正子において n+1 ステッ プの最終的なu の値が求められる。予測子では空間差分として前進差分が用いられ、 修正子では逆に後退差分が用いられていることに注意しなければならない。この差 分方向の逆転により誤差の蓄積を防ぐことができる。この MacCormack Method に線型の波動方程式を適用すれば、Lax-Wendroff スキームに一致することになる。 この他、さまざまな差分解法が双曲型方程式においても存在するが、紙面の関係 で割愛する。また、差分により全領域格子点を計算する方法以外に、双曲型方程式 の数値計算においては特性曲線法を用いた手法があるが、数値的解法に適用するこ とが困難なことが多いのでこれも省略する。参考文献 （１）小門純一、八田夏夫：「数値計算法の基礎と応用」、森北出版、1988. （２）C.R.ワイリー：「工業数学 上（富久泰明訳）」、ブレイン図書出版、1962. （３）フォーサイス、ワソー：「偏微分方程式の差分法による近似解法（上） （藤野精一訳）」、吉岡書店、1968. （４）G.D.スミス：「コンピュータによる偏微分方程式の解法（藤川洋一郎訳）」、 サイエンス社、1971. （５）戸川隼人：「数値解析とシミュレーション」、共立出版、1976. （６）日本機械工学会編：「流れの数値シミュレーション」、コロナ社、1988. （７）土木学会編：「土木工学における数値解析／基礎編」、サイエンス社、1974. （８）高橋亮一、棚町芳弘：「差分法」、培風館、1991. （９）スタンリー、ファーロウ：「偏微分方程式 科学者・技術者のための 使い方と解き方（伊理正夫、伊理由美訳）」、啓学出版、1983.
（10）Ching-Jen Chen：Computational Methods in Fluids and Heat Transfer,
Lecture Note, Department of Mechanical Engineering, The University of Iowa, 1991
演習問題（１） 【１】式(3.10)の２次の微係数の差分商の誤差が 2 ( ) O h であることを示しなさい。 【２】物理的な考察のもとに２次元の熱伝導方程式を誘導しなさい。 【３】１次元の計算領域において熱伝導方程式 2 2 u u t κ x ∂ = ∂ ∂ ∂ を数値的に解くことを考える。左端の温度がu0 =0、右端の温度がu10 =1.0に保たれ ているような境界条件のもとで、途中の等間隔に刻んだ格子点における温度 ( 1,2,....,9) i u i= はそれぞれどのような代数方程式の関係になるか求めよ。ただし、 差分解法には式（6.2）を用いるものとする。 【４】【３】の問題において差分解法として式(6.3)を採用することにすれば、代数 方程式はどのようになるか求めなさい。 【５】線型の１次元１階波動方程式 0 u u c t x ∂ + ∂ = ∂ ∂
に対してMacCormack Method を適用した結果が Lax-Wendroff Method に一致す
