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

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

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

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

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

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

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

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

メインプログラム 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); 初期値 変数の定義 説明文 初期値の出力 二分割の関数プログラムへ 結果の出力

二分割関数プログラム 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); } 関数名，引数 変数の定義 繰り返し計算 中点の計算 判定式 終了の判定 結果の出力

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

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

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

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

収束判定の改善 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); } 関数名，引数 変数の定義 繰り返し計算 中点の計算 判定式 終了の判定 結果の出力

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

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

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

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

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

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

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); 初期値 変数の定義 （）内が条件を満足しないか，０で停止 次の値 終了の判定 結果の出力

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); 関数 導関数

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

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

別の収束判定法 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

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

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

z = x + yi とおく z3 + 6 z2 + 21 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

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)

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)

{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

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

#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) 複素数のため 複素数の定義 導関数

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

if(cabs(f(x1)) < eps){break;} x0 = x1; } 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; 初期値 複素数の絶対値 複素数の実部 複素数の虚部