eigenbox代码
原eigenbox的代码如下：
int init(bool restarting) {
solver->add(f, "f");
solver->add(g, "g");
return 0;
}
int rhs(BoutReal t) {
mesh->communicate(f);
ddt(g) = D2DX2(f);
ddt(f) = g;
return 0;
}
可以看到其核心方程为
$$\pdv[2]{f}{t}=\pdv[2]{f}{x}$$
增加打印矩阵
原矩阵为shellMat类型，不能直接打印，需要转化为显式形式才能打印出来，使用如下代码进行转化：
// +++ 求解后检查矩阵状态 +++
// added by myron
output.write("myron=> Debug: Post-solve matrix state");
Mat explicit_mat;
MatConvert(shellMat, MATAIJ, MAT_INITIAL_MATRIX, &explicit_mat);
MatView(explicit_mat, PETSC_VIEWER_STDOUT_WORLD);
MatDestroy(&explicit_mat);打印的矩阵内容如下:
row 0: (0, 0) (1, 1) (2, 0) (3, 0) (4, 0) (5, 0) (6, 0) (7, 0) (8, 0) (9, 0) (10, 0) (11, 0) (12, 0) (13, 0) (14, 0) (15, 0) (16, 0) (17, 0) (18, 0) (19, 0) (20, 0) (21, 0) (22, 0) (23, 0) (24, 0) (25, 0) (26, 0) (27, 0) (28, 0) (29, 0) (30, 0) (31, 0) (32, 0) (33, 0) (34, 0) (35, 0) (36, 0) (37, 0) (38, 0) (39, 0) (40, 0) (41, 0) (42, 0) (43, 0) (44, 0) (45, 0) (46, 0) (47, 0) (48, 0) (49, 0) (50, 0) (51, 0) (52, 0) (53, 0) (54, 0) (55, 0) (56, 0) (57, 0) (58, 0) (59, 0) (60, 0) (61, 0) (62, 0) (63, 0)
row 1: (0, -3072) (1, 0) (2, 1024) (3, 0) (4, 0) (5, 0) (6, 0) (7, 0) (8, 0) (9, 0) (10, 0) (11, 0) (12, 0) (13, 0) (14, 0) (15, 0) (16, 0) (17, 0) (18, 0) (19, 0) (20, 0) (21, 0) (22, 0) (23, 0) (24, 0) (25, 0) (26, 0) (27, 0) (28, 0) (29, 0) (30, 0) (31, 0) (32, 0) (33, 0) (34, 0) (35, 0) (36, 0) (37, 0) (38, 0) (39, 0) (40, 0) (41, 0) (42, 0) (43, 0) (44, 0) (45, 0) (46, 0) (47, 0) (48, 0) (49, 0) (50, 0) (51, 0) (52, 0) (53, 0) (54, 0) (55, 0) (56, 0) (57, 0) (58, 0) (59, 0) (60, 0) (61, 0) (62, 0) (63, 0)
row 2: (0, 0) (1, 0) (2, 0) (3, 1) (4, 0) (5, 0) (6, 0) (7, 0) (8, 0) (9, 0) (10, 0) (11, 0) (12, 0) (13, 0) (14, 0) (15, 0) (16, 0) (17, 0) (18, 0) (19, 0) (20, 0) (21, 0) (22, 0) (23, 0) (24, 0) (25, 0) (26, 0) (27, 0) (28, 0) (29, 0) (30, 0) (31, 0) (32, 0) (33, 0) (34, 0) (35, 0) (36, 0) (37, 0) (38, 0) (39, 0) (40, 0) (41, 0) (42, 0) (43, 0) (44, 0) (45, 0) (46, 0) (47, 0) (48, 0) (49, 0) (50, 0) (51, 0) (52, 0) (53, 0) (54, 0) (55, 0) (56, 0) (57, 0) (58, 0) (59, 0) (60, 0) (61, 0) (62, 0) (63, 0)
...分析矩阵内容
用LaTeX表示矩阵$M$如下：
$$
\vb*M = \mqty( 0 & 1 & 0 & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & 0 \\ -3072 & 0 & 1024 & 0 & \cdots & \cdots & \cdots & \cdots & \cdots & 0 \\ 0 & 0 & 0 & 1 & 0 & \cdots & \cdots & \cdots & \cdots & 0 \\ 1024 & 0 & -2048 & 0 & 1024 & 0 & \cdots & \cdots & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \ddots & \ddots & \ddots & \vdots & \vdots & \vdots \\ 0 & \cdots & 0 & 1024 & 0 & -2048 & 0 & 1024 & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \ddots & \ddots & \ddots & \vdots & \vdots \\ 0 & \cdots & \cdots & \cdots & \cdots & \cdots & 0 & 1 & 0 & 0 \\ 0 & \cdots & \cdots & 0 & 1024 & 0 & -2048 & 0 & 1024 & 0\\ 0 & \cdots & \cdots & \cdots & \cdots & \cdots & 0 & 0 & 0 & 1 \\ 0 & \cdots & \cdots & \cdots & \cdots & 0 & 1024 & 0 & -3072 & 0 )
$$
通过观察数据，矩阵元$A_{ij}$特点如下：
- 偶数行$\qty(0, 2, 4)$
- $j=i+i$处为1,其余为零
- 奇数行$\qty(0, 2, 4)$
- 第1行第1列是-3072, 第3列1024, 其余为0
- 第63行第60列为1024, 第62列为-3072, 其余为0
- 其余行满足如下关系
- $j=i-1$时, $A_{ij}=-2048$
- $j=i-3$时, $A_{ij}=1024$
- $j=i+1$时, $A_{ij}=1024$
与原方程对比
回顾原方程为如下形式：
$$\mqty(\pdv{}{t} \\ \pdv{}{t}) \mqty(f \\ g) = \mqty(0 & 1 \\ \pdv[2]{}{x} & 0) \mqty(f \\ g)$$
将其离散化有：
$$\mqty(\pdv{}{t} \\ \pdv{}{t} \\ \vdots \\ \pdv{}{t} \\ \pdv{}{t} \\ \vdots \\ \pdv{}{t} \\ \pdv{}{t})
\mqty(f_0 \\ g_0 \\ \vdots \\ f_i \\ g_i \\ \vdots \\ f_{31} \\ g_{31}) = \vb*M * \mqty(f_0 \\ g_0 \\ \vdots \\ f_i \\ g_i \\ \vdots \\ f_{31} \\ g_{31})$$
展开有：
$$
\mqty(\pdv{}{t} \\ \pdv{}{t} \\ \vdots \\ \pdv{}{t} \\ \pdv{}{t} \\ \vdots \\ \pdv{}{t} \\ \pdv{}{t})
\mqty(f_0 \\ g_0 \\ \vdots \\ f_i \\ g_i \\ \vdots \\ f_{31} \\ g_{31}) =
\mqty( 0 & 1 & 0 & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & 0 \\ -3072 & 0 & 1024 & 0 & \cdots & \cdots & \cdots & \cdots & \cdots & 0 \\ 0 & 0 & 0 & 1 & 0 & \cdots & \cdots & \cdots & \cdots & 0 \\ 1024 & 0 & -2048 & 0 & 1024 & 0 & \cdots & \cdots & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \ddots & \ddots & \ddots & \vdots & \vdots & \vdots \\ 0 & \cdots & 0 & 1024 & 0 & -2048 & 0 & 1024 & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \ddots & \ddots & \ddots & \vdots & \vdots \\ 0 & \cdots & \cdots & \cdots & \cdots & \cdots & 0 & 1 & 0 & 0 \\ 0 & \cdots & \cdots & 0 & 1024 & 0 & -2048 & 0 & 1024 & 0\\ 0 & \cdots & \cdots & \cdots & \cdots & \cdots & 0 & 0 & 0 & 1 \\ 0 & \cdots & \cdots & \cdots & \cdots & 0 & 1024 & 0 & -3072 & 0 )
* \mqty(f_0 \\ g_0 \\ \vdots \\ f_i \\ g_i \\ \vdots \\ f_{31} \\ g_{31})
$$
可以看到二阶微分使用了中心差分格式。
结果分析
原eigen-box的运行结果如下：
k ||Ax-kx||/||kx||
================= ==================
0 1.05144e-12+6.40000e+01i (64) -6.40000e+01+1.05144e-12i (64)
1 1.05144e-12-6.40000e+01i (64) 6.40000e+01+1.05144e-12i (64)
2 9.33011e-13+6.39229e+01i (63.9229) -6.39229e+01+9.33011e-13i (63.9229)
3 9.33011e-13-6.39229e+01i (63.9229) 6.39229e+01+9.33011e-13i (63.9229)
4 9.20997e-13+6.36918e+01i (63.6918) -6.36918e+01+9.20997e-13i (63.6918)
5 9.20997e-13-6.36918e+01i (63.6918) 6.36918e+01+9.20997e-13i (63.6918)
6 7.88002e-13+6.33073e+01i (63.3073) -6.33073e+01+7.88002e-13i (63.3073)
7 7.88002e-13-6.33073e+01i (63.3073) 6.33073e+01+7.88002e-13i (63.3073)
8 2.57072e-13+6.27703e+01i (62.7703) -6.27703e+01+2.57072e-13i (62.7703)
9 2.57072e-13-6.27703e+01i (62.7703) 6.27703e+01+2.57072e-13i (62.7703)
10 6.83981e-13+6.20820e+01i (62.082) -6.20820e+01+6.83981e-13i (62.082)
11 6.83981e-13-6.20820e+01i (62.082) 6.20820e+01+6.83981e-13i (62.082)
12 6.00998e-13+6.12442e+01i (61.2442) -6.12442e+01+6.00998e-13i (61.2442)
13 6.00998e-13-6.12442e+01i (61.2442) 6.12442e+01+6.00998e-13i (61.2442)
14 3.69158e-13+6.02588e+01i (60.2588) -6.02588e+01+3.69158e-13i (60.2588)
15 3.69158e-13-6.02588e+01i (60.2588) 6.02588e+01+3.69158e-13i (60.2588)
16 -4.19456e-14+5.91283e+01i (59.1283) -5.91283e+01-4.19456e-14i (59.1283)
17 -4.19456e-14-5.91283e+01i (59.1283) 5.91283e+01-4.19456e-14i (59.1283)
18 2.72976e-14+5.78553e+01i (57.8553) -5.78553e+01+2.72976e-14i (57.8553)
19 2.72976e-14-5.78553e+01i (57.8553) 5.78553e+01+2.72976e-14i (57.8553)
20 1.57958e-12+5.64430e+01i (56.443) -5.64430e+01+1.57958e-12i (56.443)
21 1.57958e-12-5.64430e+01i (56.443) 5.64430e+01+1.57958e-12i (56.443)
22 -2.93209e-13+5.48946e+01i (54.8946) -5.48946e+01-2.93209e-13i (54.8946)
23 -2.93209e-13-5.48946e+01i (54.8946) 5.48946e+01-2.93209e-13i (54.8946)
24 -1.35050e-12+5.32141e+01i (53.2141) -5.32141e+01-1.35050e-12i (53.2141)
25 -1.35050e-12-5.32141e+01i (53.2141) 5.32141e+01-1.35050e-12i (53.2141)
26 -2.71009e-13+5.14053e+01i (51.4053) -5.14053e+01-2.71009e-13i (51.4053)
27 -2.71009e-13-5.14053e+01i (51.4053) 5.14053e+01-2.71009e-13i (51.4053)
28 1.39111e-13+4.94727e+01i (49.4727) -4.94727e+01+1.39111e-13i (49.4727)
29 1.39111e-13-4.94727e+01i (49.4727) 4.94727e+01+1.39111e-13i (49.4727)
30 3.65052e-13+4.74209e+01i (47.4209) -4.74209e+01+3.65052e-13i (47.4209)
31 3.65052e-13-4.74209e+01i (47.4209) 4.74209e+01+3.65052e-13i (47.4209)
32 -7.31137e-13+4.52548e+01i (45.2548) -4.52548e+01-7.31137e-13i (45.2548)
33 -7.31137e-13-4.52548e+01i (45.2548) 4.52548e+01-7.31137e-13i (45.2548)
34 4.86552e-13+4.29798e+01i (42.9798) -4.29798e+01+4.86552e-13i (42.9798)
35 4.86552e-13-4.29798e+01i (42.9798) 4.29798e+01+4.86552e-13i (42.9798)
36 -4.60793e-13+4.06012e+01i (40.6012) -4.06012e+01-4.60793e-13i (40.6012)
37 -4.60793e-13-4.06012e+01i (40.6012) 4.06012e+01-4.60793e-13i (40.6012)
38 1.15241e-13+3.81248e+01i (38.1248) -3.81248e+01+1.15241e-13i (38.1248)
39 1.15241e-13-3.81248e+01i (38.1248) 3.81248e+01+1.15241e-13i (38.1248)
40 8.64686e-14+3.55565e+01i (35.5565) -3.55565e+01+8.64686e-14i (35.5565)
41 8.64686e-14-3.55565e+01i (35.5565) 3.55565e+01+8.64686e-14i (35.5565)
42 -2.48124e-13+3.29026e+01i (32.9026) -3.29026e+01-2.48124e-13i (32.9026)
43 -2.48124e-13-3.29026e+01i (32.9026) 3.29026e+01-2.48124e-13i (32.9026)
44 -4.28473e-13+3.01694e+01i (30.1694) -3.01694e+01-4.28473e-13i (30.1694)
45 -4.28473e-13-3.01694e+01i (30.1694) 3.01694e+01-4.28473e-13i (30.1694)
46 2.50383e-13+2.73635e+01i (27.3635) -2.73635e+01+2.50383e-13i (27.3635)
47 2.50383e-13-2.73635e+01i (27.3635) 2.73635e+01+2.50383e-13i (27.3635)
48 -6.17683e-14+2.44917e+01i (24.4917) -2.44917e+01-6.17683e-14i (24.4917)
49 -6.17683e-14-2.44917e+01i (24.4917) 2.44917e+01-6.17683e-14i (24.4917)
50 -5.11584e-13+2.15610e+01i (21.561) -2.15610e+01-5.11584e-13i (21.561)
51 -5.11584e-13-2.15610e+01i (21.561) 2.15610e+01-5.11584e-13i (21.561)
52 -8.63979e-14+1.85782e+01i (18.5782) -1.85782e+01-8.63979e-14i (18.5782)
53 -8.63979e-14-1.85782e+01i (18.5782) 1.85782e+01-8.63979e-14i (18.5782)
54 3.87100e-13+1.55507e+01i (15.5507) -1.55507e+01+3.87100e-13i (15.5507)
55 3.87100e-13-1.55507e+01i (15.5507) 1.55507e+01+3.87100e-13i (15.5507)
56 -9.68132e-14+1.24858e+01i (12.4858) -1.24858e+01-9.68132e-14i (12.4858)
57 -9.68132e-14-1.24858e+01i (12.4858) 1.24858e+01-9.68132e-14i (12.4858)
58 1.00364e-13+9.39075e+00i (9.39075) -9.39075e+00+1.00364e-13i (9.39075)
59 1.00364e-13-9.39075e+00i (9.39075) 9.39075e+00+1.00364e-13i (9.39075)
60 9.21981e-30+6.27310e+00i (6.2731) -6.27310e+00+9.21981e-30i (6.2731)
61 9.21981e-30-6.27310e+00i (6.2731) 6.27310e+00+9.21981e-30i (6.2731)
62 2.93182e-13+3.14033e+00i (3.14033) -3.14033e+00+2.93182e-13i (3.14033)
63 2.93182e-13-3.14033e+00i (3.14033) 3.14033e+00+2.93182e-13i (3.14033)编写matlab程序如下：
n = 64;
A = zeros(n, n);
A(1, 2) = 1;
A(2, 1) = -3072;
A(2, 3) = 1024;
for i = 3:2:n-1
A(i, i+1) = 1;
end
for i = 4:2:n-2
A(i, i-3) = 1024;
A(i, i+1) = 1024;
A(i, i-1) = -2048;
end
A(n, n-3) = 1024;
A(n, n-1) = -3072;
eigenvalues = eig(A);
disp('特征值:');
disp(eigenvalues);
figure;
plot(real(eigenvalues), imag(eigenvalues), 'bo');
grid on;
title('特征值在复平面上的分布');
xlabel('实部');
ylabel('虚部');输出如下：
-0.0000 +64.0000i
-0.0000 -64.0000i
0.0000 +63.9229i
0.0000 -63.9229i
-0.0000 +63.6918i
-0.0000 -63.6918i
0.0000 +63.3073i
0.0000 -63.3073i
0.0000 +62.7703i
0.0000 -62.7703i
0.0000 +62.0820i
0.0000 -62.0820i
-0.0000 +61.2442i
-0.0000 -61.2442i
-0.0000 +60.2588i
-0.0000 -60.2588i
-0.0000 +59.1283i
-0.0000 -59.1283i
-0.0000 +57.8553i
-0.0000 -57.8553i
-0.0000 +56.4430i
-0.0000 -56.4430i
-0.0000 +54.8946i
-0.0000 -54.8946i
-0.0000 +53.2141i
-0.0000 -53.2141i
-0.0000 +51.4053i
-0.0000 -51.4053i
0.0000 +49.4727i
0.0000 -49.4727i
-0.0000 +47.4209i
-0.0000 -47.4209i
-0.0000 +45.2548i
-0.0000 -45.2548i
0.0000 +42.9798i
0.0000 -42.9798i
-0.0000 +40.6012i
-0.0000 -40.6012i
0.0000 +38.1248i
0.0000 -38.1248i
0.0000 +35.5565i
0.0000 -35.5565i
0.0000 +32.9026i
0.0000 -32.9026i
-0.0000 +30.1694i
-0.0000 -30.1694i
-0.0000 +27.3635i
-0.0000 -27.3635i
0.0000 +24.4917i
0.0000 -24.4917i
0.0000 +21.5610i
0.0000 -21.5610i
-0.0000 +18.5782i
-0.0000 -18.5782i
0.0000 +15.5507i
0.0000 -15.5507i
-0.0000 +12.4858i
-0.0000 -12.4858i
-0.0000 + 9.3908i
-0.0000 - 9.3908i
-0.0000 + 6.2731i
-0.0000 - 6.2731i
-0.0000 + 3.1403i
-0.0000 - 3.1403i对应分布图如下
eigen-box和matlab的值完全对上。
第一个式子确定不是
$$\pdv[2]{f}{t}=\pdv[2]{f}{x}$$
吗？
已改