方程式の数値解析

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

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

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

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

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. y y f(b) f(c) f(c) a a x x b b f(a) スタート 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 = -2.638 三角関数の例 X = 1.114 解

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

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

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

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

29. z3 + 6 z2 + 21 z +32 =0 実数部 虚数部 f(x, y)= g(x, y)= z = x + yi とおく =(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 x3-3xy2+6x2-y2+21x+32=0 3x2y-y3+12xy+21y=0

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

31. f(x,y) x y 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 f(x0,y0) (a,b) f(x,y)=0 g(x,y)

32. x ← x1 y ← y1 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 dy dy x0 ← x1 y0 ← y1 {f(x1,y1)}2+{g(x1,y1)}2<e2 No Yes

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

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

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

36. double _Complex x0, x1; double eps = 1.0e-6; int i, n = 10000; x0 = -1.0 + 2.0*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; } 初期値 複素数の絶対値 複素数の実部 複素数の虚部