2005年4月１3日（水）2005年度
数値解析講義 前期水曜３限
青柳 睦
aoyagi@cc.kyushu-u.ac.jp
研究室：情報基盤センター５階５００，５０２1

講義内容｛括弧｝の項目はSkipまたは参考程度

１．超越方程式（非線形方程式）の解法
１．１ ２分法 → 演習
１．２ ｛補間法｝
１．３ ニュートン法 → 演習

２．代数方程式の解法
２．１ グラーフェの方法 → （演習）２．２ ｛ベルヌーイの方法｝

３．連立一次方程式の解法
３．１ 行列計算 → 演習
３．２ ガウスの消去法 → 演習
３．３ ３重対角行列の場合の解法
３．４ LU分解法 → 演習

講義ｖｓ演習スケジュールに余裕がある日を見つけて，「数値解析」を必要とする理学工学(社会学？)分野で活用されている計算科学の2概論を紹介します（予定）．

３．５ 特異値分解法
３．６ 共役勾配法
３．７ 反復法

４．固有値と固有ベクトルの近似計算４．１ レバリエールの方法
４．２ ３重対角行列の固有値と固有ベクトルの計算
４．３ ランチョスの方法（３重対角行列への帰着）
４．４ ギブンスの方法（対称な３重対角行列への帰着）
４．５ ヤコビの方法（対角行列への帰着）

５．ラグランジュの未定乗数法
６．勾配法
６．１ 最急降下法
６．２ 最急降下加速法
６．３ ステップ幅の決め方

4章 「固有値・・」以降は，講義と演習の進捗状況によっては項目を絞った内容にする可能性があります・・．
3成績評価，その他• 出席点５割＋前期末試験５割• 講義ノートは，毎回配りますがもし欠席したら各自でDLしてください
server-500.cc.kyushu-u.ac.jp（講義ノートはその週の土曜日深夜までには公開予定）• 次回から座席指定でお願いします
３０分後には座席表から出欠を取ります• 研究室の電話は 092-642-3838• メールは，aoyagi@cc.kyushu-u.ac.jpSubjectには，「数値解析」と記入してください4非線形方程式
非線形方程式方程式

f(x)=0の解を求める。

f(x)がxに関する１次関数ならば解は簡単に求まる。

f(x)=ax+bx = -b/af(x)がxに関する１次関数以外の場合

非線形方程式f(x)=(ax2+b)logx-c5

非線形方程式
非線形方程式方程式

f(x)=0の解を求める。
４次以下の代数方程式では、 根の公式 が存在する。
例えば２次方程式では、

ax2+bx+c=0− b ± b 2 − 4acx=2a

例えば４次方程式では、フェラリの公式（Ferrari’sformula)6

x+logx=0
一般の非線形方程式(ax2+b)logx-c=0解を求める公式が存在しない。

近似的解法（コンピュ−タ）を使って強引に解く

逐次近似（ちくじきんじ）

数値計算のあらゆる分野に応用できる最も基本的で有力な手法

7逐次近似逐次近似

逐次法に共通する一般的な考え方真の値を求めるのではなく、「なるべく真の値に近い値」を求める。
（精度の良い近似値を求める。）「１つの近似値が見つかったら、それを使ってもっと精度の良い近似値を計算する公式を」作る。
まず、１つの近似値を見つけ、それにこの公式を適用し、その結果をまたこの公式に入れ、これを何度も繰り返す。
最終的に実用上十分な精度の解が得られる。
計算科学的な観点からは・・「逐次法」は 問題サイズ(n) に対する演算量の正確な予測が得にくい，一方，「直接法」は演算量が明確に決まる場合が多い・・
e.g. O(n2), O(log n) etc ．

8

５次方程式を解く
５次方程式を解く
５次方程式

x5+0.1=x4+x2+xグラフを描いて解くを解く。

まず、

f(x)=x5-x4-x2-x+0.1の値をx=0.0から0.01おきに2.0まで表を作って計算し、グラフにプロットする。
9f(x)=x5-x4-x2-x+0.1#include<stdio.h>#include<math.h>main(){
float x;float y;x=0.0から0.01おきに2.0まで計算

yの値がゼロに近づいたら、xの値を出力

for(x=0.0; x<=2.0; x=x+0.01){
y=pow(x,5)-pow(x,4)-pow(x,2)-x+0.1;
if((y<0.01)&&(y>-0.01)){
printf("x=%5.2f¥ty=%5.3f¥n",x,y);
}
}
}

10

グラフを描く
グラフを描く
x0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
1.6
1.7
1.8
1.9
2f(x)
0.100
-0.010-0.141-0.296-0.475-0.681-0.912-1.162-1.422-1.676-1.900-2.064-2.125-2.033-1.723-1.119-0.128
1.356
3.458
6.319
10.100

原始的で強引なやり方

f(x)=x5-x4-x2-x+0.1
12
108解

6
420-2-40

0.5

1

1.5

211

２分法
２分法

「aとbの間に必ず根があってその個数は一つ」

aとbの中点にcをとる。 そして根が① aとcの間にある。

y② ちょうどcのところにある。
f(b)③ cとbの間にある。

y=f(x)

f(c)を計算し、

f(c)

ac
f(a)

xbf(a)とf(c)が異符号①f(a)とf(c)が同符号③f(c)=0②根の存在する範囲を狭めていく。（領域削除法） 12

２分法
２分法

根の存在する範囲を狭めていき、範囲がある程度以上小さくなったら、y近似値が求まったものとする。
y=f(x)

ad ec

第３回
第２回

bxf(x)を計算しているが、符号しか利用していない。

第１回
13

２分法が適用できない例
２分法が適用できない例¾関数f(x)が不連続な場合¾方程式が偶数乗根をもつ場合
y=f(x)

ay=(x-a)2ba14

補間法
補間法

f(x)の絶対値も利用

PとQを直線で結び、x軸と交わる点をRとする。
R点でのf(x)を計算し、f(x)=0であれば、その点が根となる。
0でなければ、SとQを直線で結び、x軸と交わる点をTとし、これを繰り返す（次項ニュートン法の原型）。

Qb−af (b )xR = b −f (b ) − f (a )
aR

TbP

S

15

ニュ−トンニュ−トン (またはNewton-Raphson)法(またはNewton-Raphson)法
微係数（勾配）を使う。

f ′( xk )f ′′( xk )( x − xk ) +(x − xk )2 + K
f ( x ) = f ( xk ) +
1!
2!

yf(x)

高次の項は無視f ( xk )xk +1 = xk −f ′( xk )xk+1xkx16

∂f∂f
f ( x, y ) = f ( xk , yk ) +( x − xk ) +( y − yk ) + K∂x ( xk , yk )∂y ( xk , yk )g ( x, y ) = g ( xk , yk ) +∂g∂gx−x+( k)( y − yk ) + K∂x ( xk , yk )∂y ( xk , yk )17

Newton-Raphson法の例題
Newton-Raphson法の例題
2 をNewton-Raphson法で求める。（初期値x0=2）

f(x)=x2-2

yx1 = x0 −f ( x0 )2= 2 − = 1.5
f ′( x0 )4f ( x1 )0.25
x2 = x1 −= 1.5 −
= 1.4147f ′( x1 )30xx3 = 1.4147 −-20.0014= 1.4142
2.829
18

ニュートン法の復習（参考１）
１変数非線形方程式f ( x) = 0を解く方法。
ニュートン法の原理

f (x) の x = x のまわりのテイラー展開
f ( x) = f ( x ) + f ' ( x )( x − x ) + L

より

0 = f ( x) = f ( x ) + f ' ( x )h19

これを解くと
f (x)
h=−
: 補正量
f '(x)f (x)~x = x+h= x−
f '(x)
反復式
(k )
f(x)
( k +1)(k )x=x −, k = 0, 1, ...
(k )f '(x )x ( 0 ) ： 初期値20

☆収束は速い（２次収束）。しかし、初期値を解の近くにとらなければならない。
¾方程式が偶数乗根をもつ場合にも適用可例

f ( x) = x − cos x

(k )
(k )
x−cosxx ( k +1) = x ( k ) −, k = 0, 1, ...
(k )1 + sin x

21

（参考２） 多変数への応用
理学工学の分野では重要！

２変数非線形連立方程式

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となる。

−1∂f (k ) (k ) (x , y )  f ( x( k ) , y ( k ) ) ∂y ∂g (k ) (k )   g( x(k ) , y(k ) ) (x , y )  ∂y22

実際には、逆行列の計算は行わなわず、
連立一次方程式 ∂f (k ) (k )
 (x , y ) ∂x ∂g (k ) (k ) ∂x ( x , y )∂f (k ) (k ) (x , y ) (k )(k )hfxy(,)
∂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を行い、これを繰り返す。
23

n変数非線形連立方程式
f 1 ( x1 , x 2 ,..., x n ) = 0 ,f 2 ( x1 , x 2 ,..., x n ) = 0 ,M

f n ( x1 , x 2 ,..., x n ) = 0の場合は、簡単のために、記号x = ( x1 , x 2 ,..., x n ),f = ( f 1 , f 2 ,..., f n )を導入し、上の方程式をf ( x) = 0と書く。
24この場合、ニュートン法はx(k +1) = x(k ) − H ( x(k ) )−1 f ( x(k ) )と書ける。ただし、∂f ∂ f1 ( k ) ∂ f1 ( k )(x )( x ) L 1 ( x (k ) ) ∂x n∂x2
 ∂ x1 ∂f∂f 2 ( k )∂f 2 ( k ) (k )2(x )(x ) L(x )
(k )
∂x2
∂x n
H ( x ) =  ∂ x1


L


∂f n ( k )  ∂f n ( k ) ∂f n ( k ) ∂x ( x ) ∂x ( x ) L ∂x ( x ) n2
 1

25

x(k +1) , x(k ) ,f ( x(k ) ) は n次元ベクトルでであることに注意。
もちろん、実際には、連立一次方程式

H ( x(k ) ) h = − f ( x(k ) )を解き、補正x( k +1)=x(k )+hを行い、これを繰り返す。

今日はここまで・・
26

コンピュータ科学の発展とその方向
社会科学，情報科学，・・

理学，工学，・・，バイオ・・

自然の模倣

人間の模倣

人間の思考

コンピュータ

自然現象

27