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