f(x)

.
Ｃプログラミング (2014)
.
第7回
− 方程式の解：Newton 法 −
松田七美男
2014 年 11 月 13 日
. 方程式 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 )| < ε
☞上記の条件を満たさない場合に備えて，反復回
数の上限を定めておくことも必要
. ニュートン法：概要
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 (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 (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 (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 (xk )
(k = 0, 1, · · · )
. ニュートン法：原理
xk+1 − xk = h
1 とおいて，点 xk で Taylor 展開
し h の 2 次以上の項を無視すれば，
0
f (xk+1 )
= f (xk + h)
f (xk ) + f (xk )h
= f (xk ) + f (xk )(xk+1 − xk )
f (xk )
∴ xk+1 = xk −
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
π
4
0
−
π
4
-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
π
4
0
−
π
4
-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
π
4
0
−
π
4
-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
π
4
0
−
π
4
-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
π
4
0
−
π
4
-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
π
4
0
−
π
4
-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 法は初期値に敏感
π
4
0
π
−
4
-0.2
-0.4
-1.5
-1
-0.5
0
x
.
0.5
1
1.5
. 数値微分
Newton 法の改善
. 数値微分
Newton 法の改善
. 初期値に敏感 ⇒ 人間が判断するしかない？
1
. 数値微分
Newton 法の改善
. 初期値に敏感 ⇒ 人間が判断するしかない？
2. f (x) の導関数を式で与えるのは面倒 ⇒ 数値微分
1
例えば，次の近似式を用いる．
.
f (x + ∆x) − f (x)
f (x)
∆x
.
(1)
. 数値微分
Newton 法の改善
. 初期値に敏感 ⇒ 人間が判断するしかない？
2. f (x) の導関数を式で与えるのは面倒 ⇒ 数値微分
1
例えば，次の近似式を用いる．
.
f (x + ∆x) − f (x)
f (x)
∆x
.
(1)
素直なコーディング
関数へのポインタ
double df(double (*f)(), double x, double dx)
{
return (f(x+dx)-f(x))/dx ;
}
.
√
. arccos (1/ 2) の値
GUI 関数電卓ソフト（xcalc 等）を使わなくともコマンドラ
イン上で計算させることができる．
octave（MATLAB 互換）
.
octave
-q --eval ’acos(1/sqrt(2))’
.
gnuplot（グラフ作成ツールだが数学関数も豊富）
.
gnuplot
-e ’print acos(1/sqrt(2))’
.
☞
1
π
f (x) = 0 の点は，x = ± arccos √ = ± でした．
4
2
. Newton 法よりも高次の収束
Halley の式（3 次収束）
公式： xk+1 = xk −
fk /fk
1 − fk fk /2(fk )2
2 階微分の近似式：
f (x + ∆x) − f (x)
∆x
f (x + 2∆x) − 2f (x + ∆x) + f (x)
∼
(∆x)2
f (x) ∼
を用いて，数値微分を行うことができる．

Download Report