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 // Π (ｚi-ｚj) の計算
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 }