2014年4月16日（水）
数値解析 講義
前期水曜３時限
情報理学コース
青柳 睦
aoyagi@cc.九州大学のドメイン
研究室
箱崎キャンパス 情報基盤研究開発センター５階５０２，５００
伊都キャンパス ウエスト2号館 1006室右ドア側
1
講義内容
｛括弧｝の項目はSkipまたは参考程度の項目です
1．非線形方程式の数値解法
1.1 はじめに
1.2 ２分法
→ 演習
1.3 ｛補間法｝
1.4 ニュートン法 → 演習
1.5 多変数問題への応用
2．連立一次方程式の解法
2.1 序論と行列計算の基礎 → 演習
2.2 ガウスの消去法 → 演習
2.3 ３重対角行列の場合の解法
2.4 LU分解法 → 演習
2.5 ｛特異値分解法｝
2.6 共役勾配法 →演習
2.7 反復法 2.7.1 ヤコビ法 →演習
2.8 ガウス・ザイデル法
2.9 ｛ＳＯＲ法（参考）｝
2
3．固有値と固有ベクトルの数値計算
3.1 固有値問題の序論
3.2 ヤコビの方法（対角行列への帰着） →演習
3.3 ギブンスの方法（ ３重対角行列への帰着）→演習？
3.4 DKA法（ ３重対角行列の固有値の計算 )
3.5 ハウスホルダー法（３重対角行列への帰着）
4. その他の数値計算法
4.1 ・・・資料作成中・・
4.2
OFFICE HOUREは 毎週 月曜日午後（箱崎）と水曜日午後（伊都）
です．質問等はこの時間に受け付けます。
3.章4章の内容は多少 変更の可能性があります．
3
成績評価，その他の連絡事項
• 出席点５割＋前期末試験５割
• 講義ノートは，できる限り毎回配りますが
もし欠席したら各自でダウンロードしてください
http://server-500.cc.kyushu-u.ac.jp
（pdf版講義ノートはその週の木曜朝までには公開予定）
• 座席指定でお願いします
約３０分後には出欠を取ります
• 研究室の電話は 092-642-3838と4039
• メールは，aoyagi@cc.九州大学ドメイン
メールSubjectには，必ず「数値解析」と記入
してください
4
1. 非線形方程式の数値解法
1.1 はじめに
方程式
f(x)=0
の解を求める。
f(x)がxに関する１次関数ならば解は簡単に求まる。
f(x)=ax+b
x = -b/a
f(x)がxに関する１次関数以外の場合
例えば、
f(x)=(ax2+b)logx-c
非線形方程式
5
方程式
f(x)=0
の解を求める。
４次以下の代数方程式では、“根の公式”が存在する。
例えば２次方程式では、
ax2+bx+c=0
− b ± b 2 − 4ac
x=
2a
【参考】 例えば、3次方程式では、カルダノの公式、4次方程式では、
フェラリの代数的公式（係数から解を直接求める式）などがある。
【参考】5 次より高次の方程式には、代数的解法 は存在しない（ガロ
ア、ガウスにより証明済）
6
【参考0】 ３次方程式のカルダノ解法（直接法）
x + bx + cx + d =
0
3
2
x = y − b / 3とおき，y 3 + py + q = 0と変形
ただし，p =
c − b 2 / 3, q =
d − bc / 3 + 2b3 / 27
t + qt − p / 27 =
0の解をt1 , t2とし・・
2
3
3
=
u t11/=
, v t21/ 3と書くとき，
解
x1 = u + v − b / 3
x2 =uw + vw − b / 3
2
x3 = uw2 + vw − b / 3, ただしw = exp(2π i / 3)
7
x+logx=0
一般の非線形方程式
(ax2+b)logx-c=0
解を求める公式が存在しない。
近似的解法（コンピュ－タ）を使って強引に解く
逐次近似（ちくじきんじ）
数値計算のあらゆる分野に応用
できる最も基本的で有力な手法
8
逐次近似
逐次法に共通する一般的な考え方
真の値を求めるのではなく、「なるべく真の値に近い値」を求め、
１つの近似値が見つかったら、それを使ってもっと精度の良い
近似値を計算するための“公式”を作る。
xold ← 初期値
まず、１つの近似値（初期値）
を見つけ、それにこの公式を
適用し、その結果をまたこの
公式に入れ・・、これを何度も
繰り返す。
実用上、十分な精度の
近似解が得られたら繰
り返しを終える。
xold ← xnew
xnew ⇐ F ( xold )
“公式”
xnewの精度は
十分か?
Yes
No
1.2
２分法
aとbの中点にcをとる。
「aとbの間に必ず根があって
そして根が，
その個数は一つ」であるとする．
y
解はx軸とy=f(x)の交点だから・・
f(b)
f(c)を計算し、
f(c)
a
中点
f(a)
② ちょうどcのところにある。
③ cとbの間にある。
y=f(x)
c
① aとcの間にある。
x
b
f(a)とf(c)が異符号
①
f(a)とf(c)が同符号
③
f(c)=0
②
根の存在する範囲を狭めていく。（領域削除法）
10
２分法
根の存在する範囲を狭めていき、
範囲がある程度以上小さくなったら、
y
近似値が求まったものとする。
y=f(x)
a
d e
c
第３回
第２回
第１回
b
x
☆ f(x)を計算しているが、
符号しか利用していない。
f(x)の値を評価する計算コス
トが高い場合でも、符号だけ
ならば簡単に求められる場合
は計算量的に有利
11
【参考】 ２分法は１次収束 （またはそれ以下）
f ( x=
) x 2 − 2 初期区域 a=1.0,b=2.0で ４倍精度計算を行った。
n回目の近似値
誤差
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
1.5000000000000000000
1.2500000000000000000
1.3750000000000000000
1.4375000000000000000
1.4062500000000000000
1.4218750000000000000
1.4140625000000000000
1.4179687500000000000
1.4160156250000000000
1.4150390625000000000
1.4145507812500000000
1.4143066406250000000
1.4141845703125000000
1.4142456054687500000
1.4142150878906250000
1.4141998291015625000
1.4142074584960937500
1.4142112731933593750
1.4142131805419921875
1.4142141342163085938
0.0857864376269049512
0.1642135623730950488
0.0392135623730950488
0.0232864376269049512
0.0079635623730950488
0.0076614376269049512
0.0001510623730950488
0.0037551876269049512
0.0018020626269049512
0.0008255001269049512
0.0003372188769049512
0.0000930782519049512
0.0000289920605950488
0.0000320430956549512
0.0000015255175299512
0.0000137332715325488
0.0000061038770012988
0.0000022891797356738
0.0000003818311028613
0.0000005718432135449
12
【参考】
１次収束
繰り返し(n)毎に１桁、2桁、3桁、4桁、・・n桁
と近似値が真値に近づく。
− log10 | approx. − TURE |∝ n
13
２分法が適用できない例
関数f(x)が不連続な場合
方程式が偶数乗根をもつ場合
y=f(x)
a
y=(x-a)2
b
a
14
1.3 補間法(参考)
f(x)の値も利用
PとQを直線で結び、x軸と交わる点をRとする。
R点でのf(x)を計算し、f(x)=0であれば、その点
が根となる。
0でなければ、（Qと符号が異なる端点）SとQを
直線で結び、x軸と交わる点をTとし、これを繰
り返す（次項ニュートン法の原型）。
Q
b−a
xR = b −
f (b )
f (b ) − f (a )
a
R
T
b
P
S
15
1.4 ニュ－トン (またはNewton-Raphson)法
微分係数（勾配）を(も)使う。
y
f ′( xk )
f ′′( xk )
( x − xk ) +
(x − xk )2 +
f ( x ) = f ( xk ) +
1!
2!
接線の方程式
f (x ) − f (xk ) = f ′(xk )(x − xk )
f(x)
X軸との交点を求めると、
xk+1
高次の項は無視
f ( xk )
xk +1 = xk −
f ′( xk )
xk
x
16
Newton-Raphson法の例題
2 をNewton-Raphson法で求める。（初期値x0=2）
f(x)=x2-2
y
2
f ( x0 )
x1 = x0 −
= 2 − = 1.5
4
f ′( x0 )
0.25
f ( x1 )
x2 = x1 −
= 1.5 −
= 1.4147
3
f ′( x1 )
0
x
0.0014
x3 = 1.4147 −
= 1.4142
2.829
たった３回の繰り返しで４ケタの有効桁！
-2
17
☆収束は速い（２次収束）。しかし、初期値を
解の近くにとらなければならない。
☆方程式が偶数乗根をもつ場合にも適用可
例
x
f ( x) = x − cos x
( k +1)
=x
(k )
− cos x
−
, k = 0, 1, ...
(k )
1 + sin x
x
(k )
(k )
18
【参考1】 高次収束の公式
一般式
ニュートン法
２次収束
xn +=
xn −
1
f ( xn )
f ′( xn )
ハレー法
３次収束
f ( xn )
−
xn +=
x
n
1
f ( xn ) f ′′( xn )
′
f ( xn ) −
2 f ′( xn )
４次収束 x = x
n +1
n
f ( xn )
−
f ( xn ){3 f ′′( xn ) f ′( xn ) − f ′′′( xn ) f ( xn )}
f ′( xn ) −
3{2 f ′( xn ) 2 − f ′′( xn ) f ( xn )}
f ( x) = x 2 − a = 0 の近似解
xn 2 − a
xn +=
xn −
1
2 xn
2 xn ( xn 2 − a )
xn +=
xn −
1
3 xn 2 + a
( xn2 − a )(3 xn2 + a )
xn +=
xn −
1
4 xn ( xn2 + a )
19
【参考2】 高次収束の公式を用いた 2 の計算
f ( x=
) x2 − 2
繰り返し
Poorな初期値=10.0で
2次（ニュートン法）
3次（ハレー法）
４倍精度計算をしてみた。
4次の式
1
5.100000000000000000000000
3.509933774834437086092715 2.746078431372549019607843
2
2.746078431372549019607843
1.650475173253007848100583 1.444238094866231938963796
3
1.737194874379598322727878
1.415510038078370553365033 1.414213596802269328471475
4
1.444238094866231938963796
1.414213562645118340447414 1.414213562373095048801689
5
1.414525655148737741224995
1.414213562373095048801689 1.414213562373095048801689
6
1.414213596802269328471475
1.414213562373095048801689 1.414213562373095048801689
7
1.414213562373095467892568
1.414213562373095048801689 1.414213562373095048801689
8
1.414213562373095048801689
1.414213562373095048801689 1.414213562373095048801689
9
1.414213562373095048801689
1.414213562373095048801689 1.414213562373095048801689
− log10 | appx. − true |∝ 2n
繰り返し(n)毎に１桁、2桁、４桁、８桁、１６桁・・
と近似値が真値に近づく。
3n
１桁、３桁、９桁、２７桁・・
と近似値が真値に近づく。
4n
２桁、８桁、３２桁・・
と近似値が真値に近づく。
収束過程で、２分法は誤差が振動するがニュートン法はしない。
20
1.5 多変数の非線形連立方程式への応用
理学工学分野への応用では重要
２変数非線形連立方程式
f ( x, y ) = 0,
g ( x, y ) = 0
の場合は、ニュートン法は
∂f ( k ) ( k )
(x , y )
x ( k +1) x ( k ) ∂x
=
−
y ( k +1) y ( k ) ∂g ( k ) ( k )
(x , y )
∂x
となる。
∂f ( k ) ( k )
(x , y )
∂y
∂g ( k ) ( k )
(x , y )
∂y
−1
f ( x(k ) , y (k ) )
g ( x(k ) , y (k ) )
21
【証明】
∂f
f (=
x, y ) f ( xk , yk ) +
∂x
∂f
( x − xk ) +
∂y
( xk , yk )
( y − yk ) +
( xk , yk )
z = f ( x, y ) の( xk , yk )における接平面・・・①
∂g
g ( x=
, y ) g ( xk , yk ) +
∂x
∂g
( x − xk ) +
∂y
( xk , yk )
( y − yk ) +
( xk , yk )
z = g ( x, y ) の( xk , yk )における接平面・・・②
求めたい近似値は，接平面①②が xy平面( z = 0)と交わる点( xk +1 , yk +1 )
従って， 連立２元１次方程式 f ( xk , y k ) +
g ( xk , y k ) +
∂f
∂f
( x − xk ) +
∂y
∂x ( xk , yk )
∂g
∂x
( x − xk ) +
( xk , y k )
∂g
∂y
( y − yk ) = 0
( xk , y k )
( y − yk ) = 0 ，
( xk , y k )
を ( x, y ) について解けばそれが ( xk +1 , yk +1 ) となる．
22
∂f
∂x
∂f
を ( xk , yk )などと略記し，行列形式で書き直すと・・
∂x
( xk , y k )
∂f
∂x ( xk , yk )
∂g ( xk , yk )
∂x
∂f
( xk , y k )
( x − xk )
f ( xk , y k )
∂y
=
−
g(x , y )
∂g
−
(
)
y
y
k
k
k
( xk , y k )
∂y
よって，
∂f
(x , y )
xk +1 xk ∂x k k
y = y − ∂g
k +1 k ( xk , yk )
∂x
−1
∂f
( xk , y k )
f ( xk , y k )
∂y
∂g
g
x
y
(
,
)
k
k
( xk , y k )
∂y
23
実際には、逆行列の計算は行なわず、
∂f
∂x ( xk , yk )
−
∂g
∂x ( xk , yk )
連立一次方程式
∂f ( k ) ( k )
(x , y )
∂x
∂g ( k ) ( k )
∂x ( x , y )
を解いて
−1
∂f
( xk , yk )
∂y
f ( xk , yk ) h1
=
h とおき・・
g
x
y
(
,
)
∂g
k
k
2
( xk , yk )
∂y
∂f ( k ) ( k )
(x , y )
(k )
(k )
h
(
,
)
f
x
y
∂y
1 = −
∂g ( k ) ( k ) h2 g ( x ( k ) , y ( k ) )
(x , y )
∂y
h1 , h2 を求め、近似解を更新する
x ( k +1) = x ( k ) + h1 ,
y ( k +1) = y ( k ) + h2
を行い、これを繰り返す
24
n変数非線形連立方程式
f1 ( x1 , x2 ,..., xn ) = 0,
f 2 ( x1 , x2 ,..., xn ) = 0,
f n ( x1 , x2 ,..., xn ) = 0
の場合は、簡単のために、記号
x = ( x1 , x2 ,..., xn ),
f = ( f1 , f 2 ,..., f n )
を導入し、上の方程式を
f ( x) = 0
と書く。
25
この場合、ニュートン法は
x
( k +1)
=x
(k )
( k ) −1
− H (x ) f (x )
(k )
と書ける。ただし、
∂f1 ( k )
(x )
∂x1
∂f
(k )
2
x
(
)
(k )
H ( x ) = ∂x1
∂f n ( k )
∂x ( x )
1
∂f1 ( k )
∂f1 ( k )
(x )
(x )
∂x2
∂xn
∂f 2 ( k )
∂f 2 ( k )
(x )
(x )
∂x2
∂xn
∂f n ( k )
∂f n ( k )
(x )
(x )
∂x2
∂xn
ヘシアン行列またはヤコビ行列と呼ばれる．
26
x
( k +1)
(k )
, x ,
f ( x ) は n次元ベクトルで
(k )
であることに注意。
もちろん、実際には、連立一次方程式
H (x ) h = − f (x )
(k )
(k )
を解き、補正
x
( k +1)
=x
(k )
+h
を行い、これを繰り返す。
☆ヤコビ行列が正則であれば
ｈ（繰り返し毎の補正値）は一意に決まる．
27
【参考】ポテンシャル 極値問題
• ポテンシャルエネルギーがV(ｘ，ｙ，ｚ）で与え
られた力学系の平衡点（ｘ，ｙ，ｚ）をニュートン
法で求める．
v = V ( x, y, z ) f ( x, y, z ) =
∂V
∂V
∂V
， g ( x, y, z ) =
，h ( x , y , z ) =
，と書くと・・
∂x
∂y
∂z
安定点，すなわちどの方向への勾配も０となる点を求めるということは，
3変数非線形連立方程式 f ( x, y, z ) = 0，g ( x, y, z ) = 0，h( x, y, z ) = 0 の解
( x, y, z )を求めることと等価である．
ポテンシャル極値問題としての分子の「かたち」
28
Ｖ（ｒ）：２次元の場合の概念図
© Copyright 2024 ExpyDoc
