ニュートン法で平方根を求める

ある数C の「平方根」というのは、二乗するとC になる数の事です。つまり、ある数をC 平方根をx とすると、

2=C

となる数であり、このx はx2-C=0の2次方程式を解けば求まる事になります。

ニュートン法による数値計算

今回は、このx2-C=0 という式を「接線」とx軸との交点を求めていくニュートン(ニュートン・ラプソン)法というアルゴリズムで解いてみましょう。

まず、解(正)よりも右側にある適当なx 座標(x0)におけるf(x)=x2-Cの接線とx 軸(y=0)との交点を求めます。f'(x)=2x なので、この接線の傾きは2x0であり式は 2x0+α(平方根を求めるからα<0)です。

次に、傾きというのはy の変化をx の変化で割ったものなので、交点のx 座標をx1 とすると2x0=f(x0)/x0-x1 となります。この式に、f(x0)=x02-C を代入して整理すると、x1=(x0+C/x0)/2 となり、交点(x1,0)を計算することが出来ます。
さて、この交点x1 はx0 とf(x)とx 軸との交点の間にある、つまりx0 よりもより求める解(平方根)に近い事がわかるでしょう。なぜなら、今回の場合接線の傾きはx が大きくなるほど大きくなりxが小さくなるほど小さくなる上、常に正だからです。

という事は、同じようにしてx1 からx2 を、x2 からx3を、というふうにして交点を求めていけば、そのx 座標の値は、次第(実際は急激)にCの平方根に近づくわけです。

プログラム

では、ニュートン法のアルゴリズムをプログラムにしてみましょう。

ニュートン法では、x0 としてなるべくCの平方根に近い値を与えた方が良いので、まずCの平方根より大きい最小の(つまりCの平方根に最も近い)整数を探すようにしました。これは、単に変数を1づつ増やしながらその2乗がCより大きくなるのを待つだけで実現できますね。
後は、こうして求めたx0 を使って上とまったく同じ計算を何度かやるだけです。

・引数の平方根を返す関数
  double sqr(double c) { /* c の平方根を返す*/

      double xn;
      int i,n;

      n=10;
      xn=0;

      do { /* 平方根に最も近い大きい方の整数を探す */

          xn+=1;

      } while (xn*xn<c);

      for (i=0;i<n;i++) {

          xn=(xn+c/xn)/2; /* 交点のx 座標を求める */

      }

      return xn; /* 求めた近似値を返す */

  }

プログラムソース表示

例として、入力した値の平方根をニュートン法で求めるプログラムを作っておきました。収束の様子も見れるようにしたので、いろいろな値(正の整数)を入れて値の変化を見てみてください。10回計算を繰り返し、その値をグラフに表示していますが、かなり速く収束しています。