5 素朴な Newton 法の実行例

complex-newton.C

   1 /*
   2  * complex-newton.C -- Newton 法で方程式 f(x)=0 を解く
   3  *  コンパイル: g++ -o complex-newton complex-newton.C
   4  *  実行: ./complex-newton
   5  */
   6 
   7 #include <iostream.h>
   8 #include <complex.h>
   9 
  10 int main()
  11 {
  12     int i, maxitr = 100;
  13     complex<double> f(complex<double>), dfdz(complex<double>), x, dx;
  14     double eps;
  15 
  16     cout << " 初期値 x0, 許容精度ε=";
  17     cin >> x >> eps;
  18 
  19     cout.precision(16);
  20     for (i = 0; i < maxitr; i++) {
  21         dx = - f(x) / dfdz(x);
  22         x += dx;
  23         cout << "f(" << x << ")=" << f(x) << endl;
  24         if (abs(dx) <= eps) break;
  25     }
  26     return 0;
  27 }
  28 
  29 /*
  30  * この関数の例は杉原・室田『数値計算法の数理』岩波書店, p.67 による
  31  *   f(z) = z^3-2z+2
  32  *   1実根(≒-1.769292354238631),2虚根(≒0.8846461771193157±0.5897428050222054)
  33  *   を持つが、原点の近くに初期値を取ると、0 と 1 の間を振動する。
  34  */
  35 
  36 complex<double> f(complex<double> z)
  37 {
  38   /* return z * z * z - 2 * z + 2; */
  39   return z * (z * z - 2) + 2;
  40 }
  41 
  42 /* 関数 f の導関数 (df/dx のつもりで名前をつけた) */
  43 complex<double> dfdz(complex<double> z)
  44 {
  45   return 3 * z * z - 2;
  46 }

complex-newton 実行結果

oyabun%	g++ -o complex-newton complex-newton.C
oyabun%	./complex-newton 
 初期値 x0, 許容精度ε=(0.07,0.07) 1e-15
f((1.000479890500462,0.01402105439035678))=(0.9998905285720316,0.01405867910352982)
f((0.008690071144917044,0.08327942929892811))=(1.982439704951537,-0.1671175729050391)
f((0.9899874424723547,0.002680506223014176))=(0.9902658531083521,0.002520260799472426)
f((-0.06300187281112724,0.01783194417259117))=(2.125813775994601,-0.03545722093219378)
(中略)
f((-9.992007221626409e-16,0))=(2.000000000000002,0)
f((0.9999999999999999,0))=(1,0)
f((-9.992007221626409e-16,0))=(2.000000000000002,0)
f((0.9999999999999999,0))=(1,0)
f((-9.992007221626409e-16,0))=(2.000000000000002,0)
f((0.9999999999999999,0))=(1,0)
f((-9.992007221626409e-16,0))=(2.000000000000002,0)
f((0.9999999999999999,0))=(1,0)
f((-9.992007221626409e-16,0))=(2.000000000000002,0)
f((0.9999999999999999,0))=(1,0)
f((-9.992007221626409e-16,0))=(2.000000000000002,0)
f((0.9999999999999999,0))=(1,0)
f((-9.992007221626409e-16,0))=(2.000000000000002,0)
f((0.9999999999999999,0))=(1,0)
f((-9.992007221626409e-16,0))=(2.000000000000002,0)
f((0.9999999999999999,0))=(1,0)
f((-9.992007221626409e-16,0))=(2.000000000000002,0)
oyabun%