数値解析(第二章) 非線形方程式を解く

数値解析 (第二章)
非線形方程式を 解く
非線形方程式と は
線形方程式と は
線形方程式の例
• 一次方程式
• 二元連立一次方程式
ax + b = 0
(
ax + by = c
dx + ey = f
• n 元連立一次方程式

a11 x1 + a12 x2 + · · · + a1n xn = b1



a x + a x + ···+ a x = b
21 1
22 2
2n n
2

···
···



an1 x1 + an2 x2 + · · · + ann xn = bn
数値解析
非線形方程式と は
非線形方程式と は 「 線形でない方程式」 である
非線形方程式の例を 挙げれば
1. 二次方程式
a x2 + b x + c = 0
2. 三次方程式
a x3 + b x2 + c x + d = 0
3. n 次方程式
an xn + an−1 xn−1 + · · · + a1 x + a0 = 0
4. そ の他
x3 + sin(x2 + x) = 0
exp(−x) + x3 = 0
数値解析
···
非線形方程式はやっ かい
そ の理由は
1. (ほと んど の場合) 解を 求める 公式が存在し な い
2. 解が幾つ存在する かすら 分から な い場合が多い
例外は n 次方程式であ る (必ず n 個の解が存在する )
n 次方程式
an xn + an−1 xn−1 + · · · + a1 x + a0 = 0
は未知数が一つである が， 解は必ず n 個 存在する 。 だが， n ≥ 5 のと き は公式が
存在し な い。
公式が存在し ない場合は計算機を 用いた数値解法を が必要になる
数値解析
最も 簡単な 数値解法
二分法 (1)
a<b
かつ
f (a) < 0, f (b) > 0
のと き ， f (x) = 0 と な る 点 α は区間 (a, b) の間に 必ず存在する 。
y = f (x)
a
b
α
数値解析
二分法 (2)
1. 微小な 数 ε > 0 と f (a) < 0, f (b) > 0 と な る 点 a < b を 選ぶ
2. 中点 c = (a + b)/2 に おいて f (c) を 計算する
3. |f (c)| ≤ ε な ら 終了し c を 解と する 。 そ う でな いな ら
4. f (c) の符号を 調べ
• f (c) > 0 な ら ば b ← c (左下図)
• f (c) < 0 な ら ば a ← c (右下図)
と し て ， 2. へ戻る
こ の方法を 用いる と ， 反復一回毎に ， 解 α が存在する 区間の幅が半分ずつ減っ
て いく → 表 2.1
a
数値解析
c
b
a
c
b
二分法 (3)
問 次のプロ グラ ム の空白 (1), (2), (3) を 埋めよ 。 た だし ， f (a) < 0, f (b) > 0 と
する 。
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
数値解析
#include <stdio.h>
#include <math.h>
double f(double x);
main ()
{
double a=1.0,b=2.0,c,eps=0.0001;
c=(a+b)/2.0;
while (
(1)
) {
if (f(c)>0)
( (2) );
else
( (3) );
c=(a+b)/2.0;
}
printf("解は %f\n",c);
}
二分法 (4)
二分法が使え ない場合
1. 関数 f (x) が不連続な 場合 (左下図)
2. 解 α の重複度が偶数の場合 (右下図)
例え ば
f (x) = (x − α)2
数値解析
二分法 (5)
二分法を 使う に あ た っ て
1. 初めに ， a < b かつ f (a) < 0, f (b) > 0 と な る 点を 見つけ る
逆の場合， すな わち f (a) > 0, f (b) < 0 の場合， −f を f と 置き 換え る
2. 反復停止条件は |f (c)| < ε と する (ε > 0 は小さ い数)
3. ε > 0 を あ ま り 小さ く する と 止ま ら な く な る 恐れがあ る
4. 高精度計算向き ではな い
数値解析
Newton 法
Newton 法の原理 (1)
点 x0 に おけ る y = f (x) の接線の方程式は
y − f (x0 ) = f ′ (x0 ) (x − x0 )
であ る 。 こ の接線が x-軸と 交わる 点を x1 と する 。
f (x0 )
x1 = x0 − ′
f (x0 )
x1 でも 同じ こ と を する 。 すな わち ， 点 (x1 , f (x1 )) に おけ る y = f (x) の接線と
x-軸と の交点を 求め x2 と する 。
α
y = f (x)
数値解析
x2 x1
Newton 法の原理
x0
Newton 法の原理 (2)
要する に ， Newton 法と は
xk+1
f (xk )
,
= xk − ′
f (xk )
k = 0, 1, 2, . . .
と いう 漸化式 を 用いて ， 解 α の近似値 xk を 逐次更新 し て いく ア ルゴ リ ズム で
ある 。
前ページの図の様な 状態では， x0 , x1 , x2 , · · · と 急速に 解 α に 近づいて いく (収
束し て いく )。
数値解析
Newton 法のプロ グ ラ ム
問 次のプロ グラ ム の空白を 埋めよ 。
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
数値解析
#include <stdio.h>
#include <math.h>
double f(double x); /* f(x) */
double d(double x); /* f’(d) */
main ()
{
double x0=2.0,x1,eps=0.0001;
while (fabs(f(x0))>eps) {
x1=x0-f(x0)/d(x0);
(
);
}
printf("解は %f\n",x0);
}
二分法と Newton 法
Newton 法の方が圧倒的に 速い！
k
0
1
2
3
4
5
6
7
26
x2 − 2 = 0 を 解いた 例
二分法
Newton 法
|f (c)|
|f (xk )|
2.500000e-01 2.000000e+00
4.375000e-01 2.500000e-01
1.093750e-01 6.944444e-03
6.640625e-02 6.007305e-06
2.246094e-02 4.510835e-12
2.172852e-02
4.272461e-04
1.063538e-02
···
5.236811e-09
二分法では初期値 a = 1, b = 2， Newton 法では初期値 x0 = 2 と し た
数値解析
グ ラ フ で表す と
1e+00
|f (x)| 1e-02
1e-04
Newton 法
1e-06
二分法
1e-08
1e-10
1e-12
1e-14
1
10
k
数値解析
100
Newton 法と は
1. f (x) の他に f ′ (x) も 計算し な け ればな ら な い
2. (う ま く いけ ば ) 収束は速い
3. 初期値の与え 方が悪いと 暴走する こ と も あ る (例 2.3)
4. 暴走を 食い止める た め， 最大反復回数を 決めて おく
数値解析
暴走を 食い止める には (1)
問 以下のプロ グラ ム は最大 k_max 回ま でし か反復し な いよ う に 作ら れて いる
(f (x), f ′ (x) の定義は省略)。 こ のプロ グラ ム が正確に 動く よ う に 11 行目
の空白を 埋めよ (実行例は次ページ )。
数値解析
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
#include <stdio.h>
#include <math.h>
double f(double x); /* f(x) */
double d(double x); /* f’(d) */
main ()
{
int k_max=10,k=0;
double x0=1.0, x1, eps=1.e-6;
printf(" %d %f\n",k,x0);
while ((fabs(f(x0))>eps) && (
)
) {
x1=x0-f(x0)/d(x0);
x0=x1;
k=k+1;
printf(" %d %f\n",k,x0);
}
printf("\n");
if (k==k_max)
printf(" 収束し ま せんでし た !! \n");
else {
printf(" %d 回の反復で収束し ま し た 。 \n",k);
printf(" 解は %f です。 \n",x0);
}
}
暴走を 食い止める には (2)
実行結果
0 1.000000
1 0.000000
2 1.000000
3 0.000000
4 1.000000
5 0.000000
6 1.000000
7 0.000000
8 1.000000
9 0.000000
10 1.000000
収束し ま せんでし た !!
数値解析
0
1
2
3
4
-1.500000
-1.842105
-1.772827
-1.769301
-1.769292
4 回の反復で収束し ま し た 。
解は -1.769292 です。
セ カ ン ト 法 (secant method)
セカ ン ト 法 (1)
secant: 割線 (2 点， あ る いはそ れ以上の点で曲線と 交わる 直線)
Newton 法では f (x) の他に f ′ (x) も 計算し な く て はな ら な い。 式が複雑で
f ′ (x) の計算がやっ かいな と き はど う する ？
a ≈ b のと き
f (a) − f (b)
f (a) ≈
a−b
′
であ る から ， こ の関係式を 用いて f ′ (x) を 近似する 。
そ う する と ， Newton 法は， xk と xk−1 が近いと き
xk+1 = xk −
≈ xk −
f (xk )
f ′ (xk )
f (xk )
f (xk )−f (xk−1 )
xk −xk−1
と な る 。 上式で ≈ を = に 置き 換え た のがセ カ ン ト 法であ る 。
数値解析
セカ ン ト 法 (2)
実際に 計算する と き は
xk+1
xk − xk−1
= xk −
f (xk ),
f (xk ) − f (xk−1 )
k = 1, 2, . . .
と いう 式を 用いる 。 セ カ ン ト 法では初期値が 2 つ必要に な る (x0 と x1 )。
と こ ろ で， ２ 点 (x0 , f (x0 ))， (x1 , f (x1 )) を 結ぶ直線の方程式は
y=(
) であ り ， そ の直線の x–軸と の交点の座標を x2 と すれば
x2 = (
) と なる 。
数値解析
y = f (x)
x0
x2
x1
α
セ カ ン ト 法の原理
問 上の図で x3 , x4 はど こ に な る か。 図示せよ 。
数値解析
セカ ン ト 法のプロ グ ラ ム (1)
問 次のセ カ ン ト 法のプロ グラ ム が正常に 動く よ う に 空白 (1)∼(4) を 埋めよ 。
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
数値解析
#include <stdio.h>
#include <math.h>
double f(double x);
main()
{
int k;
double eps=1.0e-10, x0=1.2, x1=1.5, x2, f0, f1, f2;
k=0; f0=f(x0);
printf("%d %e
k=1; f1=f(x1);
printf("%d %e
%e \n",k,x0,f0);
%e \n",k,x1,f1);
while( fabs(f1) > eps) {
x2=x1-f1*(x1-x0)/(f1-f0);
f2=f(x2);
k=k+1;
printf("%d %e %e \n",k,x2,f2);
(1) ; (2) ;
(3) ; (4) ;
}
}
double f(double x)
{
return x*x-2.0;
}
セカ ン ト 法のプロ グ ラ ム (2)
実行結果
0
1
2
3
4
5
数値解析
1.200000e+00
1.500000e+00
1.407407e+00
1.414013e+00
1.414214e+00
1.414214e+00
-5.600000e-01
2.500000e-01
-1.920439e-02
-5.679744e-04
1.370231e-06
-9.729584e-11
まとめ
セカ ン ト 法と は
• f ′ (x) を 計算し な く て も よ い
• 初期値が 2 つ必要
• Newton 法と 同様に 暴走する こ と も あ る
• Newton 法よ り 収束が少し 遅い
• し かし 関数の評価回数は１ ス テ ッ プあ た り １ 回
数値解析
代数方程式の解法
代数方程式と は
代数方程式と は
an xn + an−1 xn−1 + · · · + a1 x + a0 = 0
と いう 方程式であ る 。
• こ の方程式の解は， x を 複素数 に 拡張する と ， 必ず n 個存在する
• し かし ， n ≥ 5 では解を 求める 公式が存在し な い
• 解を 求める 公式が存在し て も 数値解法を 用いた 方が簡単な こ と も あ る
• 数値解法は Newton 法がよ く 用いら れる
数値解析
多項式の計算法
n 次多項式
は
Pn (x) = an xn + an−1 xn−1 + · · · + a1 x + a0
(· · · (((an x + an−1 ) x + an−2 ) x + an−3 ) x + · · · + a1 ) x + a0
と 書き 表せる 。 例え ば
3 x3 + 2 x2 − 2 x + 1 = ((3 x + 2) x − 2) x + 1
であ る 。 こ のよ う な 計算法は ホーナー (Horner) 法 と 呼ばれて いる 。
問 次の多項式を ホーナー法で表せ。
4 x4 − 2 x3 + x2 − 2 x + 1,
数値解析
5 x4 + 2 x3 − 3 x + 4
ホーナ ー法の計算量
Pn (x) を ま と も ？ に 計算する と ， ai xi の計算で乗算 i 回必要と する ので， 乗算
は合計
n (n + 1)
1 + 2 + ··· + n =
回
2
と なる 。
ホーナ ー法のプロ グ ラ ム
P=a[n];
for (k=n-1; k>=0; k--)
P=P*x+a[k];
ホーナー法では乗算は n 回
数値解析
xn の計算法 (1)
xn は何回の乗算が必要か
ま と も に 計算すれば x100 は乗算が 99 回必要に な る 。 では xn は n − 1 回必
要か ?
x100 は
x100 = x64 · x32 · x4
であ る から ， x64 ， x32 ， x4 の 3 つがあ ればよ い。 こ の３ つは
x → x2 → x4 → x8 → x16 → x32 → x64
と いう よ う に 二乗演算を 6 回繰り 返せば得ら れる 。 し たがっ て 合計 6+2=8 回 で
x100 が計算でき る こ と に な る 。
問 こ のよ う な 方法を 用いる と ， 下のそ れぞれを 計算する のに 何回の乗算が必
要か。
数値解析
x38 ,
x78 ,
x47
xn の計算法 (2)
こ の計算法は ロ シア の農民算法 と 呼ばれて いる 。 こ の方法は， 本来はベキ xn の
計算でな く ， 乗算の計算法 と し て 考え 出さ れた 。
例え ば
37 × 45 = 1665
と いう 計算は， 37=32+4+1 であ る から
(32 + 4 + 1) × 45 = 32 · 45 + 4 · 45 + 45
と な り ， 4 · 45 と 32 · 45 は 45 を
45 → 90 → 180 → 360 → 720 → 1440
のよ う に ， 次々 と 2 倍に し て いく (同じ も のを 加え て いく ) と いう 操作を 5 回繰
り 返し ， 必要な と こ ろ (赤い部分) だけ を 保存し て おけ き ， そ れら を 掛け れ足せ
ばよ い。 し た がっ て 足し 算だけ で， し かも た っ た 5+2=7 回の足し 算で 37 × 45
が求ま る 。
数値解析
ロ シ ア の農民算法のプロ グ ラ ム
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
/*
xˆn の計算
*/
#include <stdio.h>
int pwr(int x, int n);
main()
{
int x,n;
}
scanf("%d",&x);
scanf("%d",&n);
printf("%dˆ%d=%d \n",x,n,pwr(x,n));
int pwr(int x, int m)
{
int p,y=x;
int mm;
if (m == 0) return (1);
p=1.0; mm=m;
do {
if ((mm % 2)==1) p=p*y;
y=y*y;
mm=mm/2;
} while (mm >0);
return (p) ;
}
3
13
3ˆ13=1594323
数値解析
組み立て 除法
ホーナー法の計算は， 手計算では 組み立て 除法 を 用いる と わかり やすい：
例 P (x) = 3 x3 − 2 x2 + x − 1 のと き ， ホーナー法に よ っ て P (2) を 計算する と
P (2) = ((3 · 2 − 2) · 2 + 1) · 2 − 1 = 17
であ る から ， 以下のよ う な 組み立て 除法
3
+
3
-2
1
-1
6
8
18
4
9
17
2
P (2)
を 行っ て いる こ と に な る (ր は ×2 を 表し て いる )。
問 1 P (x) = 2 x4 + x3 − 2 x2 + x + 1 のと き P (−1) を 求めよ 。
問 2 P (x) = 3 x4 − x3 + 2 x − 4 のと き P (2) を 求めよ 。
数値解析
組み立て 除法と 割り 算の関係
P (x) を x − a で割っ た と き ， 商を Q(x)， 余り を R と する と
P (x) = (x − a) Q(x) + R
であ る から ， P (a) = R と な る → P (a) は P (x) を x − a で割っ た 余り
よ く 見れば， 組み立て 除法は割り 算と 同じ こ と を 行っ て いる 。
P (x) = 3 x3 − 2 x2 + x − 1 のと き ， こ れを x − 2 で割れば
3 x2 +4 x +9
x − 2 3 x3 −2 x2 +x −1
3 x3 −6 x2
4 x2
+x
4 x2 −8 x
9x
−1
9 x −18
17
数値解析
組み立て 除法を 用いて P ′ (x) も 計算でき る (1)
前ページから ， 組み立て 除法に よ っ て 余り R だけ でな く ， 商 Q(x) の係数 も 求
ま っ て いる こ と がわかる 。
と こ ろ で， P (x) = (x − a) Q(x) + R を 微分する と
P ′ (x) = Q(x) + (x − a) Q′ (x)
であ る から ， P ′ (a) = Q(a) であ る → Q(x) に も 組み立て 除法を 適用すれば
P ′ (a) も 求ま る
数値解析
組み立て 除法を 用いて P ′ (x) も 計算でき る (2)
P (x) = 3 x3 − 2 x2 + x − 1 のと き ， P ′ (2) = 29 は， 以下のよ う に 組み立て 除法
を 2 回 実行すれば求ま る 。
3
+
Q(x) の係数
3
+
3
-2
1
-1
6
8
18
4
9
17
6
20
10
29
2
P (2)
P ′ (2)
問 P (x) = 2 x4 − 3 x3 + x2 − 2 x − 1 のと き P (−1) と P ′ (−1) を 求めよ 。
数値解析
P (x) と P ′ (x) を 同時に計算す る プロ グ ラ ム
P (x) と P ′ (x) を 同時に計算す る プロ グ ラ ム
1: #include <stdio.h>
2: #define n 3
3: main ()
4: {
5:
int i;
6:
double x,a[n+1],R,T;
7:
8:
a[3]=3.0; a[2]=-2.0; a[1]=1.0; a[0]=-1.0;
9:
scanf("%lf",&x);
10:
11:
R=a[n]; T=a[n];
12:
for (i=n-1; i>0; i--) {
13:
R=R*x+a[i];
14:
T=T*x+R;
15:
}
16:
R=R*x+a[0];
17:
printf("P(%3.2f)=%5.2f \n", x,R);
18:
printf("P’(%3.2f)=%5.2f \n", x,T);
19: }
2.0
P(2.00)=17.00
P’(2.00)=29.00
数値解析
組み立て 除法を 用いた Newton 法のプロ グ ラ ム
/*
/*
数値解析
P_n(x)=a_n xˆn+ ・・・ +a_0=0 の解を 求める
Newton 法のプロ グラ ム */
main ()
{
......
k=0;
eps=....;
x0= ....;
R= P_n(x0) を 計算する か， eps よ り 大き な 値を 設定する
R= ....;
while ( fabs(R) > eps ) {
組み立て 除法のプロ グラ ム
（ 前ページ 11 行---16 行）
x1=x0-R/T;
x0=x1;
k=k+1;
printf(" k= %d x1 = %f \n",k,x1);
}
.....
}
*/
プロ グ ラ ム (1)
x4 − x3 − x2 − x − 1 = 0 の解を 求める Newton 法のプロ グラ ム
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
数値解析
#include <stdio.h>
#include <math.h>
#define n 4
main () {
double a[n+1],R,T,x0,x1,eps=1.0e-10;
int i,k=0;
a[4]=1.0; a[3]=-1.0; a[2]=-1.0; a[1]=-1.0; a[0]=-1.0;
x0=1.5;
R=1000.0;
while (fabs(R)>eps) {
R=a[n]; T=a[n];
for (i=n-1; i>0; i--) {
R=R*x0+a[i];
T=T*x0+R;
}
R=R*x0+a[0];
x1=x0-R/T;
x0=x1;
k=k+1;
printf(" %d %f %e \n",k,x1,R);
}
}
プロ グ ラ ム (2)
実行結果
1
2
3
4
5
6
7
8
数値解析
2.613636
2.202741
1.992124
1.932199
1.927588
1.927562
1.927562
1.927562
-3.062500e+00
1.836513e+01
4.799914e+00
8.829103e-01
5.896841e-02
3.310245e-04
1.062681e-08
1.304512e-15
Newton 法の誤差 (1)
f (x) = 0 の解を α と する 。 すな わち f (α) = 0 と する 。 こ のと き Newton 法の
漸化式
f (xk )
xk+1 = xk − ′
f (xk )
の両辺から α を 引く と
f (xk )
xk+1 − α = xk − α − ′
f (xk )
と な る 。 xk − α は xk と 真値 α と の差 (誤差) であ り ， そ れを ek と 表す。 そ う
する と
f (α + ek )
ek+1 = ek − ′
f (α + ek )
と な る 。 こ こ で f (α + ek ) は Taylor 展開よ り
1 ′′
1 ′′
′
2
f (α + ek ) ≈ f (α) + f (α) ek + f (α) ek = f (α) ek + f (α) e2k
2
2
′
数値解析
Newton 法の誤差 (2)
さ ら に ， f ′ (α) は
と な る 。 こ れら よ り
f ′ (α + ek ) ≈ f ′ (α) + ek f ′′ (α)
ek+1
ek f ′ (α) + 12 e2k f ′′ (α)
≈ ek −
f ′ (α) + ek f ′′ (α)
f ′′ (α) 2
e
≈
2 f ′ (α) k
を 得る 。
こ れよ り ， Newton 法では
次回の誤差 ≈ 定数 × (今回の誤差) 2
な ので， 収束が速い。 た だし f ′ (α) 6= 0 のと き 。 ⇒ 二次収束
数値解析
Newton 法はいかに凄いか
f (x) = x2 − 2 = 0 に Newton 法を 適用し た 場合 (x0 = 2)
√
k ek (= xk − 2 ) ek /e2k−1
0 5.8578 e-0001
1 8.5786 e-0002 0.25000
2 2.4531 e-0003 0.33333
3 2.1239 e-0006 0.35294
4 1.5949 e-0012 0.35355
5 8.9929 e-0025 0.35355
6 2.8593 e-0049 0.35355
7 2.8905 e-0098 0.35355
8 2.9539 e-0196 0.35355
9 3.0849 e-0392 0.35355
10 3.3647 e-0784 0.35355
11 4.0026 e-1568 0.35355
12 5.6641 e-3136 0.35355 ⇐ (1/2)f ′′ (α)/f ′ (α)
‘e-’ の後の数値がほぼ倍々 で増加し て いる ！
数値解析