next up previous contents
Next: B.6.2 コンパイル&実行例 Up: B.6 DK 法のサンプル・プログラム Previous: B.6 DK 法のサンプル・プログラム

B.6.1 ソースプログラム DK3.C


   1 // DK3.C -- DKA 法で代数方程式の根を求める
   2 //
   3 //   数学科 6701 号室のワークステーション環境でのコンパイル
   4 //     (グラフィックス・ライブラリィ GLSC を使っているので
   5 //      情報科学センターではコンパイルできない。悪しからず。)
   6 //     g++ -o DK3 DK3.C -I/usr/local/include -lglscd -lX11 -lsocket
   7 //   あるいは
   8 //     ccmg+ DK3.C
   9 //
  10 //  Win32 環境下の glscwin でも (若干の修正の下に) 動作させられることが
  11 //  分かったので、画面に表示するメッセージを変更しました (2000/10/18)。
  12 
  13 #include <iostream.h>
  14 #include <iomanip.h>   // setprecision() のため
  15 #include <math.h>
  16 #include <complex.h>   // complex<double> のため
  17 
  18 // GLSC のヘッダーファイルをインクルード
  19 extern "C" {
  20 #define G_DOUBLE
  21 #include <glsc.h>
  22 }
  23 
  24 // プロトタイプ宣言
  25 complex<double> polynomial(int, complex<double> *, complex<double>);
  26 complex<double> bunbo(int, complex<double> *, int);
  27 
  28 // 解きたい代数方程式の次数の最高値
  29 #define MAXN 100
  30 
  31 int main()
  32 {
  33   int i, k, n;
  34   complex<double> a[MAXN+1], x[MAXN], newx[MAXN], dx[MAXN];
  35   complex<double> g, I(0,1);
  36   double r0, R, max, pi;
  37 
  38   // 数学定数 (円周率) の準備
  39   pi = 4 * atan(1.0);
  40 
  41   // 表\示の桁数を 16 桁にする
  42   cout << setprecision(16);
  43 
  44   // 方程式の入力
  45   cout << "次数 n を入力してください (1≦n≦" << MAXN << "): ";
  46   cin >> n;
  47   if (n > MAXN || n <= 0) {
  48     cerr << "次数は" << MAXN << "以下の自然数として下さい。" << endl;
  49     exit(0);
  50   }
  51   for (i = 0; i <= n; i++) {
  52     cout << (n-i) << "次の係数を入力してください: ";
  53     cin >> a[i];
  54   }
  55 
  56   // 多項式を最高次の係数で割って monic にする
  57   cout << "monic にします。" << endl;
  58   for (i = 1; i <= n; i++)
  59     a[i] /= a[0];
  60   a[0] = 1;
  61   cout << "修正した係数" << endl;
  62   for (i = 0; i <= n; i++)
  63     cout << "a[" << i << "]=" << a[i] << endl;
  64 
  65   // Aberth の初期値を配置する円の決定
  66   g = - a[1] / n;
  67   cout << "根の重心" << g << endl;
  68   max = 0;
  69   for (i = 1; i <= n; i++)
  70     if (abs(a[i]) > max)
  71       max = abs(a[i]);
  72   cout << "max|a_i|=" << max << endl;
  73   r0 = abs(g) + 1 + max;
  74   cout << "根は重心中心で、半径 r0=" << r0 << "の円盤内にある" << endl;
  75 
  76   cout << "円の半径 (分からなければ上の値を指定してください): ";
  77   R = r0;
  78   cin >> r0;
  79   if (r0 > R)
  80     R = r0;
  81   cout << "図は根の重心を中心として半径 " << R
  82 #ifdef __CYGWIN__
  83        << "の円が表\示できるようにします。" << endl;
  84 #else
  85        << "の円が表示できるようにします。" << endl;
  86 #endif
  87 
  88   // Aberth の初期値
  89   cout << "初期値" << endl;
  90   for (i = 0; i < n; i++) {
  91     double theta;
  92     theta = 2 * i * pi / n + pi / (2 * n);
  93     x[i] = g + r0 * exp(I * theta);
  94     cout << x[i] << endl;
  95   }
  96 
  97   // グラフィックス・ライブラリィ GLSC の初期化
  98   g_init("DKGRAPH", 120.0, 120.0);
  99   g_device(G_BOTH);
 100   // 座標系の指定
 101   g_def_scale(0,
 102 	      real(g) - 1.1 * R, real(g) + 1.1 * R,
 103 	      imag(g) - 1.1 * R, imag(g) + 1.1 * R,
 104 	      10.0, 10.0, 100.0, 100.0);
 105   g_sel_scale(0);
 106   // 線種、マーカーの指定
 107   g_def_line(0, G_BLACK, 0, G_LINE_SOLID);
 108   g_sel_line(0);
 109   g_def_marker(0, G_RED, 2, G_MARKER_CIRC);
 110   g_sel_marker(0);
 111 
 112   // 初期値と初期値を置く円を描く
 113   g_circle(real(g), imag(g), r0, G_YES, G_NO);
 114   for (i = 0; i < n; i++)
 115     g_marker(real(x[i]), imag(x[i]));
 116 
 117   // 反復
 118   for (k = 1; k <= 1000; k++) {
 119     double error;
 120     cout << "第" << k << "反復" << endl;
 121     for (i = 0; i < n; i++) {
 122       dx[i] = polynomial(n, a, x[i]) / bunbo(n, x, i);
 123       newx[i] = x[i] - dx[i];
 124       // 図示する
 125       g_move(real(x[i]), imag(x[i]));
 126       g_plot(real(newx[i]), imag(newx[i]));
 127       g_marker(real(newx[i]), imag(newx[i]));
 128     }
 129     // 更新
 130     for (i = 0; i < n; i++) {
 131       x[i] = newx[i];
 132       cout << x[i] << endl;
 133     }
 134     // 変化量を計算する (ここは非常に素朴)
 135     error = 0.0;
 136     for (i = 0; i < n; i++)
 137       error += abs(dx[i]);
 138     cout << "変化量=" << error << endl;
 139     // 変化量が小さければ収束したと判断して反復を終了する。
 140     if (error < 1e-12)
 141       break;
 142   }
 143 
 144   // GLSC 終了の処理
 145   cout << "終了には" << endl;
 146   cout << " グラフィックスのウィンドウをクリック (X の場合)。" << endl;
 147   cout << " グラフィックスのウィンドウをクローズ (Windows の場合)。" << endl;
 148   g_sleep(-1.0);
 149   g_term();
 150 
 151   return 0;
 152 }
 153 
 154 // 多項式の値の計算 (Horner 法)
 155 complex<double> polynomial(int n,
 156 			   complex<double> *a,
 157 			   complex<double> z)
 158 {
 159   int i;
 160   complex<double> w;
 161   w = a[0];
 162   for (i = 1; i <= n; i++)
 163     w = (w * z + a[i]);
 164   return w;
 165 }
 166 
 167 //  Π (zi-zj) の計算
 168 // j≠i
 169 complex<double> bunbo(int n,
 170 		      complex<double> *z,
 171 		      int i)
 172 {
 173   int j;
 174   complex<double> w(1,0);
 175   for (j = 0; j < n; j++)
 176     if (j != i)
 177       w *= (z[i] - z[j]);
 178   return w;
 179 }


next up previous contents
Next: B.6.2 コンパイル&実行例 Up: B.6 DK 法のサンプル・プログラム Previous: B.6 DK 法のサンプル・プログラム
Masashi Katsurada
平成21年7月9日