-
Notifications
You must be signed in to change notification settings - Fork 12
/
MattssonAlmquistVanDerWeide2018Minimal_dev.jl
127 lines (114 loc) · 2.79 KB
/
MattssonAlmquistVanDerWeide2018Minimal_dev.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
using LinearAlgebra
# order 4
h = [
2.6864248295847e-01,
1.0094667153500e+00,
9.9312068011715e-01,
1, 1, 1, 1, 1
]
e = [1//1, 0, 0, 0, 0, 0, 0, 0]
H = Diagonal(h)
q1 = [0,
6.1697245625434e-01,
-1.1697245625434e-01,
0.0000000000000e+00,
0, 0, 0, 0]'
q2 = [-q1[2], 0,
7.0030578958767e-01,
-8.3333333333333e-02,
0.0000000000000e+00,
0, 0, 0]'
q3 = [-q1[3], -q2[3], 0,
2//3, -1//12, 0, 0, 0]'
q4 = [-q1[4], -q2[4], -q3[4], 0,
2//3, -1//12, 0, 0]'
q5 = [-q1[5] -q2[5] -q3[5] -q4[5] 0 2//3 -1//12 0]
q6 = [0 0 0 1//12 -2//3 0 2//3 -1//12]
q7 = [0 0 0 0 1//12 -2//3 0 2//3]
q8 = [0 0 0 0 0 1//12 -2//3 0]
Q = vcat(q1, q2, q3, q4, q5, q6, q7, q8)
D1 = H \ (Q - 1//2*e*e'); display(D1)
# order 6
h = [
1.2740260779883e-01,
6.1820981002054e-01,
9.4308973897679e-01,
1.0093019060199e+00,
9.9884825610465e-01,
1, 1, 1, 1
]
e = [1//1, 0, 0, 0, 0, 0, 0, 0, 0]
H = Diagonal(h)
q1 = [0,
6.3217364546846e-01,
-1.6411963429825e-01,
3.6495407984639e-02,
-4.5494191548490e-03,
0, 0, 0, 0]'
q2 = [-q1[2], 0,
8.0515625504417e-01,
-2.0755653563249e-01,
3.4573926056780e-02,
0, 0, 0, 0]'
q3 = [-q1[3], -q2[3], 0,
7.9402676057785e-01,
-1.6965680649860e-01,
1.6666666666667e-02,
0, 0, 0]'
q4 = [-q1[4], -q2[4], -q3[4], 0,
7.5629896626333e-01,
-3//20, 1//60, 0, 0]'
q5 = [-q1[5] -q2[5] -q3[5] -q4[5] 0 3//4 -3//20 1//60 0]
q6 = [0 0 -1//60 3//20 -3//4 0 3//4 -3//20 1//60]
q7 = [0 0 0 -1//60 3//20 -3//4 0 3//4 -3//20]
q8 = [0 0 0 0 -1//60 3//20 -3//4 0 3//4]
q9 = [0 0 0 0 0 -1//60 3//20 -3//4 0]
Q = vcat(q1, q2, q3, q4, q5, q6, q7, q8, q9)
D1 = H \ (Q - 1//2*e*e'); display(D1)
# order 8
h = [
1.4523997892351e-01,
7.6864793350174e-01,
9.9116487068535e-01,
9.9992473335107e-01,
1.0002097054636e+00,
9.9996591555866e-01,
1, 1, 1, 1
]
e = [1//1, 0, 0, 0, 0, 0, 0, 0, 0, 0]
H = Diagonal(h)
q1 = [0,
6.6697342753834e-01,
-2.2919342278749e-01,
7.4283116457276e-02,
-1.2020661178873e-02,
-4.2460029252999e-05,
0, 0, 0, 0]'
q2 = [-q1[2], 0,
8.8241196934163e-01,
-2.6653314104602e-01,
5.5302527504316e-02,
-4.2079282615860e-03,
0, 0, 0, 0]'
q3 = [-q1[3], -q2[3], 0,
8.2904844081126e-01,
-2.1156614214635e-01,
3.9307676460659e-02,
-3.5714285714286e-03,
0, 0, 0]'
q4 = [-q1[4], -q2[4], -q3[4], 0,
8.0305501223679e-01,
-2.0078040553808e-01,
3.8095238095238e-02,
-3.5714285714286e-03,
0, 0]'
q5 = [-q1[5], -q2[5], -q3[5], -q4[5], 0,
8.0024692689207e-01,
-1//5, 4//105, -1//280, 0]'
q6 = [-q1[6] -q2[6] -q3[6] -q4[6] -q5[6] 0 4//5 -1//5 4//105 -1//280]
q7 = [0 0 1//280 -4//105 1//5 -4//5 0 4//5 -1//5 4//105]
q8 = [0 0 0 1//280 -4//105 1//5 -4//5 0 4//5 -1//5]
q9 = [0 0 0 0 1//280 -4//105 1//5 -4//5 0 4//5]
q10 = [0 0 0 0 0 1//280 -4//105 1//5 -4//5 0]
Q = vcat(q1, q2, q3, q4, q5, q6, q7, q8, q9, q10)
D1 = H \ (Q - 1//2*e*e'); display(D1)