非線形方程式の数値解法〜ニュートン法〜荻田　武史目的 •  非線形方程式 f(x) = 0 の解を数値計算で求める。（例） f(x) = x2 -‐ 4 であれば、x=-‐2, 2のとき f(x) = 0 f(x)は、多項式から三角関数、指数関数などを含んで良い。 •  数値計算によって、厳密に f(x＊) = 0 を満たす真の解 x＊を求めることは難しいが、x＊に近い値（近似解）を求めることは可能。 •  本講義では、非線形方程式の反復解法であるニュートン法によって、f(x) = 0 の近似解を求める方法を学習し、そのアルゴリズムを理解した上で、プログラミングを行う。インライン関数を定義する方法＜例＞ f(x) = cos x -‐ x2 の場合 >> f=inline('cos(x) -‐ x^2') f = Inline func8on: f(x) = cos(x) -‐ x^2 >> f(0) ans = 1 >> ezplot(f,[-‐1,1]) ← f(x)のグラフを -‐1 ≦ x ≦ 1の範囲で見る。練習問題以下をインライン関数で定義して、グラフを描いてみよう（xの範囲も適当に変更してみよう）。 (1) f(x) = x2 - 2 (2) f(x) = x2 - 3x + 2 (3) f(x) = x - sin x ! 2π $ (4) f (x) = sin # x & " 1+ e % 反復解法について・f(x) = 0 について、x0を初期値として与えた後、決められた手順を繰り返して（反復）、x1 , x2 , … , xn と真の解x* へ近づける（収束）。・反復の停止条件： f(xn)が、ある程度0に近づいたら反復を止めることにする。 f (xn ) ≤ ε ε：十分に小さい正の数（たとえば、10-‐6など。プログラムでは、1e-‐6 のように書く）ニュートン法 •  連続で1回微分可能な関数 y = f(x) を考える。 •  ある点xnでの導関数f’(xn)が与えられるとき、 (xn,f(xn))を通り、傾きがf'(xn)の直線方程式は y = f '(xn )(x − xn ) + f (xn ) •  これをf(x)の近似解とみなすと、f(x) = 0の近似解は y = 0 のとき（x軸との交点）なので f (xn ) xn+1 = xn − f '(xn ) ニュートン法（つづき） •  このように、１次式で近似することを繰り返し、 •  以下を満たすまで行う。 f (xn ) ≤ ε x1, x2, x３と反復する毎に真の解x*に近づいていく x* 初期値【実習】ニュートン法のプログラム •  プログラム作成のポイント –  関数f(x)と導関数f'(x)をインライン関数で定義する。（f'(x)は手計算で求める） –  初期値x0を引数として与える。 –  for文で反復を行う。（収束しない場合もあるので、反復回数に上限を設ける） –  以下の反復条件を満たしたら、反復をbreakで止めて、解を出力する。 f (xn ) ≤ ε ニュートン法のプログラム func8on [x,k] = mynewton1(f,fd,x0) eps1 = 1e-‐6; % 許容誤差 kmax = 50; % 最大反復回数 x = x0; % xの初期値 for k=1:kmax fv = f(x); fg = fd(x); if abs(fv) < eps1 ← 反復停止条件 break; end if (abs(fg) <= eps * abs(fv)) ←割る数が0でないか（小さ過ぎないか）チェック break; end d = fv / fg; x = x -‐ d; ← xを更新 end 演習問題 •  以下の問題について、作成したニュートン法で解いてみよう。 (1) f(x) = x2 - 2 (2) f(x) = x2 - 3x + 2 (3) f(x) = x - sin x ! 2π $ (4)f (x) = sin # x & " 1+ e % 入力する初期値を色々変えてみて、収束速度を調べてみよう