解析的には解が得られない 方程式を数値的に求める。 例：３次方程式

Presentation on theme: "解析的には解が得られない 方程式を数値的に求める。 例：３次方程式"— Presentation transcript:

1 解析的には解が得られない 方程式を数値的に求める。 例：３次方程式
方程式の数値解析 解析的には解が得られない 方程式を数値的に求める。 例：３次方程式 f(x) = x3 + x = 0

2 電卓を用いて f(x) = x3 + x - 1 = 0 f(0) = - 1.0 f(1) = 1.0 f(0.5) = - 0.375
0.75 0.5 1

3 右にいくか，左に行くかはどうやって判断するか？
電卓を用いて f(x) = x3 + x - 1 = 0 右にいくか，左に行くかはどうやって判断するか？ ある幅に達したとき終了 0.5 1

4 ２分法 関数が負と正となる間に必ず解が有る。 f(a)* f(b)<0 x y a b f(a) f(b)

5 中点の関数値を求める。 f(a)*f(c)<0 a-c 間に解 y f(c) f(b) a x x c b f(a)

6 中点cを新たに端点bとする。 f(a)*f(b)<0 a-b 間に解 y f(b) a x x c b f(a)

7 x y a b f(a) f(c) x y a b f(b) f(c) f(c) f(a)*f(c)<0 b ← c a ← c
スタート f(c) Yes No f(a)*f(c)<0 b ← c a ← c Yes No |b-a|<e x = c

8 テキストのプログラム #include <stdio.h> #include <math.h>
#define EPS .0001 double nibun(double,double); double func_y(double); int main() { return 0; } Ｃプログラムの開始 数学関数を利用 数値の設定 関数定義（内容は プログラム後に記載） メインプログラム

9 メインプログラム double a=0.0, b=1.0; double x; printf("x^3 + x -1\n");
printf("%6.3lf\n",a); printf("%6.3lf\n",b); x = nibun(a,b); printf("%6.3lf\n",x); 初期値 変数の定義 説明文 初期値の出力 二分割の関数プログラムへ 結果の出力

10 二分割関数プログラム double nibun(double a, double b) { double c; do{
c=(a+b)/2.0; if((func_y(c)*func_y(a))<0) b=c; else a=c; }while(fabs(a-b)>EPS); return(c); } 関数名，引数 変数の定義 繰り返し計算 中点の計算 判定式 終了の判定 結果の出力

11 関数の定義 例：f(x) = x3 + x - 1 double func_y(double x) {
return(pow(x,3)+x-1.0); }

12 例題 関数： f(x) = x3 + 6x2 + 21x +32 =0 初期値： a = -5 b = -1
関数： f(x) = x sin x -1 =0 初期値： a = 1 b = 2

13 解 ３次関数の例 X = 三角関数の例 X = 1.114

14 収束判定について c = 0.0012345のとき fabs(b-a) < 0.0001 a b c = 0.0123???
ある幅に達したとき終了 c = のとき fabs(b-a) < a b c = ??? fabs((b-a)/c) < c = ??

15 収束判定の改善 double nibun(double a, double b) { double c; do{ c=(a+b)/2.0;
if((func_y(c)*func_y(a))<0) b=c; else a=c; }while(fabs((a-b)/c)>EPS); return(c); } 関数名，引数 変数の定義 繰り返し計算 中点の計算 判定式 終了の判定 結果の出力

16 二分法の注意点 二つの初期値に対して， 関数が正負反転していること．

17 ニュートン法 関数の傾き（微係数）を用いて求析。 y = f’(x0)(x-x0) + f(x0) y f(x0) x x0
x1 = x0 - f(x0)/f’(x0)

18 ニュートン法 y x x3 x2 x1 x0

19 xi+1 = xi - f(xi)/f’(xi) xi ← xi+1 |xi+1-xi|<e x = xi+1 スタート x0 y x
No x0 y x x2 x1 Yes x = xi+1

20 多項式の例 y(x) = x3 + x‐1 =0 の解をニュートン法で求めよ。

21 テキストのプログラム #include <stdio.h> #include <math.h>
#define EPS .0001 double func_y(double); double func_z(double); int main() { return 0; } Ｃプログラムの開始 数学関数を利用 数値の設定 関数定義 導関数の定義 メインプログラム

22 c = a - func_y(a)/func_z(a); if(fabs(a-c)>=EPS) a=c; else break; }
double a=1.0; double c; while(1) { c = a - func_y(a)/func_z(a); if(fabs(a-c)>=EPS) a=c; else break; } printf("近似式x=%6.3lf\n",c); 初期値 変数の定義 （）内が条件を満足しないか，０で停止 次の値 終了の判定 結果の出力

23 double func_y(double x) { return(pow(x,3.0) + x - 1.0); }
double func_z(double x) return(3.0*pow(x,2.0) + 1.0); 関数 導関数

24 章末問題 y = x ln x – 1 = 0 の解をニュートン法で求めよ。 ln x → log(x)

25 double func_y(double x) { return(x*log(x) - 1.0); }
double func_z(double x) return(log(x) + 1.0); 関数 導関数

26 別の収束判定法 xi+1 = xi - f(xi)/f’(xi) xi ← xi+1 |f(xi+1)|<e x = xi+1
スタート xi+1 = xi - f(xi)/f’(xi) xi ← xi+1 |f(xi+1)|<e No x0 y x x2 x1 Yes x = xi+1

27 ニュートン法の注意点 初期値によっては，予想外の解が出る． y a0 x x0

28 の複素解をニュートン・ラフソン法で求める。
複素数の解を求める f(z) = z3 + 6 z z +32 = 0 の複素解をニュートン・ラフソン法で求める。 二次元ニュートン・ラフソン法

29 z = x + yi とおく z3 + 6 z z +32 =0 =(x+yi)3+6(x+yi)2+21(x+yi)+32 (x3+3x2yi-3xy2-y3i)+6(x2-2xyi-y2)+21(x+yi)+32=0 (x3-3xy2+6x2-y2+21x+32)+(3x2y-y3+12xy+21y)i=0 実数部 虚数部 f(x, y)= g(x, y)= x3-3xy2+6x2-y2+21x+32=0 3x2y-y3+12xy+21y=0

30 f(x,y) = f(x0,y0) + ―――――― (x-x0) + ―――――― (y-y0) dx dy = 0
df(x0,y0) df(x0,y0) f(x,y) = f(x0,y0) + ―――――― (x-x0) + ―――――― (y-y0) dx dy = 0 x y f(x,y) f(x0,y0) f(x,y)=0 f(x,y)

31 f(x,y) = f(x0,y0) + ―――――― (x-x0) + ―――――― (y-y0) = 0 dx dy
df(x0,y0) df(x0,y0) f(x,y) = f(x0,y0) + ―――――― (x-x0) + ―――――― (y-y0) = 0 dx dy dg(x0,y0) dg(x0,y0) g(x,y) = g(x0,y0) + ―――――― (x-x0) + ―――――― (y-y0) = 0 dx dy x y f(x,y) f(x0,y0) (a,b) f(x,y)=0 g(x,y)

32 {f(x1,y1)}2+{g(x1,y1)}2<e2
df(x0,y0) dg(x0,y0) df(x0,y0) dg(x0,y0) J0 = ―――― ―――― － ―――― ―――― dx dy dy dx スタート dg(x0,y0) df(x0,y0) x1 = x0 - {f(x0,y0) ――――― － g(x0,y0) ―――――}/J0 dy dy df(x0,y0) dg(x0,y0) y1 = y0 - {g(x0,y0) ――――― － f(x0,y0) ―――――}/J0 x0 ← x1 y0 ← y1 {f(x1,y1)}2+{g(x1,y1)}2<e2 No Yes x ← x1 y ← y1

33 の複素解をニュートン・ラフソン法で求める。
複素数の解を求める y = x x x +32 の複素解をニュートン・ラフソン法で求める。 複素変数を用いて！

34 #include <stdio.h> #include <math.h>
#include <complex.h> double _Complex f(double _Complex x) { return ((x +6.0)*x )*x +32.0; } double _Complex fd(double _Complex xd) return (3.0*xd )*xd +21.0; int main(void) 複素数のため 複素数の定義 導関数

35 double _Complex x0,x1; x0 = *I; cabs(x0) creal(x0)→実数部 cimag(x0)→虚数部

36 if(cabs(f(x1)) < eps){break;} x0 = x1; }
double _Complex x0, x1; double eps = 1.0e-6; int i, n = 10000; x0 = *I; for(i = 0; i < n; i ++){ x1 = x0 - f(x0)/fd(x0); if(cabs(f(x1)) < eps){break;} x0 = x1; } printf(“x=%f+%fi\n”,creal(x1),cimag(x1)); return 0; 初期値 複素数の絶対値 複素数の実部 複素数の虚部