数値流体力学

第 3 回 方程式の代数化と連立一次方程式の解法

筑波大学システム情報工学科
構造エネルギー工学専攻
田中聖三

基礎方程式の離散化 => 連立１次方程式 (flow)
Q. 数値解析とは？
A. 基礎方程式を離散化して最終的に得られる連立１次方程式を解く．
基礎方程式（ 1 次元 Laplace 方程式）

差分法で離散化

差分方程式

差分方程式の代数表示

連立１次方程式
計算領域モデル
境界条件

基礎方程式の離散化 => 連立１次方程式 ( 実作業 )差分方程式

境界条件の処理が必要

境界条件の処理が必要

差分方程式の代数表示

基礎方程式の離散化 => 連立１次方程式 ( 実作業 )問題設定

差分方程式の代数表示

境界条件の導入 ( 古典的、正統派なやり方 )連立 1 次方程式の解法
直接法 (Direct Method)
行列の変形や、逆行列に相当するものを計算する。
●Gauss の消去法、 Gauss-Jordan 法、 LU 分解等●
長所：● 安定 ( とりあえず解ける )●疎行列、密行列いずれにも適用可能●
短所：●
反復法よりも記憶容量、計算時間が必要●
桁落ち誤差、丸め誤差が大きくなる可能性●* 数値誤差は？
丸め誤差：
桁落ち誤差：
情報落ち誤差：連立 1 次方程式の解法
反復法 (Iterative Method)
適当な初期解から繰り返し計算により真の解へ収束させる。
●
長所：●
直接法よりもメモリ使用量、計算量が少ない●
並列計算に適している●
短所：● 収束性が問題の影響を受けやすい ( 収束しない場合もある )● 収束性が重要 ( 前処理の導入 )●定常法： 反復計算中に、解ベクトル以外の変数は変化しない。
Jacobi 法、 Gauss-Seidel 法、 SOR 法など。
単純だが、収束性能が悪い
非定常法：拘束、最適化条件が加わる。
CG (Conjugate Gradient: 共役勾配法 )Bi-CG 法 (Bi-Conjugate Gradient: 双共役勾配法 )GMRES(Generalized Minimal RESidual: 一般化残差最小法 )直接法 1 ： Gauss の消去法
連立 1 次方程式を解を変えないように変形し、以下のような形 ( 上三角行列 ) に変形する。
( この変形を前進消去 (Forward Elimination) と言う。 )解は、後退代入 (Backward Substitution) により求められる。

直接法 1 ： Gauss の消去法
例題：

直接法 1 ： Gauss の消去法
プログラム
do i = 1, ndof-1 !Frontward Eliminationai = 1.d0 / a(i,i)
do j = i+1, ndof
cc = a(j,i) * aia(j,i) = 0.0d0
do k = i+1, ndofa(j,k) = a(j,k) - a(i,k) * ccenddox(j) = x(j) - x(i) * cc
enddo
enddo
!
do i = 1, ndof
!Backward Substitution
i1 = ndof + 1 - i
do j = i1+1, ndof
x(i1) = x(i1) - a(i1,j) * x(j)a(i1,j) = 0.0d0
enddo
x(i1) = x(i1) / a(i1,i1)
a(i1,i1) = 0.0d0
enddo直接法 2 ： Gauss-Jordan 法
連立 1 次方程式

を解を変えないように変形し、以下のような形 ( 単位対角行列 ) に変形する。

解は、 b' である。

直接法 2 ： Gauss-Jordan 法
例題：

直接法 2 ： Gauss-Jordan 法
プログラム
do i = 1, ndof!
ai = 1.0d0 / a(i,i)x(i) = x(i) * ai
do j = 1, ndofa(i,j) = a(i,j) * ai
enddo!
do j = 1, ndof
if( i == j ) cycle
cc = a(j,i)
do k = 1, ndofa(j,k) = a(j,k) - cc * a(i,k)enddo
x(j) = x(j) - cc * x(i)
enddo
!
enddo

直接法 3 ：完全 LU 分解
連立 1 次方程式を解を変えないように変形し、次頁のような形 ( 下三角、上三角行列 ) に変形する。

直接法 3 ：完全 LU 分解
簡単のため 4x4 の正方行列とする．となる LU 行列を求める． L ， U はそれぞれ下・上三角行列である．

L, U の成分の計算

直接法：完全ＬＵ分解
例題：

直接法 3 ：完全 LU 分解
L, U の成分の計算

do k = 1, ndofdtmp = 1.d0 / a(k,k)
do i = k+1, ndofa(i,k) = a(i,k) * dtmp
enddo
do j = k+1, ndof
dakj = a(k,j)
do i = k+1, ndofa(i,j) = a(i,j) - a(i,k) * dakj
enddo
enddo
enddo

kk

* a の下半分に L( 対角は省略 ) ，上半分に U が格納されている．

直接法 3 ：完全 LU 分解
解くべき連立 1 次方程式と置くと，となり、について前進代入により解く。

前進代入プログラム：
do i = 1, ndof !Forward substitution for Ly=b
dtmp = 0.d0
do j = 1, i-1
dtmp = dtmp + a(i,j) * x(j)enddox(i) = x(i) - dtmp
enddo

直接法 3 ：完全 LU 分解
解くべき連立 1 次方程式と置くと，が求まれば，となり、となり、について前進代入により解く。

について後退代入により解く。

後退代入プログラム：
do k = 1, ndof !Backward substitution for Ux=yi = ndof - k + 1
dtmp = 0.d0
do j = i+1, ndof
dtmp = dtmp + a(i,j) * x(j)enddo
x(i) = ( x(i) - dtmp ) / a(i,i)
enddo反復法 1 ： Jacobi 法
連立 1 次方程式
L: 下三角、 D: 対角、 U: 上三角

Jacobi 法では、以下のように k 回目の反復解を用いて、 k+1 回目の推定値を求める。

反復法 1 ： Jacobi 法k+1 ステップの推定値

プログラム
x(1:ndof) = 0.0d0
do k = 1, kmax
do i = 1, ndof
dtmp = 0.0d0
do j = 1, ndof
if(i==j) cycle
dtmp = dtmp + a(i,j) * x(j)enddo
xk(i) = (b(i) - dtmp) / a(i,i)
enddox(1:ndof) = xk(1:ndof)
enddo

反復法 2 ： Gauss-Seidel 法
連立 1 次方程式
L: 下三角、 D: 対角、 U: 上三角

Gauss-Jordan 法では、以下のように k 回目の反復解を用いて、ただし、 Jacobi 法と違い、すでに求めた k+1 回目の推定値を使用する。

反復法 2 ： Gauss-Seidel 法
Jacobi 法

Gauss-Seidel 法

プログラム (Jacobi)x(1:ndof) = 0.0d0
do k = 1, kmax
do i = 1, ndof
dtmp = b(i)do j = 1, ndof
if(i==j) cycle
dtmp = dtmp - a(i,j) * x(j)enddo
xk(i) = dtmp / a(i,i)
enddox(1:ndof) = xk(1:ndof)
enddo

プログラム (Gauss-Seidel)x(1:ndof) = 0.0d0
do k = 1, kmax
do i = 1, ndof
dtmp = b(i)do j = 1, ndof
if(i==j) cycle
dtmp = dtmp - a(i,j) * x(j)enddo
x(i) = dtmp / a(i,i)
enddo
enddo反復法 3 ： SOR 法
SOR(Successive Over-Relaxation) 法

プログラム (SOR)
x(1:ndof) = 0.0d0
do k = 1, kmax
do i = 1, ndof
dtmp = b(i)do j = 1, ndof
if(i==j) cycle
dtmp = dtmp - a(i,j) * x(j)enddo
xt = dtmp / a(i,i)x(i) = x(i)+omega*(xt - x(i))
enddo
enddo

Gauss-Seidel 法の修正量に加速パラメータ ω を乗じて修正量を拡大する。 ωは 1 以上の値となるが、大きくしすぎると収束性能が悪くなる。
1.1 ~ 1.3 程度。
問題によっては、 Gauss-Seidel 法(ω=1.0) が良かったりもする。

反復法 4 ： CG 法
CG(Conjugate Gradient, 共役勾配 ) 法を厳密解とするとき、以下の式を最小とする
つまり、下記のを最小とするを求める。

CG 法は任意のから始めて、探索方向 が決まったとすると、を最小とするには：残差を求める：は以下の式により更新できる：の最小値を逐次探索する。

反復法 4 ： CG 法
CG(Conjugate Gradient, 共役勾配 ) 法 (2/3)
探索方向は次の漸化式によって求める：＊探索方向の更新を決める β を決めたい。
ここで、収束した解が求まっているとすると、以下のような直交条件がある：

よって以下が成り立つ：とはとが行列に関して共役反復法 4 ： CG 法
CG(Conjugate Gradient, 共役勾配 ) 法
Initial guess

for k = 0, 1, …, do:ifendexit連立 1 次方程式の解法まとめ
直接法 (Gauss の消去法 ,Gauss-Jordan, LU 分解 )使用制限：
対角項がゼロではない。 (Pivoting による回避が必要 )「定常」反復法 (Jacobi, Gauss-Seidel, SOR)
使用制限：
対角優位 ( 第 i 行の対角項の絶対値がそれ以外の成分の絶対値の和より大きい )「非定常」反復法 (CG)(CG 法の ) 使用制限：
対称正定値行列 (Symmetric Positive Definite)* 非対称行列では多項漸化式を用いる Bi-CG 法
残差を Krylov 部分空間内で最小化する GMRES 法などがある。

課題 ST1( 提出日 :2016 年 11 月 11 日 )1. 上の連立 1 次方程式を Gauss の消去法、 Gauss-Jordan 法、 LU 分解法で解け。
2. 上の連立 1 次方程式を Jacobi 法、 Gauss-Seidel 法、 SOR 法で解け。
3. 本講義で解説した解法の中から一つ選び、プログラミングし動作確認せよ。
4. 下の問題を 3. のプログラムを利用し解け。
基礎方程式（ 1 次元 Laplace 方程式）

レポートは計算過程を詳細に記すこと。 (Gauss の消去法なら、上三角化の過程など )3. のプログラミングの言語はしていしない。コードと動作確認根拠を示すこと。