-
Notifications
You must be signed in to change notification settings - Fork 12
/
MattssonAlmquistCarpenter2014Extended_dev.jl
79 lines (66 loc) · 2.7 KB
/
MattssonAlmquistCarpenter2014Extended_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
# order 2
h = [5//12, 7//6, 11//12, 1, 1]
e = [1//1, 0, 0, 0, 0]
H = Diagonal(h)
q1 = [0 7//12 -1//12 0 0]
q2 = [-q1[2] 0 7//12 0 0]
q3 = [-q1[3] -q2[3] 0 1//2 0]
q4 = [0 0 -1//2 0 1//2]
q5 = [0 0 0 -1//2 0]
Q = vcat(q1, q2, q3, q4, q5)
D1 = H \ (Q - 1//2*e*e'); display(D1)
x=1//1:length(h)
ks = 0:5
for k in ks
u=x.^k
println("k = ", k)
println(D1*u - k*x.^(k-1))
end
# order 4
h = [511//1600, 3971//2880, 929//1440, 587//480, 2659//2880, 14551//14400, 1, 1, 1]
e = [1//1, 0, 0, 0, 0, 0, 0, 0, 0]
H = Diagonal(h)
q1 = [0 436061//680400 -312391//5443200 -100109//907200 54107//5443200 263//15552 0 0 0]
q2 = [-q1[2] 0 238831//544320 290357//1088640 -125//5184 -220357//5443200 0 0 0]
q3 = [-q1[3] -q2[3] 0 3467//7776 -74827//1088640 3859//907200 0 0 0]
q4 = [-q1[4] -q2[4] -q3[4] 0 342799//544320 -149959//5443200 0 0 0]
q5 = [-q1[5] -q2[5] -q3[5] -q4[5] 0 428789//680400 -1//12 0 0]
q6 = [-q1[6] -q2[6] -q3[6] -q4[6] -q5[6] 0 2//3 -1//12 0]
q7 = [0 0 0 0 1//12 -2//3 0 2//3 -1//12]
q8 = [0 0 0 0 0 1//12 -2//3 0 2//3]
q9 = [0 0 0 0 0 0 1//12 -2//3 0]
Q = vcat(q1, q2, q3, q4, q5, q6, q7, q8, q9)
D1 = H \ (Q - 1//2*e*e'); display(D1)
x=1//1:length(h)
ks = 0:5
for k in ks
u=x.^k
println("k = ", k)
println(D1*u - k*x.^(k-1))
end
# order 6
h = [7497391//25401600, 1106227//725760, 105317//403200, 260179//145152, 303631//725760,
513973//403200, 671171//725760, 25631999//25401600, 1, 1, 1, 1]
e = [1//1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
H = Diagonal(h)
q1 = [0 14085349//21168000 -205189//8064000 -7583609//38102400 -89573//4233600 9//80 -39279943//1524096000 -13//2000 0 0 0 0]
q2 = [-q1[2] 0 56626691//304819200 97579603//152409600 1214239//16934400 -246914711//762048000 23380423//304819200 3//200 0 0 0 0]
q3 = [-q1[3] -q2[3] 0 22302157//101606400 -6048877//60963840 9737351//169344000 -508309//25401600 780937//304819200 0 0 0 0]
q4 = [-q1[4] -q2[4] -q3[4] 0 8120941//30481920 3467971//6350400 -42123859//304819200 -2075677//152409600 0 0 0 0]
q5 = [-q1[5] -q2[5] -q3[5] -q4[5] 0 29205107//152409600 15263287//304819200 -81127//3386880 0 0 0 0]
q6 = [-q1[6] -q2[6] -q3[6] -q4[6] -q5[6] 0 321139631//508032000 -49594423//762048000 1//60 0 0 0]
q7 = [-q1[7] -q2[7] -q3[7] -q4[7] -q5[7] -q6[7] 0 30841499//43545600 -3//20 1//60 0 0]
q8 = [-q1[8] -q2[8] -q3[8] -q4[8] -q5[8] -q6[8] -q7[8] 0 3//4 -3//20 1//60 0]
q9 = [0 0 0 0 0 -1//60 3//20 -3//4 0 3//4 -3//20 1//60]
q10 = [0 0 0 0 0 0 -1//60 3//20 -3//4 0 3//4 -3//20]
q11 = [0 0 0 0 0 0 0 -1//60 3//20 -3//4 0 3//4]
q12 = [0 0 0 0 0 0 0 0 -1//60 3//20 -3//4 0]
Q = vcat(q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12)
D1 = H \ (Q - 1//2*e*e'); display(D1)
x=1//1:length(h)
ks = 0:7
for k in ks
u=x.^k
println("k = ", k)
println(D1*u - k*x.^(k-1))
end