求根アルゴリズムについて
1
14年10月22日水曜日
求根アルゴリズム (Root Finding Algorithms)
Idea
与えられた一次元の関数 f (x) に対して
⇤
0
の f (x ) = となる点を求める．
例
f (x)
2
14年10月22日水曜日
注意事項（Observations）
二つの根の場合 f (x) = x2 x 1
根無しの場合
f (x) = tanh(x)
3
実験からのデータはどうするのか？
3
14年10月22日水曜日
アルゴリズム的なアプローチ
何が必要か？
• 初期値: x0
xn+1 , xn , ...
• の関係
{x0 , x1 , ...}
• 終了基準（stopping criteria）
1. |xn+1
xn | < tol?
2. |f (xn+1 )
f (xn )| < tol?
3. Both?
4
14年10月22日水曜日
基本的な手法 (standard methods)
• Bisection (二分法)
– 単純で使い易い反復法の一つ
• Newton’s Method (ニュートン法)
– 関数の微分が必要
– ２次収束する (非常に早く解に収束)
• Secant Method （セカント法）
– 微分の代用としてセカント(secant)の線を使う
– 超線形収束する (Converges super-linearly)
5
14年10月22日水曜日
(二分法) Bisection Method
• 中間値の定理 (intermediate value theorem)
f (x)
⇤
⇤
9x s.t. f (x ) = 0
6
14年10月22日水曜日
(二分法) Bisection Method
Step 1
f (x)
+
-
7
14年10月22日水曜日
(二分法) Bisection Method
Step 2
+
左辺は他の解があるが，
中間値の定理により，右
それらを必要としない．
辺の間にはf(x)＝０とな
る点があるので，この区
間に着目する．
14年10月22日水曜日
8
(二分法) Bisection Method
Step 3
Step 4
+
-
+
-
ここ！
繰り返す
9
14年10月22日水曜日
Exercise
p
1. 二分法を用いての値を近似せよ．
2
また，各反復回数ごとに以下の誤差も表示せよ．
i.
ii.
Note: |ak
|ak
bk |
|f (ak )
bk | =
f (bk )|
1 k
2
|a0
b0 |
10
14年10月22日水曜日
二分法の疑似コード（Pseudo code）
fa
f (a)
for i = 1 to M AX do
m
a + (b a)/2;
fm
f (m)
if fa fm > 0 then
a
m
fa
fm
else
b
m
end
end
return m
11
14年10月22日水曜日
実験からのデータ場合はどうするのか？
データファイルを読み込む！
=) 配列 (arrays)
12
14年10月22日水曜日
配列 (Arrays)
配列の宣言(declaration)
double a[7];
a:
a[0]
a[1]
a[2]
a[3]
a[4]
a[5]
a[6]
初期化(initialization)
a[3] = 1.3;
1.3
13
14年10月22日水曜日
配列 (Arrays)
初期化(initialization)
a:
0.7
1.1
2.1
2.7
2.4
2.2
1.8
Exercise
2. 上記の配列を作り，各要素の値を画面に出力せよ．
ヒント
double a[7];
a[2]=2.1;
14
14年10月22日水曜日
配列 (Arrays)
Initialization (もう一つのアプローチ)
double a[7] = {0.7,1.1,2.1,2.7,2.4,2.2,1.8};
a:
0.7
1.1
2.1
2.7
2.4
2.2
1.8
Exercise
3. 上記の配列を作り，各要素の値を画面に出力せよ．
15
14年10月22日水曜日
Exercise
4. 与えられた配列の要素の値を二乗にする
関数を書きなさい．
a:
0.7
1.1
2.1
2.7
2.4
2.2
1.8
a:
0.49 1.21 4.41 7.29 5.76 4.84 3.24
ヒント
#define NUMDATA 7
...
void sqr(double *arr){
int i;
for(i=0;i<NUMDATA;i++){
arr[i]= ...
}
}
14年10月22日水曜日
16
5. 自分が好きな関数の値を配列に入
れ，それらの値を画面に出せ，更に
gnuplotで確認しなさい．
x
a:
17
14年10月22日水曜日
6. 先ほど自分で選んだ関数の微分を画面に表示せよ．
ヒント
f (x + h)
f (x) = lim
h!0
h
0
f (x)
f (x + h)
⇡
h
f (x)
(0 < h ⌧ 1)
double f[100];
...
double df[99];
...
例
f (x) = sin(x) + cos(10x)/10
f (x)
14年10月22日水曜日
0
f (x)
18
ファイルからのデータを配列に入れてみよう
7. データを取得し，Cプログラムをセットアップする．
i.
http://mmc01.es.hokudai.ac.jp/~eginder/Data/sensor.txt
から
sensor.txtをダウンロードする．
ii. データの長さを求めよ．
iii.データを配列に保存する．
19
14年10月22日水曜日
Black Box
mydata.dat
20
14年10月22日水曜日
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define NUMDATA 7
void readdata(double *arr){
int i;
FILE *fp;
fp = fopen("mydata.dat","r");
if (fp == NULL) {
// report an error if unable to open the file.
printf("Unable to open mydata.dat for reading.\n");
exit(0); // requires the stdlib.h header file
}
for(i=0;i<NUMDATA;i++){
// READ DATA FROM FILE
fscanf(fp,"%lf",&arr[i]);
}
fclose(fp);
}
int main(void){
int i;
double a[NUMDATA];
readdata(a);
for(i=0;i<NUMDATA;i++){
// WRITE DATA TO SCREEN
printf("%lf\n",a[i]);
}
return 0;
}
14年10月22日水曜日
exit(0) . . .
データの数
ファイルを開く
r = “read”
数を一つずつ読み込む
ファイルを閉じる
配列を定義する
readdata()を呼ぶ
配列のデータを画面
に表示する．
21
ニュートン法 (Newton’s Method)
(for some ⇠k between xk and x⇤ )
テイラー展
0 = f (x⇤ ) = f (xk ) + (x⇤
開:
xk )f 0 (xk ) +
(x⇤
xk ) 2
2
f 00 (⇠k ),
従って，
x⇤ = xk
f (xk )
f 0 (xk )
(x⇤ xk )2 00
f (⇠k )
0
k
2(f (x ))
ニュートン法の近似
x = xk
f (xk )
.
0
f (xk )
xk+1 = xk
f (xk )
f 0 (xk )
⇤
繰り返す:
(k = 1, 2, ...)
22
14年10月22日水曜日
ニュートン法 (Newton’s Method)
f (x)
(x1 , f (x1 ))
x1
(x2 = x1
x
f (x1 )/f 0 (x1 ), 0)
23
14年10月22日水曜日
ニュートン法 (Newton’s Method)
zoom in
f (x)
(x2 , f (x2 ))
x2
(x3 = x2
(x1 , f (x1 ))
x1
x
0
f (x2 )/f (x2 ), 0)
24
14年10月22日水曜日
ニュートン法の問題点
• 必ずしも収束するとは言えない．
• 関数の微分を計算する事が不可欠である．
• 微分によって飛ぶ事がある．
例
f (x) = atan(x)
x0
x0
25
14年10月22日水曜日
Exercise
8. ニュートン法のCプログラムを書きなさい．
ただし，f (x)
= sin(x)とする．また，プログラムの
出力を以下のように表示せよ．
Remember:
xk+1 = xk
f (xk )
f 0 (xk )
26
14年10月22日水曜日
問題点
9. 下記の関数に対するニュートン法の性質を
調べなさい．
f (x) = x3
x
f (x) = x2
3
3, x0 = 0
(a cycle)
(slow convergence)
2
f (x) = x + x + x + 1 + cos(10x)
27
14年10月22日水曜日
10.*
newton
approximate
secant
Secant Method
f (xk )
f 0 (xk )
xk+1 = xk
f (xk )
f (xk ) ⇡
xk
0
xk+1 = xk
f (xk
xk 1
1)

xk
f (xk )
f (xk )
f (x)
f (xk
xk 1
f (xk
1)
1)
f (xk )
xk+1
14年10月22日水曜日
28

Download Report