ce-12. ニュートン法による方
程式の求解
1
金子邦彦
C プログラミング応用)(全14回)
URL: https://www.kkaneko.jp/pro/c/index.html
非線形方程式 f(x) = 0 の近似解を求めるニュートン法の原
理とCプログラムによる実装、および台形則による数値積分
の基礎。
学習内容の構成
1. ニュートン法の原理:初期近似値から反復公式 x_{i+1}
= x_i - f(x_i)/f'(x_i) により解に収束させる手法
2. 収束の判定|f(x_i)| < δ となった時点で計算を終了す
る方法
3. ニュートン法の限界:初期値の選び方による収束速度の
変化、f'(x) = 0 となる点での発散、虚数解が求まらな
い制約
4. 台形則による数値積分:区間を n 等分し、小台形の面積
の和で定積分を近似する手法
前提:Cプログラミングの基礎、微分の概念(導関数、接
線)、定積分の定義
意義:解析的に解けない方程式や積分を数値的に解く手法
の習得
2
ニュートン法による求解
3
ニュートン法
ニュートン法による非線形方程式 f (x) = 0 の求
4
出力:f (x) = 0
近似解の1
入力:関数 f と,
初期近似値 x
0
再帰による繰返し計算により,x
1
x
2
... を求める
(手順は次ページ)
収束すれば,解の近似値が得られる
(注)収束しない場合もありえる
ニュートン法
初期近似値 x0 から出発
反復公式
を繰り返す
x1x2x3 ... は,徐々に解に近づいていく(と
期待できる)
5
)('
)(
1
i
i
ii
xf
xf
xx =
+
ニュートン法
初期近似値 x0 から出発
反復公式
を繰り返す
x1x2x3 ... は,徐々に解に近づいていく(と
期待できる)
6
)('
)(
1
i
i
ii
xf
xf
xx =
+
これは,点(x
i
, f (x
i
) )における 接線と,
x (y = 0) との交点 (x
i+1
, 0) を求めてい
-8
-6
-4
-2
0
2
4
6
8
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5
x
x
3
6x
2
+11x 6
7
-8
-6
-4
-2
0
2
4
6
8
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5
x
関数上の点を1つ選ぶ
関数上の点
8
-8
-6
-4
-2
0
2
4
6
8
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5
x
関数上の点
接線を引
接線
9
-8
-6
-4
-2
0
2
4
6
8
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5
x
接線
関数上の点
解の1
この場合,接線と x
軸の交点は解の1つに
近づいている
接線と x
の交点
10
-8
-6
-4
-2
0
2
4
6
8
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5
x
関数上の点 (x
0
, f (x
0
))
接線と x 軸の交
(x
0
f (x
0
)/f
'(x
0
), 0)
接線 y = f '(x
0
)(x x
0
) + f (x
0
)
関数上の点 から
接線と x 軸の交点の x 座標
を求めることを繰返す
( )
)(,
ii
xfx
)(/)('
1 iiii
xfxfxx =
+
11
-8
-6
-4
-2
0
2
4
6
8
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5
x
関数上の点 (0, 6)
x
0
= 0
接線と x 軸の交点 (0.54545454..., 0)
x
1
= 0.54545454...
12
-8
-6
-4
-2
0
2
4
6
8
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5
x
関数上の点 (0.54545454..., 1.62283996...)
x
1
= 0.54545454...
接線と x 軸の交点 (0.84895321..., 0)
x
2
= 0.84895321...
13
-8
-6
-4
-2
0
2
4
6
8
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5
x
関数上の点 (0.84895321..., 0.37398512...)
x
2
= 0.84895321...
接線と x 軸の交点 (0.97467407..., 0)
x
3
= 0.97467407...
14
-8
-6
-4
-2
0
2
4
6
8
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5
x
関数上の点 (0.97460407..., 0.052592310...)
x
3
= 0.97467407...
接線と x 軸の交点 (0.99909154..., 0)
x
3
= 0.999909154...
15
ニュートン法での収束の判定
16
(解そのものが分かっていないから)
今日の授業では次の方法で行ってみる
ある小さな正の数δに対して
)(
i
xf
となった時点で計算を終了
ニュートン法では,現在の x
i
の誤差(ど
れだけ真の解に近いか)は,正確には分か
らない
-8
-6
-4
-2
0
2
4
6
8
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5
x
f (x)
幅2δ
)(
i
xf
x
i
の値がこの範囲に入ったら計算を終了
17
ニュートン法の能力と限界
18
初期近似値 x
0
で十分近い
解を指定できれば,収束
が早い.
関数 f (x) が単調でなくて変曲点を持つ
(つまり f '(x) の符号が変わる)とき
例)上の図の1から開始
2 (f '(x) = 0となる点)が選ばれる
3 y 軸との交点が求まらない
(負の無限大に発散)
収束しない
初期近似値 x
0
の選び方
によっては,収束が遅い
ことがある
収束しないこともありえ
(右図)
x
f (x)
0
ニュートン法の注意点
虚数解は求まらない
f (x) = 0 の解が複数あっても,1回に求まる解は
1つだけ
初期近似値 x0
x0 は,コンピュータでなく,人間が決める
x0 の値によっては,収束しないこともありえる
x0 の値によって,求まる解が変わってくる(f (x) = 0
が複数の解を持つ場合)
求まる解は近似解
19
例題1.ニュートン法のプログラム
f(x)=x2 2 をニュートン法で解くプログラム
20
#include "stdio.h"
#include <math.h>
/* f(x) */
double f(double x)
{
return pow(x,2) - 2;
}
/* f(x)の導関数 */
double g(double x)
{
return 2 * x;
}
int main()
{
double x;
double new_x;
double delta;
int i;
int ch;
/* 初期値 */
x = 10;
/* 収束判定のための delta */
delta = 0.000001;
printf("繰り返し回数¥tnew_x¥t¥tf(x)¥t¥tg(x)¥n");
/* for ... になっているのは 100 回繰り返しても収束しなかったら計算を終わりたいから */
for(i = 0 ; i < 100 ; i++) {
new_x = x - f(x) / g(x);
printf("%2d¥t¥t%lf¥t%lf¥t%lf¥n",i,new_x,f(x),g(x));
/* fabs double 型の変数について絶対値を求める関数 */
if( fabs(f(new_x)) < delta ) {
printf( "収束した¥n" );
printf( "解は %lf¥n", new_x);
break;
}
x = new_x;
}
printf( "Enter キーを1,2回押してください. プログラムを終了します¥n");
ch = getchar();
ch = getchar();
return 0;
}
21
収束したかの判定を
行っている部分
反復公式
実行結果の例
22
台形則による数値積分
23
定積分
区間 [a, b] で,連続関数 f, x, x=a, x=b
囲まれた面積
24
f(x)
a
b
x
定積分:
区間 [a, b] の小区間への分割
n 個の等間隔な小区間に分割
幅: h = (b-a) / n
小区間: [x0, x1], [x1, x2], ..., [xn-1, xn]
但し,x0 = a, xi = x0 + i × h
25
f(x)
x
0
= a
x
n
= b
x
x
1
x
2
x
n-1
小台形
小台形の面積は
26
f(x)
x
0
= a
x
n
= b
x
x
i
x
i+1
小台形
f(x
i
)
f(x
i+
1
)
))()((
2
1+
+
ii
xfxf
h
台形則 trapezoidal rule
小台形の面積の和は
定積分 I を,この和 Sn で近似 台形則という
27
f(x)
x
0
= a
x
n
= b
x
))()((
2
1
1
0
+
=
+=
ii
n
i
n
xfxf
h
S
x
1
x
2
x
n-1
x
3
x
n-2
台形則
両端 x0 = a xn = b を除いて,f(xi) 2度出
28
( )
)()(
2
1
1
0
+
=
+=
ii
n
i
n
xfxf
h
S
( )
)())()((2)(
2
110 nn
xfxfxfxf
h
++++=
但し h = (b - a) / n
2回現れる部分
++=
=
)(
2
1
)()(
2
1
1
1
0 n
n
i
i
xfxfxfh
+++=
=
)(
2
1
)()(
2
1
1
1
bfihafafh
n
i
台形則による数値積分
区間[ab]を n 等分 (1区間の幅h=(b-a)/n)
n 個の台形を考え, その面積の和 Sn で,定積分 I
を近似
f(x) が連続関数のときは,n を無限大に近づければ,
Sn I に近づく
但し,単純に「n を大きくすればよい」とは言え
ない
n を大きくすると 計算時間の問題,丸め誤差の
問題が発生
29
n
n
ii
n
i
n
Sxfxf
h
I
+
=
=+=
lim))()((
2
lim
1
1
0