B.5 素朴な Newton 法の実行例

complex-newton.C

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

complex-newton 実行結果

% g++ complex-newton.cpp
% ./a.out 
 初期値 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((1.00568580973205,-0.003615972418678629))=(1.005743530042432,-0.00373963422208404)
(中略)
f((0.9999999999999999,1.517769664024309e-319))=(1,1.517769664024309e-319)
f((-9.992007221626409e-16,9.106617984145856e-319))=(2.000000000000002,-1.821323596829171e-318)
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)
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)
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)
%