.
Ｃプログラミング (2014)
.
第6回
−方程式の解：2 分法，はさみうち法 −
松田七美男
2014 年 11 月 6 日
. 方程式 f (x) = 0 の解法
解析的な公式がない場合反復法のような近似計算
で求める他ない
.
反復解法
.
初期値 x0 を設定し，x1 , x2 , · · · を関係式
xk+1 = g(xk ) (k = 0, 1, 2, · · · )
により，逐次計算する方法．
.
. 方程式 f (x) = 0 の解法
解析的な公式がない場合反復法のような近似計算
で求める他ない
.
反復解法
.
初期値 x0 を設定し，x1 , x2 , · · · を関係式
xk+1 = g(xk ) (k = 0, 1, 2, · · · )
により，逐次計算する方法．
.
2 分法，はさみうち法，Newton 法
. 収束判定
ある微小な値 ε > 0 を設定して，例えば以下のよ
うな条件を満たしたときに収束したとみなして，
計算を打ちきる．
.
収束判定条件の例
.
|xk+1 − xk | < ε
.
|xk+1 − xk | < ε|xk |
|f (xk+1 )| < ε
. 収束判定
ある微小な値 ε > 0 を設定して，例えば以下のよ
うな条件を満たしたときに収束したとみなして，
計算を打ちきる．
.
収束判定条件の例
.
|xk+1 − xk | < ε
.
|xk+1 − xk | < ε|xk |
|f (xk+1 )| < ε
+上記の条件を満たさない場合に備えて，反復回
数の上限を定めておくことも必要
. 2 分法：概要
k 回の操作の後の区間 [ xLk : xUk ] の幅を δk とすると，
xU0
xL0
δ0
.
. 2 分法：概要
k 回の操作の後の区間 [ xLk : xUk ] の幅を δk とすると，
xL0
xU0
.
. 2 分法：概要
k 回の操作の後の区間 [ xLk : xUk ] の幅を δk とすると，
xU0
xL0
xL1
xU1
δ1= 1δ0
2
.
. 2 分法：概要
k 回の操作の後の区間 [ xLk : xUk ] の幅を δk とすると，
xU0
xL0
xL1
xU1
.
. 2 分法：概要
k 回の操作の後の区間 [ xLk : xUk ] の幅を δk とすると，
xU0
xL0
xL1
xU1
xL2
xU2
δ 2= 1 δ1= 1 δ0
2
4
.
. 2 分法：概要
k 回の操作の後の区間 [ xLk : xUk ] の幅を δk とすると，
xU0
xL0
xL1
xU1
xL2
xU2
.
. 2 分法：概要
k 回の操作の後の区間 [ xLk : xUk ] の幅を δk とすると，
xU0
xL0
xL1
xU1
xL2
xL3
xU2
xU3
.
. 2 分法：概要
k 回の操作の後の区間 [ xLk : xUk ] の幅を δk とすると，
xU0
xL0
xL1
xU1
xL2
xL3
xU2
xU3
( )k
1
δk =
δ0
2
.
. 2 分法：概要
k 回の操作の後の区間 [ xLk : xUk ] の幅を δk とすると，
xU0
xL0
xL1
xU1
xL2
xL3
xU2
xU3
( )k
1
δk =
δ0
2
.
.いずれ必ず収束する
.
. 2 分法：コーディング
do {
dx = xU - xL;
x = (xU + xL)/2;
if (f(x)*f(xU) > 0) xU = x; else xL = x;
} while (dx > eps);
. 2 分法：コーディング
.
do {
dx = xU - xL;
x = (xU + xL)/2;
if (f(x)*f(xU) > 0) xU = x; else xL = x;
} while (dx > eps);
.
真値の両側に区間の端がある条件を保つ
.
f (xL )f (xU ) < 0
.
(1)
. y = f (x) の大域的な挙動
0.4
f (x) =
1
tan(x) − x
2
y
0.2
( √ )
arccos 1/ 2
( √ )
− arccos 1/ 2
0
-0.2
-0.4
-1.5
-1
-0.5
0
x
.
0.5
1
1.5
. y = f (x) の大域的な挙動
0.4
f (x) =
1
tan(x) − x
2
y
0.2
( √ )
arccos 1/ 2
( √ )
− arccos 1/ 2
0
解のおおよその位置
-0.2
-0.4
-1.5
-1
-0.5
0
x
.
0.5
1
1.5
. y = f (x) の大域的な挙動
0.4
f (x) =
1
tan(x) − x
2
y
0.2
( √ )
arccos 1/ 2
( √ )
− arccos 1/ 2
0
解のおおよその位置
-0.2
-0.4
-1.5
-1
-0.5
0
0.5
1
1.5
x
.
この図の確認の方法：gnuplot -p -e ’plot [-1.5:1.5][-0.5:0.5] 0.5*tan(x)-x’
. はさみうち法：概要
y = f (x)
x1
.
.
xk+1 = xk −
x
x2 x3 xk k+1
x0
x0 − xk
xk f (x0 ) − x0 f (xk )
f (xk ) =
f (x0 ) − f (xk )
f (x0 ) − f (xk )
. はさみうち法：コーディング
f0 = f(x0);
fp = f(xp);
x = (xp*f0 - x0*fp)/(f0 - fp);
fx = f(x);
if ( fabs(x-xp) < DBL EPSILON
|| fabs(fx) < DBL EPSILON ) break;
if ( fx*f0 > 0 ) x0 = xp;
xp = x;
.
. はさみうち法：コーディング
f0 = f(x0);
fp = f(xp); x 軸との交点の算出
x = (xp*f0 - x0*fp)/(f0 - fp);
fx = f(x);
if ( fabs(x-xp) < DBL EPSILON
|| fabs(fx) < DBL EPSILON ) break;
if ( fx*f0 > 0 ) x0 = xp;
xp = x;
.
. はさみうち法：コーディング
f0 = f(x0);
fp = f(xp); x 軸との交点の算出
x = (xp*f0 - x0*fp)/(f0 - fp);
fx = f(x);
if ( fabs(x-xp) < DBL EPSILON
|| fabs(fx) < DBL EPSILON ) break;
if ( fx*f0 > 0 ) x0 = xp;
xp = x;
真値の両側に区間の端がある条件を保つ
.
. 修正はさみうち法：概要
y = f (x)
1 f (x )
0
2
1 f (x )
0
22
x1
x2
x3
x4
x0
.
傾きを小さくする工夫
.
1
1
f0 → f0 → 2 f0 → · · ·
2
2
.
. 修正はさみうち法：コーディング
fp = f(xp);
x = (xp*f0 - x0*fp)/(f0 - fp);
fx = f(x);
if ( fabs(x-xp) < DBL EPSILON
|| fabs(fx) < DBL EPSILON ) break;
if ( fx*f0 > 0 ) {
x0 = xp;
f0 = f(x0);
} else f0 = f0/2;
xp = x;
.
. 修正はさみうち法：コーディング
fp = f(xp);
x = (xp*f0 - x0*fp)/(f0 - fp);
fx = f(x);
if ( fabs(x-xp) < DBL EPSILON
|| fabs(fx) < DBL EPSILON ) break;
if ( fx*f0 > 0 ) {
x0 = xp;
f0 = f(x0);
} else f0 = f0/2;
xp = x;
傾きを小さくする工夫
.
. ニュートン法：概要
y = f(x)
f (xk)
xk+2
xk+1
xk
f (xk)
f ’(xk) = x - x
k
k+1
. ニュートン法：概要
y = f(x)
f (xk)
xk+2
xk+1
xk
f (xk)
f ’(xk) = x - x
k
k+1
. ニュートン法：概要
y = f(x)
f (xk)
xk+2
xk+1
xk
f (xk)
f ’(xk) = x - x
k
k+1
.
.
xk+1 = xk −
f (xk )
f 0 (xk )
(k = 0, 1, · · · )
. ニュートン法：概要
y = f(x)
f (xk)
f (xk+1)
xk+2
xk+1
xk
f (xk)
f ’(xk) = x - x
k
k+1
.
.
xk+1 = xk −
f (xk )
f 0 (xk )
(k = 0, 1, · · · )
. ニュートン法：概要
y = f(x)
f (xk)
f (xk+1)
xk+2
xk+1
f (xk+1)
f ’(xk+1) = x - x
k+1
k+2
xk
.
.
xk+1 = xk −
f (xk )
f 0 (xk )
(k = 0, 1, · · · )
. ニュートン法：概要
y = f(x)
f (xk)
xk+2
xk+1
xk
f (xk)
f ’(xk) = x - x
k
k+1
.
.
xk+1 = xk −
f (xk )
f 0 (xk )
(k = 0, 1, · · · )
. ニュートン法：原理
xk+1 − xk = h 1 とおいて，点 xk で Taylor 展開
し h の 2 次以上の項を無視すれば，
0 ' f (xk+1 )
= f (xk + h)
' f (xk ) + f 0 (xk )h
= f (xk ) + f 0 (xk )(xk+1 − xk )
f (xk )
∴ xk+1 = xk − 0
f (xk )
. ニュートン法：アルゴリズム
for (i = 0; i < N; i++) {
dx = -f(x)/df(x);
x += dx;
if (fabs(dx) <= eps) {
printf("x = %.12e", x);
}
printf("Can’t solve\n");
.
. ニュートン法：アルゴリズム
N:繰り返しの上限
for (i = 0; i < N; i++) {
dx = -f(x)/df(x);
x += dx;
if (fabs(dx) <= eps) {
printf("x = %.12e", x);
}
printf("Can’t solve\n");
.
. ニュートン法：アルゴリズム
N:繰り返しの上限
for (i = 0; i < N; i++) {
dx = -f(x)/df(x); df(x):微分値を返す関数
x += dx;
if (fabs(dx) <= eps) {
printf("x = %.12e", x);
}
printf("Can’t solve\n");
.
. ニュートン法：アルゴリズム
N:繰り返しの上限
for (i = 0; i < N; i++) {
dx = -f(x)/df(x); df(x):微分値を返す関数
x += dx;
if (fabs(dx) <= eps) { 打ち切りの判定
printf("x = %.12e", x);
}
printf("Can’t solve\n");
.
. ニュートン法：アルゴリズム
N:繰り返しの上限
for (i = 0; i < N; i++) {
dx = -f(x)/df(x); df(x):微分値を返す関数
x += dx;
if (fabs(dx) <= eps) { 打ち切りの判定
printf("x = %.12e", x); 成功：結果の表示
}
printf("Can’t solve\n");
.
. ニュートン法：アルゴリズム
N:繰り返しの上限
for (i = 0; i < N; i++) {
dx = -f(x)/df(x); df(x):微分値を返す関数
x += dx;
if (fabs(dx) <= eps) { 打ち切りの判定
printf("x = %.12e", x); 成功：結果の表示
}
printf("Can’t solve\n"); 失敗の旨の表示
.
. ニュートン法：初期値の適正範囲
0.4
f (x) =
1
tan(x) − x
2
y
0.2
( √ )
arccos 1/ 2
( √ )
− arccos 1/ 2
0
-0.2
-0.4
-1.5
-1
-0.5
0
x
.
0.5
1
1.5
. ニュートン法：初期値の適正範囲
0.4
f (x) =
1
tan(x) − x
2
y
0.2
( √ )
arccos 1/ 2
( √ )
− arccos 1/ 2
0
-0.2
-0.4
-1.5
-1
-0.5
0
x
.
0.5
1
1.5
. ニュートン法：初期値の適正範囲
0.4
f (x) =
1
tan(x) − x
2
y
0.2
( √ )
arccos 1/ 2
( √ )
− arccos 1/ 2
0
-0.2
-0.4
-1.5
-1
-0.5
0
x
.
0.5
1
1.5
. ニュートン法：初期値の適正範囲
0.4
f (x) =
1
tan(x) − x
2
y
0.2
( √ )
arccos 1/ 2
( √ )
− arccos 1/ 2
0
-0.2
-0.4
-1.5
-1
-0.5
0
x
.
0.5
1
1.5
. ニュートン法：初期値の適正範囲
0.4
f (x) =
1
tan(x) − x
2
y
0.2
( √ )
arccos 1/ 2
( √ )
− arccos 1/ 2
0
-0.2
-0.4
-1.5
-1
-0.5
0
x
.
0.5
1
1.5
. ニュートン法：初期値の適正範囲
0.4
f (x) =
1
tan(x) − x
2
y
0.2
( √ )
arccos 1/ 2
( √ )
− arccos 1/ 2
0
-0.2
-0.4
-1.5
-1
-0.5
0
x
.
0.5
1
1.5
. ニュートン法：初期値の適正範囲
0.4
f (x) =
1
tan(x) − x
2
0.2
y
Newton 法は初期値に敏感
( √ )
arccos 1/ 2
( √ )
− arccos 1/ 2
0
-0.2
-0.4
-1.5
-1
-0.5
0
x
.
0.5
1
1.5