方程式の反復解法 (2)

C プログラミング
第４回
— 方程式の反復解法 (2) —
基幹理工学部
1/12
復習
ファイル名：041.c
ニュートン法を用いて方程式 f (x) = cos x − x2 = 0 の
0 < x < 1 の範囲にある解の近似値（小数点以下６位まで）
を求めよ．初期値 x(0) は 1.0 とし，表記は以下のようにせよ．
initial value: 1
x(001) = ...
x(002) = ...
...
answer = ...
▶
▶
▶
右は f (x) の図．
f ′ (x) = − sin x − 2x．
前回作ったプログラム
のどこを直せばよいか
を考えよ．
2/12
ニュートン法のプログラム
▶
方程式 f (x) = 0 に対して，二つの関数
double func(double x)
→ 関数 f (x)
double func_d(double x) → 導関数 f ′ (x)
を定義している．
▶
▶
▶
⇓
関数 f (x) が複雑な場合，導関数 f ′ (x) の計算が面倒に
なる．
導関数の記述を間違える可能性がある．
⇓
数値微分を行うことにより，f (x) から自動的に微分値
f ′ (x) の近似を計算する．
3/12
ニュートン法のプログラム
▶
方程式 f (x) = 0 に対して，二つの関数
double func(double x)
→ 関数 f (x)
double func_d(double x) → 導関数 f ′ (x)
を定義している．
▶
▶
▶
⇓
関数 f (x) が複雑な場合，導関数 f ′ (x) の計算が面倒に
なる．
導関数の記述を間違える可能性がある．
⇓
数値微分を行うことにより，f (x) から自動的に微分値
f ′ (x) の近似を計算する．
3/12
ニュートン法のプログラム
▶
方程式 f (x) = 0 に対して，二つの関数
double func(double x)
→ 関数 f (x)
double func_d(double x) → 導関数 f ′ (x)
を定義している．
▶
▶
▶
⇓
関数 f (x) が複雑な場合，導関数 f ′ (x) の計算が面倒に
なる．
導関数の記述を間違える可能性がある．
⇓
数値微分を行うことにより，f (x) から自動的に微分値
f ′ (x) の近似を計算する．
3/12
数値微分
１変数関数 f (x) に対しての微分の定義は
f (x + h) − f (x)
h→0
h
f ′ (x) = lim
である．
この代わりに，有限の小さい正数 h を用いて
f ′ (x) ≈
f (x + h) − f (x)
h
により微分値の近似を計算する方法を数値微分という．
4/12
数値微分
▶
数値微分を行う際の刻み幅 h の値を適当な微小値（例
えば h = 1e − 4）として定義することで，ニュートン法
のプログラムは例えば，以下のように変更できる．
int main(void){
…
i=0;
while (i<max){
f = func(x);
df = func_d(x);
x_n= x- f/df;
…
}
…
return 0;
}
int main(void){
…
h=1e-4;
i=0;
while (i<max){
f = func(x);
df = (func(x+h)-func(x))/h;
x_n= x- f/df;
…
}
…
return 0;
⇒
}
5/12
演習
ファイル名：042.c
041.c を修正して f ′ (x) の箇所を数値微分を用いたものに書
き変えよ．初期値 x(0) = 1.0 とし，表記は以下のようにせよ．
initial value: 1
x(001) = ...
x(002) = ...
...
answer = ...
▶
右の図は，
f (x) = cos x − x2 のグ
ラフである．
6/12
前項で作成したプログラムでは…
double func(double x){
return x - cos(x);
}
double NewtonRaphson(double x0, double e, int max){
…
i=0;
while(i<max){
f = func(x);
df = (func(x+h)-func(x))/h;
…
}
▶
▶
方程式を表す関数 func( ) をそのままニュートン法を行
う関数 newton( ) の中で利用している．
⇓
方程式を表す関数が別の名前（例えば，func2( )) で与
えられたとき，関数 newton( ) の中身も修正する必要が
ある．
7/12
前項で作成したプログラムでは…
double func(double x){
return x - cos(x);
}
double NewtonRaphson(double x0, double e, int max){
…
i=0;
while(i<max){
f = func(x);
df = (func(x+h)-func(x))/h;
…
}
▶
▶
方程式を表す関数 func( ) をそのままニュートン法を行
う関数 newton( ) の中で利用している．
⇓
方程式を表す関数が別の名前（例えば，func2( )) で与
えられたとき，関数 newton( ) の中身も修正する必要が
ある．
7/12
前項で作成したプログラムでは…
double NewtonRaphson(double x0, double e, int max){
…
i=0;
while(i<max){
f = func2(x);
df = (func2(x+h)-func2(x))/h;
…
}
double func2(double x){
return x - cos(x);
}
⇓
関数へのポインタを利用することで．方程式を表す関数名
を newton( ) への引数として渡すことができる．
例： newton (func2,1.0,1e-6,100);
8/12
関数へのポインタ
▶
▶
▶
▶
変数や配列と同様に，関数もメモリ上のどこかに配置
される．つまり，変数や配列と同様，関数にもアドレス
が与えられる．
配列名は，プログラム中では配列の先頭要素へのポイ
ンタに読み替えられていた．
同様に，関数名はプログラム中では関数へのポイン
タへと読み替えられる．
例えば，関数
double function(double x){
…
}
が定義されているとき，関数名 function は関数
function( ) の先頭アドレスを指す．
9/12
関数へのポインタ（プログラム例）
#include <stdio.h>
double add(double x, double y){
return x+y;
}
double sub(double x, double y){
return x-y;
}
int main(void){
double (*func_p)(double,double);
/* func_p は２つの double 型の引数 */
/* 関数へのポインタとして定義． */
func_p = add;
/* func_p に関数 add へのポインタを代入 */
printf("%fY
=n",func_p(3.0,2.0));
func_p = sub;
/* func_p に関数 sub へのポインタを代入 */
printf("%fY
=n",func_p(3.0,2.0)); return 0;
}
10/12
演習
ファイル名: 043.c
単位円 x2 + y 2 = 1 と直線 y = x + 1/2 の２つの交点の座標の近似解
（小数点以下第６位まで）を求めるプログラムを次の事実を元に作りな
√
さい: 単位円 x2 + y 2 = 1 は y = ± 1 − x2 という２つの関数．点の x
座標は，以下の２つの関数の零点である．
√
√
f (x) = x + 1/2 − 1 − x2 ,
g(x) = x + 1/2 + 1 − x2
2
2
x +y =1
1.5
1
y
0.5
0
-0.5
-1
-1.5
-1.5
-1
-0.5
0
x
0.5
1
1.5
11/12
演習
ファイル名：043.c
▶
▶
▶
▶
前頁のことを踏まえて，単位円 x2 + y 2 = 1 と直線
y = x + 1/2 の交点の近似解を求めるプログラムを作れ．
数値微分，関数へのポインタは必ず使うこと．
ニュートン・ラプソン法を用いる際の ε, IM AX はそれぞ
れ ε = 10−6 , IM AX = 20 とせよ．初期値は f (x) に対し
ては x0 = 0，g(x) に対しては x0 = −1 とせよ．
結果の表示は以下のようにせよ．
Newton-Raphson method:
Upper part: ( ?????, ????? )
Lower part: ( ?????, ????? )
12/12

Download Report