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 から出発
• 反復公式
を繰り返す
• x1,x2,x3 ... は,徐々に解に近づいていく(と
期待できる)
•5
)('
)(
1
i
i
ii
xf
xf
xx −=
+
ニュートン法
• 初期近似値 x0 から出発
• 反復公式
を繰り返す
• x1,x2,x3 ... は,徐々に解に近づいていく(と
期待できる)
•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
で十分近い
解を指定できれば,収束
が早い.
1
2
3
関数 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
台形則による数値積分
• 区間[a,b]を 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