-
Notifications
You must be signed in to change notification settings - Fork 12
/
Mattsson2017.jl
241 lines (228 loc) · 13.6 KB
/
Mattsson2017.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
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
using LinearAlgebra
# order 2
Qp = [
-1//4 5//4 -1//2 0 0 0
-1//4 -5//4 2 -1//2 0 0
0 0 -3//2 2 -1//2 0
0 0 0 -3//2 2 -1//2
0 0 0 0 -5//4 5//4
0 0 0 0 -1//4 -1//4
]
Qm = Matrix(-Qp')
B = zero(Qp); B[1,1] = -1; B[end,end] = 1
M = Diagonal([1//4, 5//4, 1, 1, 5//4, 1//4])
Dp = M \ (Qp + B / 2)
Dm = M \ (Qm + B / 2)
# order 3
Qp = [
-1//12 3//4 -1//6 0 0 0
-5//12 -5//12 1 -1//6 0 0
0 -1//3 -1//2 1 -1//6 0
0 0 -1//3 -1//2 1 -1//6
0 0 0 -1//3 -5//12 3//4
0 0 0 0 -5//12 -1//12
]
Qm = Matrix(-Qp')
B = zero(Qp); B[1,1] = -1; B[end,end] = 1
M = Diagonal([5//12, 13//12, 1, 1, 13//12, 5//12])
Dp = M \ (Qp + B / 2)
Dm = M \ (Qm + B / 2)
# order 4
Qp = [
-1//48 205//288 -29//144 1//96 0 0 0 0 0 0 0
-169//288 -11//48 33//32 -43//144 1//12 0 0 0 0 0 0
11//144 -13//32 -29//48 389//288 -1//2 1//12 0 0 0 0 0
1//32 -11//144 -65//288 -13//16 3//2 -1//2 1//12 0 0 0 0
0 0 0 -1//4 -5//6 3//2 -1//2 1//12 0 0 0
0 0 0 0 -1//4 -5//6 3//2 -1//2 1//12 0 0
0 0 0 0 0 -1//4 -5//6 3//2 -1//2 1//12 0
0 0 0 0 0 0 -1//4 -13//16 389//288 -43//144 1//96
0 0 0 0 0 0 0 -65//288 -29//48 33//32 -29//144
0 0 0 0 0 0 0 -11//144 -13//32 -11//48 205//288
0 0 0 0 0 0 0 1//32 11//144 -169//288 -1//48
]
Qm = Matrix(-Qp')
B = zero(Qp); B[1,1] = -1; B[end,end] = 1
m = [49//144, 61//48, 41//48, 149//144]
M = Diagonal(vcat(m, ones(Int, size(Qp,1)-2*length(m)), m[end:-1:1]))
Dp = M \ (Qp + B / 2)
Dm = M \ (Qm + B / 2)
# order 5
Qp = [
-1//120 941//1440 -47//360 -7//480 0 0 0 0 0 0 0
-869//1440 -11//120 25//32 -43//360 1//30 0 0 0 0 0 0
29//360 -17//32 -29//120 1309//1440 -1//4 1//30 0 0 0 0 0
1//32 -11//360 -661//1440 -13//40 1 -1//4 1//30 0 0 0 0
0 0 1//20 -1//2 -1//3 1 -1//4 1//30 0 0 0
0 0 0 1//20 -1//2 -1//3 1 -1//4 1//30 0 0
0 0 0 0 1//20 -1//2 -1//3 1 -1//4 1//30 0
0 0 0 0 0 1//20 -1//2 -13//40 1309//1440 -43//360 -7//480
0 0 0 0 0 0 1//20 -661//1440 -29//120 25//32 -47//360
0 0 0 0 0 0 0 -11//360 -17//32 -11//120 941//1440
0 0 0 0 0 0 0 1//32 29//360 -869//1440 -1//120
]
Qm = Matrix(-Qp')
B = zero(Qp); B[1,1] = -1; B[end,end] = 1
# m = [251//720, 299//240, 41//48, 149//144]
m = [251//720, 299//240, 211//240, 739//720]
M = Diagonal(vcat(m, ones(Int, size(Qp,1)-2*length(m)), m[end:-1:1]))
Dp = M \ (Qp + B / 2)
Dm = M \ (Qm + B / 2)
# order 6
Qptmp = [
-265//128688 1146190567//1737288000 -1596619//18384000 -55265831//579096000 26269819//3474576000 2464501//144774000 0 0 0 0 0 0
-1116490567//1737288000 -8839//214480 190538869//347457600 102705469//694915200 413741//9651600 -191689861//3474576000 0 0 0 0 0 0
1096619//18384000 -135385429//347457600 -61067//321720 45137333//57909600 -253641811//694915200 70665929//579096000 -1//60 0 0 0 0 0
66965831//579096000 -208765789//694915200 -17623253//57909600 -18269//45960 410905829//347457600 -477953317//1158192000 2//15 -1//60 0 0 0 0
-49219819//3474576000 293299//9651600 26422771//694915200 -141938309//347457600 -346583//643440 2217185207//1737288000 -1//2 2//15 -1//60 0 0 0
-2374501//144774000 142906261//3474576000 -3137129//579096000 -29884283//1158192000 -630168407//1737288000 -3559//6128 4//3 -1//2 2//15 -1//60 0 0
0 0 0 0 1//30 -2//5 -7//12 4//3 -1//2 2//15 -1//60 0
0 0 0 0 0 1//30 -2//5 -7//12 4//3 -1//2 2//15 -1//60
]
bnd = size(Qptmp, 1)-2
Qp = zeros(eltype(Qptmp), 2*bnd+2, 2*bnd+2)
Qp[1:size(Qptmp,1), 1:size(Qptmp,2)] = Qptmp
for i in 1:bnd
Qp[end+1-i, end+1-size(Qptmp,1):end] = Qptmp[end:-1:1, i]
end
Qm = Matrix(-Qp')
B = zero(Qp); B[1,1] = -1; B[end,end] = 1
m = [13613//43200, 12049//8640, 535//864, 1079//864, 7841//8640, 43837//43200]
M = Diagonal(vcat(m, ones(Int, size(Qp,1)-2*length(m)), m[end:-1:1]))
Dp = M \ (Qp + B / 2)
Dm = M \ (Qm + B / 2)
# order 7
Qptmp = Rational{Int128}[
-265//300272 1587945773//2432203200 -1926361//25737600 -84398989//810734400 48781961//4864406400 3429119//202683600 0 0 0 0 0 0 0
-1570125773//2432203200 -26517//1501360 240029831//486440640 202934303//972881280 118207//13512240 -231357719//4864406400 0 0 0 0 0 0 0
1626361//25737600 -206937767//486440640 -61067//750680 49602727//81073440 -43783933//194576256 51815011//810734400 -1//140 0 0 0 0 0 0
91418989//810734400 -53314099//194576256 -33094279//81073440 -18269//107240 440626231//486440640 -365711063//1621468800 1//15 -1//140 0 0 0 0 0
-62551961//4864406400 799//35280 82588241//972881280 -279245719//486440640 -346583//1501360 2312302333//2432203200 -3//10 1//15 -1//140 0 0 0 0
-3375119//202683600 202087559//4864406400 -11297731//810734400 61008503//1621468800 -1360092253//2432203200 -10677//42896 1 -3//10 1//15 -1//140 0 0 0
0 0 0 -1//105 1//10 -3//5 -1//4 1 -3//10 1//15 -1//140 0 0
0 0 0 0 -1//105 1//10 -3//5 -1//4 1 -3//10 1//15 -1//140 0
0 0 0 0 0 -1//105 1//10 -3//5 -1//4 1 -3//10 1//15 -1//140
]
bnd = size(Qptmp, 1)-3
Qp = zeros(eltype(Qptmp), 2*bnd+3, 2*bnd+3)
Qp[1:size(Qptmp,1), 1:size(Qptmp,2)] = Qptmp
for i in 1:bnd
Qp[end+1-i, end+1-size(Qptmp,1):end] = Qptmp[end:-1:1, i]
end
Qm = Matrix(-Qp')
B = zero(Qp); B[1,1] = -1; B[end,end] = 1
m = [19087//60480, 84199//60480, 18869//30240, 37621//30240, 55031//60480, 61343//60480]
M = Diagonal(vcat(m, ones(Int, size(Qp,1)-2*length(m)), m[end:-1:1]))
Dp = M \ (Qp + B / 2)
Dm = M \ (Qm + B / 2)
function print_coefficients(D, m, bnd, suffix = "")
println(" left_boundary$(suffix) = (")
for i in 1:bnd
start = findfirst(!iszero, D[i, :])
len = findfirst(iszero, D[i, :]) - start
println(" DerivativeCoefficientRow{T,", start, ",", len, "}(SVector(T(", D[i, 1], "),")
for j in (start+1):(start+len-2)
println(" T(", D[i, j], "),")
end
println(" T(", D[i, start+len-1], "), )),")
end
println(" )")
println(" right_boundary$(suffix) = (")
for i in 1:bnd
start = findfirst(!iszero, D[end+1-i, end:-1:1])
len = findfirst(iszero, D[end+1-i, end:-1:1]) - start
println(" DerivativeCoefficientRow{T,", start, ",", len, "}(SVector(T(", D[end+1-i, end], "),")
for j in (start+1):(start+len-2)
println(" T(", D[end+1-i, end+1-j], "),")
end
println(" T(", D[end+1-i, end+1-(start+len-1)], "), )),")
end
println(" )")
start = findfirst(!iszero, D[bnd+1, :])
len = findfirst(iszero, D[bnd+1, start:end]) - 1
upper_coef = D[bnd+1, (bnd+2):(start+len-1)]
println(" upper_coef$(suffix) = SVector( T(", upper_coef[1], "),")
for i in 2:(length(upper_coef)-1)
println(" T(", upper_coef[i], "),")
end
println(" T(", upper_coef[end], "), )")
central_coef = D[bnd+1, bnd+1]
println(" central_coef_plus = T(", central_coef, ")")
lower_coef = D[bnd+1, bnd:-1:start]
println(" lower_coef$(suffix) = SVector( T(", lower_coef[1], "),")
for i in 2:(length(lower_coef)-1)
println(" T(", lower_coef[i], "),")
end
println(" T(", lower_coef[end], "), )")
println(" left_weights = SVector( T(", m[1], "),")
for i in 2:(length(m)-1)
println(" T(", m[i], "),")
end
println(" T(", m[end], "), )")
println(" right_weights = left_weights")
println(" left_boundary_derivatives = Tuple{}()")
println(" right_boundary_derivatives = left_boundary_derivatives")
end
# order 8
Qptmp = Rational{Int128}[
-16683//63018560 29345822969987//43354248537600 -2734625426591//40644608004000 -113480208109603//780376473676800 -830250230261//26012549122560 2500519492033//32515686403200 2235718279643//390188236838400 -388481888477//26543417472000 #==# 0 0 0 0 0 0 0 0
-29227665839987//43354248537600 -493793//63018560 8302717120817//26543417472000 3739408501537//9290196115200 2684481534461//13935294172800 -4450185662513//18580392230400 -1221838279381//37160784460800 90595000956023//1950941184192000 #==# 0 0 0 0 0 0 0 0
2505689537591//40644608004000 -7312922392817//26543417472000 -69332623//1323389760 10994933811709//18580392230400 -9270952411151//18580392230400 3191238635141//20644880256000 4442211176987//92901961152000 -940661365031//32515686403200 #==# 0 0 0 0 0 0 0 0
118016946570403//780376473676800 -4173878828737//9290196115200 -7990503962509//18580392230400 -207799621//1323389760 2044021560341//2477385630720 511197701761//18580392230400 1237681717213//13935294172800 -7784834666617//130062745612800 #==# 1//280 0 0 0 0 0 0 0
68609076271//2364777192960 -2235651762161//13935294172800 6527681584751//18580392230400 -1115980068821//2477385630720 -55386253//189055680 3208334350649//3716078446080 -407569013461//844563283200 136474842626653//780376473676800 #==# -1//28 1//280 0 0 0 0 0 0
-2487637785013//32515686403200 4244231077313//18580392230400 -1550378843141//20644880256000 -5726967564961//18580392230400 -1017898941929//3716078446080 -526653889//1323389760 45241297077547//37160784460800 -2279608411897//5080576000500 #==# 1//6 -1//28 1//280 0 0 0 0 0
-2164019088443//390188236838400 1263196075861//37160784460800 -6600697610987//92901961152000 556610591687//13935294172800 926842346471//9290196115200 -18757693936747//37160784460800 -584765899//1323389760 204646287449//168431424000 #==# -1//2 1//6 -1//28 1//280 0 0 0 0
387091928477//26543417472000 -90231551688023//1950941184192000 1032404418251//32515686403200 3502353445417//130062745612800 -15385068876253//780376473676800 262499068919//10161152001000 -867004691939//1852745664000 -85017967//189055680 #==# 5//4 -1//2 1//6 -1//28 1//280 0 0 0
# regular stencil starting in row 9
0 0 0 0 0 -1//168 1//14 -1//2 -9//20 5//4 -1//2 1//6 -1//28 1//280 0 0
0 0 0 0 0 0 -1//168 1//14 -1//2 -9//20 5//4 -1//2 1//6 -1//28 1//280 0
0 0 0 0 0 0 0 -1//168 1//14 -1//2 -9//20 5//4 -1//2 1//6 -1//28 1//280
]
bnd = size(Qptmp, 1)-3
Qp = zeros(eltype(Qptmp), 2*bnd+3, 2*bnd+3)
Qp[1:size(Qptmp,1), 1:size(Qptmp,2)] = Qptmp
for i in 1:bnd
Qp[end+1-i, end+1-size(Qptmp,1):end] = Qptmp[end:-1:1, i]
end
Qm = Matrix(-Qp')
B = zero(Qp); B[1,1] = -1; B[end,end] = 1
m = [7489399//25401600, 5537831//3628800, 103373//403200,
261259//145152, 298231//725760, 515917//403200,
3349159//3628800, 25639991//25401600]
M = Diagonal(vcat(m, ones(Int, size(Qp,1)-2*length(m)), m[end:-1:1]))
Dp = M \ (Qp + B / 2)
Dm = M \ (Qm + B / 2)
print_coefficients(Dp, m, bnd, "_plus")
print_coefficients(Dm, m, bnd, "_minus")
# order 9
Qptmp = Rational{Int128}[
-5561//47263920 4186300102421//6193464076800 -377895002003//5806372572000 -16485548951749//111482353382400 -113245973003//3716078446080 355360297339//4645098057600 321012170669//55741176691200 -388397049437//26543417472000 #==# 0 0 0 0 0 0 0 0 0
-4178798062421//6193464076800 -493793//141791760 725405227507//2413037952000 3904159533697//9290196115200 2483046570341//13935294172800 -4336328670953//18580392230400 -1258688487061//37160784460800 12931584852209//278705883456000 #==# 0 0 0 0 0 0 0 0 0
363359390003//5806372572000 -7539548734577//26543417472000 -69332623//2977626960 9994352248429//18580392230400 -8195655811631//18580392230400 7361486640463//61934640768000 5539855071347//92901961152000 -12898722943//422281641600 #==# 0 0 0 0 0 0 0 0 0
16773595838149//111482353382400 -372477950627//844563283200 -8659050093229//18580392230400 -207799621//2977626960 1734921317461//2477385630720 2530020015841//18580392230400 441856623253//13935294172800 -115132773073//2654341747200 #==# 1//630 0 0 0 0 0 0 0 0
108449122763//3716078446080 -2283566671541//13935294172800 6976424333231//18580392230400 -440819477447//825795210240 -55386253//425375280 2479572560009//3716078446080 -40258468963//120651897600 11808221047099//111482353382400 #==# -1//56 1//630 0 0 0 0 0 0 0
-32231128289//422281641600 4244793299753//18580392230400 -5173673584463//61934640768000 -4848139955041//18580392230400 -1506045711689//3716078446080 -526653889//2977626960 36411368691307//37160784460800 -825434105779//2903186286000 #==# 2//21 -1//56 1//630 0 0 0 0 0 0
-316459841069//55741176691200 1277069729941//37160784460800 -6499182375347//92901961152000 355606625147//13935294172800 1519272420551//9290196115200 -2240079855137//3378253132800 -584765899//2977626960 2301241355533//2382101568000 #==# -1//3 2//21 -1//56 1//630 0 0 0 0 0
387779289437//26543417472000 -12908508708209//278705883456000 147710908133//4645098057600 534025841911//18580392230400 -4119981443899//111482353382400 279819152779//2903186286000 -1510324515533//2382101568000 -85017967//425375280 #==# 1 -1//3 2//21 -1//56 1//630 0 0 0 0
# regular stencil starting in row 9
0 0 0 0 1//504 -1//42 1//7 -2//3 -1//5 1 -1//3 2//21 -1//56 1//630 0 0 0
0 0 0 0 0 1//504 -1//42 1//7 -2//3 -1//5 1 -1//3 2//21 -1//56 1//630 0 0
0 0 0 0 0 0 1//504 -1//42 1//7 -2//3 -1//5 1 -1//3 2//21 -1//56 1//630 0
0 0 0 0 0 0 0 1//504 -1//42 1//7 -2//3 -1//5 1 -1//3 2//21 -1//56 1//630
]
bnd = size(Qptmp, 1)-4
Qp = zeros(eltype(Qptmp), 2*bnd+4, 2*bnd+4)
Qp[1:size(Qptmp,1), 1:size(Qptmp,2)] = Qptmp
for i in 1:bnd
Qp[end+1-i, end+1-size(Qptmp,1):end] = Qptmp[end:-1:1, i]
end
Qm = Matrix(-Qp')
B = zero(Qp); B[1,1] = -1; B[end,end] = 1
m = [1070017//3628800, 5537111//3628800, 103613//403200, 261115//145152,
298951//725760, 515677//403200, 3349879//3628800, 3662753//3628800]
M = Diagonal(vcat(m, ones(Int, size(Qp,1)-2*length(m)), m[end:-1:1]))
Dp = M \ (Qp + B / 2)
Dm = M \ (Qm + B / 2)
e = ones(eltype(Dp), size(Dp, 2))
print_coefficients(Dp, m, bnd, "_plus")
print_coefficients(Dm, m, bnd, "_minus")