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

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

ｘ²=C

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

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

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

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

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

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

プログラム

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

ニュートン法では、x0 としてなるべくＣの平方根に近い値を与えた方が良いので、まずＣの平方根より大きい最小の（つまりＣの平方根に最も近い）整数を探すようにしました。これは、単に変数を１づつ増やしながらその２乗が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; /* 求めた近似値を返す */

  }

プログラムソース表示

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