-
Notifications
You must be signed in to change notification settings - Fork 13
/
MattssonSvärdNordström2004_dev.jl
102 lines (84 loc) · 1.59 KB
/
MattssonSvärdNordström2004_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
using SymPy
N = 11
b = [symbols("b_$i", real=true) for i in 1:N]
B = diagm(b)
# 2nd order dissipation
D1 = zeros(Int, N, N)
for i in 1:N
D1[i,i] = 1
end
for i in 1:N-1
D1[i+1,i] = -1
end
D1[1,2] = 1
display(D1)
B1 = eye(Int, N); B1[1,1] = 0
display(B1)
display(D1' * B * D1)
display(D1' * B1 * D1)
# 4th order dissipation
D2 = zeros(Int, N, N)
for i in 1:N
D2[i,i] = -2
end
for i in 1:N-1
D2[i+1,i] = 1
D2[i,i+1] = 1
end
boundary_width = 1
for i in 1:boundary_width
D2[i,:] = D2[boundary_width+1,:]
D2[end-i+1,:] = D2[end-boundary_width,:]
end
display(D2)
B2 = eye(Int, N); B2[1,1] = B2[end,end] = 0
display(B2)
display(D2' * B * D2)
display(D2' * B2 * D2)
# 6th order dissipation
D3 = zeros(Int, N, N)
for i in 1:N
D3[i,i] = -3
end
for i in 1:N-1
D3[i+1,i] = 3
D3[i,i+1] = 1
end
for i in 1:N-2
D3[i+2,i] = -1
end
boundary_width = 2
for i in 1:boundary_width
D3[i,:] = D3[boundary_width+1,:]
end
for i in 1:boundary_width-1
D3[end-i+1,:] = D3[end-boundary_width+1,:]
end
display(D3)
B3 = eye(Int, N); B3[1,1] = B3[2,2]= B3[end,end] = 0
display(B3)
display(D3' * B * D3)
display(D3' * B3 * D3)
# 8th order dissipation
D4 = zeros(Int, N, N)
for i in 1:N
D4[i,i] = 6
end
for i in 1:N-1
D4[i+1,i] = -4
D4[i,i+1] = -4
end
for i in 1:N-2
D4[i+2,i] = 1
D4[i,i+2] = 1
end
boundary_width = 2
for i in 1:boundary_width
D4[i,:] = D4[boundary_width+1,:]
D4[end-i+1,:] = D4[end-boundary_width,:]
end
display(D4)
B4 = eye(Int, N); B4[1,1] = B4[2,2]= B4[end,end] = B4[end-1,end-1] = 0
display(B4)
display(D4' * B * D4)
display(D4' * B4 * D4)