● 講義資料

● 講義資料

▼ 講義予定

• 数値積分

• Strumの二分法

● 前回の講義のまとめ

▼ 消去法

★ LU 分解

• 現実に連立一次方程式を解く際に,同一のA, 異なるb に対して,Ax=b を何度も解く場合 がある. 特に,あるb1に対してAx=b1 を解き,そのxから得られるb2 に対して解を求め る場合もある. このような場合, Gaussの消去法を何度も適用することは, O(n³)の計算を何 度も行うことになり,非常に効率が悪い.

もし,Aに対して,ある下三角行列L, 上三角行列U を用いて,事前にA=LU と分解され ていれば,Ax=b を解く手順は,LU x=b より,

Ly=b, 前進代入で解く, U x=y, 後退代入で解く

として,ともにO(n²)のアルゴリズムで解を得ることができる. ここで,前進代入とは,後退 代入と同じ手順を下三角行列に適用したものである.

• Aを正則行列としたとき,次の定理が成り立つ.

任意の正則行列Aに対して,Aの行を適当に入れ替えれば,入れ替えたあとの行列を再びA と書いたとき,ある上三角行列U,下三角行列Lが存在して,A=LU と書ける. これをLU 分解と呼ぶ.

この証明はGauss の消去法のアルゴリズムから導かれる. ここでは, Gauss の消去法で枢軸 選択が必要ない（対角成分が 0にはならない）場合のみを考える. この時, ある基本変形行 列の列{Mi}^N_i=1と上三角行列 U が存在して,

MN· · ·M1A=U

と書ける. いま,Mi には行交換の行列が含まれないので,すべてのMi は下三角行列である.

よって,

A= (M₁⁻¹· · ·M_N⁻¹)U

となる. ここで, 下三角行列の逆行列は再び下三角となり,下三角行列の積は下三角なので, M₁⁻¹· · ·M_N⁻¹は下三角行列となる.

• 実際のLU分解は,L= (ℓij),U = (uij)と書いたとき, aij=

j

X

k=1

ℓikukj, i > j,

aii=

i

X

k=1

ℓikuki=

i

X

k=1

ℓijukj, i=j,

aij=

i

X

k=1

ℓikukj, i < j,

(1)

となるが,n²+n個の未知数{ℓij, uij}に対して, n本の方程式があり,このままではL,U を定めることができない. そこで,

– Lの対角成分を1 にとる(Doolittle Type) – U の対角成分を 1にとる(Crout Type) のいずれか一方の条件をおくことが多い.

• いま, Crout Type (uii= 0)を仮定すると, aij=

j−1

X

k=1

ℓikukj+ℓij, i≥j,

aij=

i−1

X

k=1

ℓikukj+ℓiiuij, i < j, となるので,これを書き直すと,

ℓij =aij−

j−1

X

k=1

ℓikukj, i≥j,

uij = 1 ℓii

aij−

i−1

X

k=1

ℓikukj

!

, i < j,

(2)

となる.

• もっとも標準的なアルゴリズムは,Crout Algorithmと呼ばれる次の方法である. （Crout typeを仮定している.）

i= 1, . . . , nに対して,以下を行う.

1. j= 1, . . . , iに対して,

ℓij=aij−

j−1

X

k=1

ℓikukj

を計算する.

2. j=i+ 1, . . . , nに対して,

uij = 1 ℓii

aij−

i−1

X

k=1

ℓikukj

!

を計算する.

• この方法では（多くのLU分解のアルゴリズムでは）, L,U のデータはAの上に上書きで きる. 逆に,Aの上にL,U を上書きした方が,プログラムが書きやすい.

▼ 反復法

★ Jacobiの反復法

• Aは,n×n正方行列で,その対角成分は0 を含まないと仮定する. この時,連立一次方程式 Ax=b を解くために,以下の反復解法を考えることができる.

x⁽⁰⁾ を適当に与え,

x^(k+1)₁ = 1 a11

b1−a12x^(k)₂ − · · · −a1nx^(k)_n , ...

x^(k+1)n = 1 ann

bn−an2x^(k)₂ − · · · −ann−1x^(k)_n−1 , これをJacobiの反復法と呼ぶ.

• Aを対角部分D,（対角部分を含まない）上三角部分U,（対角部分を含まない）下三角部 分Lを使ってA=D+L+U と書いたとき, Jacobiの反復法は

x^(k+1)=AJx^(k)+D⁻¹b, AJ =−D⁻¹(L+U) とかきあらわすことができる.

★ Gauss-Seidel法

• Jacobiの反復法で,x1 から順に計算しているとき,すでに計算済みのx^(k+1) の成分を利用

して反復計算を行った方が,収束が速いと予想できる. そのように改良した以下の反復法を Gauss-Seidel 法と呼ぶ.

x^(k+1)₁ = 1 a11

b1−a12x^(k)₂ − · · · −a1nx^(k)_n , ...

x^(k+1)_ℓ = 1 aℓℓ

bℓ−aℓ1x^(k+1)_ℓ − · · · −aℓℓ−1x^(k+1)_ℓ−1 −aℓℓ+1x^(k)_ℓ+1− · · · −ann−1x^(k)_n−1 , ...

x^(k+1)_n = 1 ann

bn−an2x^(k+1)₂ − · · · −ann−1x^(k+1)_n−1 ,

• Gauss-Seidel法は,

(D+L)x^(k+1)=−U x^(k)+b と書けているので,

x^(k+1)=AGSx^(k)+ (D+L)⁻¹b, AGS=−(D+L)⁻¹U

と書きあらわすことができる.

★ 反復法の収束

• 一般に反復計算は,任意の初期値x⁽⁰⁾ に対して収束するわけではないが,B を n×n行列と し,その絶対値最大の固有値ρ(B)がρ(B)<1をみたすならば,反復計算

x^(k+1)=Bx^(k)+c は,任意の初期値に対して収束し,収束値x^∞ は

x^∞=Bx^∞+c

を満たすことが証明できる. そのときの収束の速さはO(ρ(B)^k)となる.

• 具体的なAに関して,AJ,AGS がρ(AJ)<1,ρ(AGS)<1をみたすことを確かめるのは,一 般には容易ではないが,以下の定理が知られている.

– Aが正定値実対称行列ならば,ρ(AGS)<1が成り立つ.

– Aが対角優位,すなわち, すべての行に対して

|aii|>X

i6=j

|aij| が成り立つならばρ(AJ)<1

▼ 数値積分

• 区間 [0,1] 上の連続関数 f に対する定積分I = Z 1

0

f(x)dx の近似値を計算することを考 える.

• [0,1] 上に相異なる n 個の点 {xi}ⁿ_i=1 と, その点での関数値 {yi}ⁿ_i=1 が与えられたとき, yi=f(xi)をみたすn−1次多項式f がただひとつ定まる. （ラグランジュ補間）

• 代表的な数値積分公式（閉じた公式）

I∼f(1) +f(0)

2 ,

I∼ 1

6(f(0) + 4f(1/2) +f(1)), I∼1

8(f(0) + 3f(1/3) + 3f(2/3) +f(1))

• 次の母関数で定義される数をBernoulli数とよぶ x

e^x−1 =X

n=0

Bn

n!xⁿ. このとき,

Bn(x) =

n

X

j=0

n j

Bjx^n−j

によって定義される多項式{Bn(x)}を Bernoulli多項式と呼ぶ. Bernoulli数, Bernoulli多 項式は次の性質をみたす.

1. B0(x) =B0= 1,B1=−1/2,B2k+1= 0,k≥1.

2. B_k^′(x) =kBk−1(x)

• 台形公式の誤差は,O(h²)である. 一般に2n−1点の公式はO(h²ⁿ), 2n点の公式はO(h²ⁿ) となる. この結果は,オイラー・マクローリンの公式から得ることができる.

• f を区間[0, n]上のなめらかな関数としたとき,

n

X

i=0

f(i) = Z n

0

f(x)dx+1

2(f(0) +f(n)) +

k

X

j=2

Bj

j!(f^(j−1)(n)−f^(j−1)(0)) +Rk, Rk =(−1)^k−1

k!

Z ⁿ

0

B˜k(x)f^(k)(x)dx

が成り立つ. （オイラー・マクローリンの公式）ただし, ˜Bk(x)は Bk(x)を周期1 で拡張し た関数である.

• 区間[0,1]上にn個の点{xi} を取ると,区間[0,1]上の任意のn−1 次以下の多項式f に 対して,

Z 1

0

f(x)dx=Ajf(xj) となる{Aj}を一意的に定めることができる.

• 区間[0,1]上のn次多項式Pnであって,以下の条件をみたす多項式をLegrandre polynomial という.

(Pn, Pm) = Z 1

0

Pn(x)Pm(x)dx=δij

• Legrandre polynomialは,以下の性質を持つ.

1. (n+ 1)Pn+1(x)−(2n+ 1)xPn(x) +nPn−1(x) = 0が成り立つ.

2. Pn は n次多項式である.

3. Pn は,区間[0,1]上に相異なる n個の零点を持つ.

• 区間[0,1]上のn個の点{xi}を,Pn の零点と取ると,区間[0,1]上の任意の2n−1次以下 の多項式f に対して,

Z 1

0

f(x)dx=Ajf(xj) となる{Aj}を一意的に定めることができる. (Gaussの積分法)

▼ Strum の二分法

• 次の条件を満たす関数列{fk}ⁿ_k=0 をStrum 列と呼ぶ.

1. f0(x)は0 でない定数関数である.

2. fk(x)とfk−1(x)は等しい根をもたない.

3. fk(t) = 0 ならばfk+1(t)fk−1(t)<0 となる.

4. fn(t) = 0 ならば,fn^′(t)fn−1(t)<0となる

• 関数列{fk(λ)}ⁿ_k=0 が Strum 列をなすと仮定する. この時, fn(a)fn(b) 6= 0をみたす区間 [a, b]内に存在するfn(x) = 0の実根の数pは,

p=N(b)−N(a) に等しい. ただしN(t)は,

fn(t), fn−1(t), . . . , f0(t)

の符号の変化の数である. ただし, fk−1(t)<0,fk(t) = 0, fk+1(t)>0 となるもので１回の 符号変化と数える. (Strumの定理)

• f(x)は重根を持たない多項式であると仮定する. さらに,あるStrum列{fk(x)}ⁿ_k=0 であっ て,f(x) =fn(x)となるものが存在すると仮定する. また,区間[α, β]内にすべてのf(x)の 根が存在すると仮定する. この時,f(x)の小さい方からk番目の根をλk を求めるための二 分法は以下の通りである.

λk∈[αj, βj],

[αj+1, βj+1]⊂[αj, βj], βj+1−αj+1= βj−αj

2 .

が成り立つ. ただし,αj,βj は,次によって定められる.

µj= αj+βj

2 ,

αj+1=

(αj if N(µj)≥k+ 1, µj if N(µj)< k+ 1 , βj+1=

(µj if N(µj)≥k+ 1, βj if N(µj)< k+ 1 . (Strumの二分法)