Merge 7ae4991 into 44e810d

ranocha · Mar 25, 2020 · abe70bc · abe70bc
2 parents 44e810d + 7ae4991
commit abe70bc
Show file tree

Hide file tree

Showing 9 changed files with 1,887 additions and 501 deletions.
diff --git a/dev/Mattsson2017.jl b/dev/Mattsson2017.jl
@@ -0,0 +1,123 @@
+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)