プログラミング III 資料 4 方程式の解法;ニュートン法 ◎ 4 ニュートン法を使って 2 1.41421356 を誤差 10 の精度で求めることにする。 つまり f ( x) x 2 2 0 のプラス解を求める。 ※この問題は前回 2 分法で解いている。計算回数と真値に近づく収束の早さにも注意する 今日の内容と達成目標 ・ 初期値(初期近似解)を与える 「ニュートン法アルゴリズムの理解」 「接線との交点が解に近づく」方法を幾何学的に理解する ◎ 2 分法は簡単なアルゴリズムで方程式を解ける代わりに求め方が粗く、計算回数がかかる。それは収束速度が 1 次 であるからである。ここでは、一般的に「方程式の数値解法」といえばどのテキストでも扱う「ニュートン法(ニュ ートン・ラフソン法) 」について触れる。この方法は 2 次で収束するので計算回数も少なくてすむ。 数式的には理解の難易度が多少上がるので、どちらかというと、図を見て幾何学的に理解して欲しい。 まずは x0 =2.0 から計算を始めることにする ------------------------------ 1 回目の計算 ----------------------------------------- ・ 初期近似値 x0 =2.0 における f( x0 )の接線と、x 軸との交点、次の 2 f(x)=x –2 近似解 x1 を計算する(図 1) ここで、関数 f (x ) の微分関数 f ' ( x ) が必要になる f ' ( x) 2 x ニュートン法(ニュートン‐ラフソン法) f(x) 関数 f (x ) がある時、 f ( x) 0 なる方程式の解を求める。 図 1 に示すように、 まず初期値 x0 において f(x0)の接線を求める。 わる。この交点を x1 とすると、x1 は x0 よりも解(つまり、x 軸 との交点)に近づいていることがわかる。 x1 さらに x1 における f(x1)の接線を求めると、今度は x2 と交わる。 x x0 ◎これらを数式で考える との交点を x1 とする。 ・関数 f(x0)の接線の傾きは微分係数 f ' ( x0 ) で与えられる。 f(x) x1 周辺の拡大図 交点のxがf(x)=0の解 x x2 x1 xは解に近づく 図 2 x1 から接線を引き、x軸との交点 x2 ◎ニュートン法における収束条件 f (x ) ことである。また、上記の例は計算区間で下に凸であり、 xk 1 が いは「2 乗収束する」という。 ある次の近似値が真値 error x1 x0 0.5 ゆえに、まだ計算は続行する 2 回目の計算 ------------------------------------- f (1.5) 1.5 1.5 2.0 0.25 f ' (1.5) 2.0 1.5 3.0 ( f ( x) 0 なる x)に近づ いていることがわかる。 x2 x1 f ( x1 ) 0.25 1.5 1.41666 f ' ( x1 ) 3 解周辺の拡大図 ゆえに、まだ計算は続行する。ただ、この時点で 1.41 とすでに小数 3 回目の計算 ----------------------------------------- 1.40 1.42 1.44 1.46 1.48 1.50 x 図 2. 2 回目の近似値計算。 ・ x2 =1.4166 における f( x2 )の接線との交点、次の近似解 x3 を計算する 2 回目の計算ですでに関数 f ( x2 ) 1.414216 f ' ( x2 ) と接線がほとんど一致して f (1.4166) 0.0069425 f ' (1.4166) 2.83333 x3 x2 ・収束(精度)の判定をする error x3 x2 0.00245 ゆえに、まだ計算は続行する。 いる。幾何学的にはこれから 先はかなり拡大してもほぼ 一致しているように見えて しまうレベルである に四捨五入して小数以下 6 桁で表示するので注意! ------------------------- 以下、繰り返し ----------------------------------------- 関数 f(x) は方程式 f(x)=0 の解 xt を挟んで この後の 4 回目の計算ですでに 1.414213562374690 であり、真値の 2 と小数以下 11 桁目まで一致している。2 必ずプラスとマイナスの値をとる(中間値の 分法では 4 回目の計算ではまだ誤差が 0.125 のレベルであり、小数以下 11 桁まで一致するには 41 回の計算が必要と 定理)ので、区間(x1,x2)に方程式の解が存在 なる。ニュートン法は 2 分法の約 1/10 の計算で済んだことになる。ニュートン法の、解に収束する速度(つまり 2 次 するならば の収束)が早いことが理解してもらえたと思う。 f ( x1) f ( x 2) 0 ニュートン法のような計算法は一般に「2 次の収束速度」ある (無限ループにして、条件を満たせば break でも良い) printf で途中経過をチェックするときは%15.15lf じゃないと勝手 ニュートン法における絶対的な収束条件は f ' ( x ) ≠0 である りして収束しないことがわかるだろう。 接線からの傾きとの交点で ------------------------------ とする。さらに解に近くなっている 図 3 のような関数では次の近似解が解を挟んで行ったり来た 図 1. 1 回目の近似値計算。 以下 2 桁まで一致していることに注意して欲しい f ( xk ) xk 1 xk f ' ( xk ) よっては収束しない場合もある。 x error x2 x1 0.08333 一般に、ニュートン法における xk の次の近似値 xk+1 は xk よりも小さくなるので収束する。しかし関数の形と初期値に 2.0 ・収束(精度)の判定をする f ' ( x0 ) で表される。今の値を使って次の値を計算していく。 1.5 ・ x1 =1.5 における f( x1 )の接線と x 軸との交点、次の近似解 x2 を計算する 図 1 初期値 x0 において接線を引き、x軸 x0 x1 f ( x0 ) 2 2.0 1.5 f ' ( x0 ) 4 ニュートン法では、前の値と今計算した値の誤差や桁の一致を検討する ------------------------------ この接線の傾きは x0における微分係数 この x2 はさらに真値に近づいていることがわかる。 ・ ゆえに、 f ' ( x ) f ( x0 ) より x x f ( x0 ) となる。 1 0 0 x1 x0 ・収束(精度)の判定をする この接線は傾きが 0 でなければ必ずどこか 1 点で必ず x 軸と交 f ( x0 ) ・ 一方で傾きは とも表される x0 x1 f ( 2.0) 2.0 * 2.0 2.0 2.0 f ' ( 2.0) 2.0 * 2.0 4.0 となる(ただし、x1,x2 のどちらかが解だった ※プログラム技術としては、ニュートン法は①次の近似値を計算して②計算条件を満たしているか判定 の 2 段階で ならば f ( x1) 図f (3x 2収束しない例 ) 0 となる)。 済む。中間値を計算したり、判定して範囲を入れ替えたりなどの面倒な操作が必要ないので 2 分法よりも簡単である。 (林英輔ら著「数値計算」森北出版 より) x1 x 2 なお、区間の中間値は で与えられる。 2
© Copyright 2024