Mercurial > repos > public > sbplib

comparison +sbp/+implementations/d2_noneq_variable_12.m @ 1325:1b0f2415237f feature/D2_boundary_opt

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add variable coefficient boundary-optimized second derivatives.
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Sun, 13 Feb 2022 19:32:34 +0100
parents
children c2d716c4f1ed
comparison
equal deleted inserted replaced
1301:8978521b0f06 1325:1b0f2415237f
1 function [H, HI, D1, D2, DI] = d2_noneq_variable_12(N, h, options)
2 % N: Number of grid points
3 % h: grid spacing
4 % options: struct containing options for constructing the operator
5 % current options are:
6 % options.stencil_type ('minimal','nonminimal','wide')
7 % options.AD ('upwind', 'op')
8
9 % BP: Number of boundary points
10 % order: Accuracy of interior stencil
11 BP = 12;
12 order = 12;
13
14 %%%% Norm matrix %%%%%%%%
15 P = zeros(BP, 1);
16 P0 = 1.0000000000011e-01;
17 P1 = 5.9616216757547e-01;
18 P2 = 9.9065699844442e-01;
19 P3 = 1.2512548713913e+00;
20 P4 = 1.3316678989403e+00;
21 P5 = 1.2093375037721e+00;
22 P6 = 1.0236491595704e+00;
23 P7 = 9.9685258909811e-01;
24 P8 = 1.0004766563923e+00;
25 P9 = 9.9993617879146e-01;
26 P10 = 1.0000063122914e+00;
27 P11 = 9.9999966373260e-01;
28
29 for i = 0:BP - 1
30 P(i + 1) = eval(['P' num2str(i)]);
31 end
32
33 Hv = ones(N, 1);
34 Hv(1:BP) = P;
35 Hv(end - BP + 1:end) = flip(P);
36 Hv = h * Hv;
37 H = spdiags(Hv, 0, N, N);
38 HI = spdiags(1 ./ Hv, 0, N, N);
39 %%%%%%%%%%%%%%%%%%%%%%%%%
40
41 %%%% Q matrix %%%%%%%%%%%
42
43 % interior stencil
44 d = [1/5544, -1/385, 1/56, -5/63, 15/56, -6/7, 0, 6/7, -15/56, 5/63, -1/56, 1/385, -1/5544];
45 d = repmat(d, N, 1);
46 Q = spdiags(d, -order / 2:order / 2, N, N);
47
48 % Boundaries
49 Q0_0 = -5.0000000000000e-01;
50 Q0_1 = 6.7597560728423e-01;
51 Q0_2 = -2.6859785384416e-01;
52 Q0_3 = 1.4850302678903e-01;
53 Q0_4 = -8.7976689586154e-02;
54 Q0_5 = 4.1833336322613e-02;
55 Q0_6 = -2.2216684976993e-03;
56 Q0_7 = -1.5910034062022e-02;
57 Q0_8 = 1.1296706376589e-02;
58 Q0_9 = -3.1823678285130e-03;
59 Q0_10 = 2.4843594063649e-04;
60 Q0_11 = 3.1501105449828e-05;
61 Q0_12 = 0.0000000000000e+00;
62 Q0_13 = 0.0000000000000e+00;
63 Q0_14 = 0.0000000000000e+00;
64 Q0_15 = 0.0000000000000e+00;
65 Q0_16 = 0.0000000000000e+00;
66 Q0_17 = 0.0000000000000e+00;
67 Q1_0 = -6.7597560728423e-01;
68 Q1_1 = 0.0000000000000e+00;
69 Q1_2 = 9.5424013647146e-01;
70 Q1_3 = -4.3389334603464e-01;
71 Q1_4 = 2.4285669347653e-01;
72 Q1_5 = -1.1443465137214e-01;
73 Q1_6 = 8.5942765682435e-03;
74 Q1_7 = 4.0290424215772e-02;
75 Q1_8 = -2.9396383714543e-02;
76 Q1_9 = 8.5601827834256e-03;
77 Q1_10 = -7.8128092862319e-04;
78 Q1_11 = -6.0444181254875e-05;
79 Q1_12 = 0.0000000000000e+00;
80 Q1_13 = 0.0000000000000e+00;
81 Q1_14 = 0.0000000000000e+00;
82 Q1_15 = 0.0000000000000e+00;
83 Q1_16 = 0.0000000000000e+00;
84 Q1_17 = 0.0000000000000e+00;
85 Q2_0 = 2.6859785384416e-01;
86 Q2_1 = -9.5424013647146e-01;
87 Q2_2 = 0.0000000000000e+00;
88 Q2_3 = 9.7065114311923e-01;
89 Q2_4 = -4.3205328628292e-01;
90 Q2_5 = 1.9549970932735e-01;
91 Q2_6 = -2.4406885385172e-02;
92 Q2_7 = -5.5737279079895e-02;
93 Q2_8 = 4.3772338637753e-02;
94 Q2_9 = -1.3727655130726e-02;
95 Q2_10 = 1.6271304373071e-03;
96 Q2_11 = 1.7066984372933e-05;
97 Q2_12 = 0.0000000000000e+00;
98 Q2_13 = 0.0000000000000e+00;
99 Q2_14 = 0.0000000000000e+00;
100 Q2_15 = 0.0000000000000e+00;
101 Q2_16 = 0.0000000000000e+00;
102 Q2_17 = 0.0000000000000e+00;
103 Q3_0 = -1.4850302678903e-01;
104 Q3_1 = 4.3389334603464e-01;
105 Q3_2 = -9.7065114311923e-01;
106 Q3_3 = 0.0000000000000e+00;
107 Q3_4 = 9.5375878629204e-01;
108 Q3_5 = -3.6113954384951e-01;
109 Q3_6 = 6.9749289223875e-02;
110 Q3_7 = 6.5161366516465e-02;
111 Q3_8 = -6.0325702283960e-02;
112 Q3_9 = 2.1188913621662e-02;
113 Q3_10 = -3.2632650250470e-03;
114 Q3_11 = 1.3097937809499e-04;
115 Q3_12 = 0.0000000000000e+00;
116 Q3_13 = 0.0000000000000e+00;
117 Q3_14 = 0.0000000000000e+00;
118 Q3_15 = 0.0000000000000e+00;
119 Q3_16 = 0.0000000000000e+00;
120 Q3_17 = 0.0000000000000e+00;
121 Q4_0 = 8.7976689586154e-02;
122 Q4_1 = -2.4285669347653e-01;
123 Q4_2 = 4.3205328628292e-01;
124 Q4_3 = -9.5375878629204e-01;
125 Q4_4 = 0.0000000000000e+00;
126 Q4_5 = 8.8676146394834e-01;
127 Q4_6 = -2.1292503103800e-01;
128 Q4_7 = -4.6037018833218e-02;
129 Q4_8 = 7.4338719466734e-02;
130 Q4_9 = -3.1217656663809e-02;
131 Q4_10 = 6.1239492854797e-03;
132 Q4_11 = -4.5892226603067e-04;
133 Q4_12 = 0.0000000000000e+00;
134 Q4_13 = 0.0000000000000e+00;
135 Q4_14 = 0.0000000000000e+00;
136 Q4_15 = 0.0000000000000e+00;
137 Q4_16 = 0.0000000000000e+00;
138 Q4_17 = 0.0000000000000e+00;
139 Q5_0 = -4.1833336322613e-02;
140 Q5_1 = 1.1443465137214e-01;
141 Q5_2 = -1.9549970932735e-01;
142 Q5_3 = 3.6113954384951e-01;
143 Q5_4 = -8.8676146394834e-01;
144 Q5_5 = 0.0000000000000e+00;
145 Q5_6 = 7.7461223007026e-01;
146 Q5_7 = -1.0609547334165e-01;
147 Q5_8 = -4.4853791547749e-02;
148 Q5_9 = 3.2436468405486e-02;
149 Q5_10 = -8.4387621360184e-03;
150 Q5_11 = 8.5964292632428e-04;
151 Q5_12 = 0.0000000000000e+00;
152 Q5_13 = 0.0000000000000e+00;
153 Q5_14 = 0.0000000000000e+00;
154 Q5_15 = 0.0000000000000e+00;
155 Q5_16 = 0.0000000000000e+00;
156 Q5_17 = 0.0000000000000e+00;
157 Q6_0 = 2.2216684976993e-03;
158 Q6_1 = -8.5942765682435e-03;
159 Q6_2 = 2.4406885385172e-02;
160 Q6_3 = -6.9749289223875e-02;
161 Q6_4 = 2.1292503103800e-01;
162 Q6_5 = -7.7461223007026e-01;
163 Q6_6 = 0.0000000000000e+00;
164 Q6_7 = 7.4758103262966e-01;
165 Q6_8 = -1.5730779067906e-01;
166 Q6_9 = 2.6517620342970e-02;
167 Q6_10 = -4.3175367549700e-03;
168 Q6_11 = 1.1092605832824e-03;
169 Q6_12 = -1.8037518037522e-04;
170 Q6_13 = 0.0000000000000e+00;
171 Q6_14 = 0.0000000000000e+00;
172 Q6_15 = 0.0000000000000e+00;
173 Q6_16 = 0.0000000000000e+00;
174 Q6_17 = 0.0000000000000e+00;
175 Q7_0 = 1.5910034062022e-02;
176 Q7_1 = -4.0290424215772e-02;
177 Q7_2 = 5.5737279079895e-02;
178 Q7_3 = -6.5161366516465e-02;
179 Q7_4 = 4.6037018833218e-02;
180 Q7_5 = 1.0609547334165e-01;
181 Q7_6 = -7.4758103262966e-01;
182 Q7_7 = 0.0000000000000e+00;
183 Q7_8 = 8.0975719267918e-01;
184 Q7_9 = -2.3568822398349e-01;
185 Q7_10 = 6.9373143801571e-02;
186 Q7_11 = -1.6606121869177e-02;
187 Q7_12 = 2.5974025974031e-03;
188 Q7_13 = -1.8037518037522e-04;
189 Q7_14 = 0.0000000000000e+00;
190 Q7_15 = 0.0000000000000e+00;
191 Q7_16 = 0.0000000000000e+00;
192 Q7_17 = 0.0000000000000e+00;
193 Q8_0 = -1.1296706376589e-02;
194 Q8_1 = 2.9396383714543e-02;
195 Q8_2 = -4.3772338637753e-02;
196 Q8_3 = 6.0325702283960e-02;
197 Q8_4 = -7.4338719466734e-02;
198 Q8_5 = 4.4853791547749e-02;
199 Q8_6 = 1.5730779067906e-01;
200 Q8_7 = -8.0975719267918e-01;
201 Q8_8 = 0.0000000000000e+00;
202 Q8_9 = 8.4765775072084e-01;
203 Q8_10 = -2.6369594097148e-01;
204 Q8_11 = 7.8759594625702e-02;
205 Q8_12 = -1.7857142857146e-02;
206 Q8_13 = 2.5974025974031e-03;
207 Q8_14 = -1.8037518037522e-04;
208 Q8_15 = 0.0000000000000e+00;
209 Q8_16 = 0.0000000000000e+00;
210 Q8_17 = 0.0000000000000e+00;
211 Q9_0 = 3.1823678285130e-03;
212 Q9_1 = -8.5601827834256e-03;
213 Q9_2 = 1.3727655130726e-02;
214 Q9_3 = -2.1188913621662e-02;
215 Q9_4 = 3.1217656663809e-02;
216 Q9_5 = -3.2436468405486e-02;
217 Q9_6 = -2.6517620342970e-02;
218 Q9_7 = 2.3568822398349e-01;
219 Q9_8 = -8.4765775072084e-01;
220 Q9_9 = 0.0000000000000e+00;
221 Q9_10 = 8.5631774953989e-01;
222 Q9_11 = -2.6769768119702e-01;
223 Q9_12 = 7.9365079365093e-02;
224 Q9_13 = -1.7857142857146e-02;
225 Q9_14 = 2.5974025974031e-03;
226 Q9_15 = -1.8037518037522e-04;
227 Q9_16 = 0.0000000000000e+00;
228 Q9_17 = 0.0000000000000e+00;
229 Q10_0 = -2.4843594063649e-04;
230 Q10_1 = 7.8128092862319e-04;
231 Q10_2 = -1.6271304373071e-03;
232 Q10_3 = 3.2632650250470e-03;
233 Q10_4 = -6.1239492854797e-03;
234 Q10_5 = 8.4387621360184e-03;
235 Q10_6 = 4.3175367549700e-03;
236 Q10_7 = -6.9373143801571e-02;
237 Q10_8 = 2.6369594097148e-01;
238 Q10_9 = -8.5631774953989e-01;
239 Q10_10 = 0.0000000000000e+00;
240 Q10_11 = 8.5712580212095e-01;
241 Q10_12 = -2.6785714285718e-01;
242 Q10_13 = 7.9365079365093e-02;
243 Q10_14 = -1.7857142857146e-02;
244 Q10_15 = 2.5974025974031e-03;
245 Q10_16 = -1.8037518037522e-04;
246 Q10_17 = 0.0000000000000e+00;
247 Q11_0 = -3.1501105449828e-05;
248 Q11_1 = 6.0444181254875e-05;
249 Q11_2 = -1.7066984372933e-05;
250 Q11_3 = -1.3097937809499e-04;
251 Q11_4 = 4.5892226603067e-04;
252 Q11_5 = -8.5964292632428e-04;
253 Q11_6 = -1.1092605832824e-03;
254 Q11_7 = 1.6606121869177e-02;
255 Q11_8 = -7.8759594625702e-02;
256 Q11_9 = 2.6769768119702e-01;
257 Q11_10 = -8.5712580212095e-01;
258 Q11_11 = 0.0000000000000e+00;
259 Q11_12 = 8.5714285714289e-01;
260 Q11_13 = -2.6785714285718e-01;
261 Q11_14 = 7.9365079365093e-02;
262 Q11_15 = -1.7857142857146e-02;
263 Q11_16 = 2.5974025974031e-03;
264 Q11_17 = -1.8037518037522e-04;
265
266 for i = 1:BP
267
268 for j = 1:BP
269 Q(i, j) = eval(['Q' num2str(i - 1) '_' num2str(j - 1)]);
270 Q(N + 1 - i, N + 1 - j) = -eval(['Q' num2str(i - 1) '_' num2str(j - 1)]);
271 end
272
273 end
274
275 %%%%%%%%%%%%%%%%%%%%%%%%%%%
276
277 %%% Undivided difference operators %%%%
278 % Closed with zeros at the first boundary nodes.
279 m = N;
280
281 DD_6 = (diag(ones(m - 3, 1), 3) - 6 * diag(ones(m - 2, 1), 2) + 15 * diag(ones(m - 1, 1), 1) - 20 * diag(ones(m, 1), 0) + 15 * diag(ones(m - 1, 1), -1) - 6 * diag(ones(m - 2, 1), -2) + diag(ones(m - 3, 1), -3));
282 DD_6(1:9, 1:12) = [0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0; 0.70504538217624485042e1 -0.15481773512478815978e2 0.15439628908176388584e2 -0.11355022328644135767e2 0.62553652361133177256e1 -0.23565505827757565576e1 0.44789845784655348861e0 0 0 0 0 0; 0 0.84178325750049836493e0 -0.30776567832262484696e1 0.55480476621608890820e1 -0.66451562589578036671e1 0.54681670852508283306e1 -0.26873907470793164527e1 0.55220578435115281189e0 0 0 0 0; 0 0 0.36124410255434790612e0 -0.18841869963752696615e1 0.49068751071627186058e1 -0.79710602635488243814e1 0.75771246506939812399e1 -0.36661050678180486359e1 0.67610846733109492706e0 0 0 0; 0 0 0 0.31883568286286899671e0 -0.22216567686182446798e1 0.72341822309820977868e1 -0.12215760254444741754e2 0.10698736280837039597e2 -0.46222617037529123987e1 0.80792453213389245238e0 0 0; 0 0 0 0 0.45451605753051033631e0 -0.36739526096585069959e1 0.11306936594922967126e2 -0.16769946247690595362e2 0.13179014442979029716e2 -0.54150410306154005497e1 0.91847279253199572926e0 0; 0 0 0 0 0 0.77355004999207313207e0 -0.54143198515959566716e1 0.14230335022999000858e2 -0.19303948318184081752e2 0.14605034473743418451e2 -0.58729419174319386168e1 0.98229054047748459948e0; ];
283 DD_6(m - 8:m, m - 11:m) = [0.98229054047748459948e0 -0.58729419174319386168e1 0.14605034473743418451e2 -0.19303948318184081752e2 0.14230335022999000858e2 -0.54143198515959566716e1 0.77355004999207313207e0 0 0 0 0 0; 0 0.91847279253199572926e0 -0.54150410306154005497e1 0.13179014442979029716e2 -0.16769946247690595362e2 0.11306936594922967126e2 -0.36739526096585069959e1 0.45451605753051033631e0 0 0 0 0; 0 0 0.80792453213389245238e0 -0.46222617037529123987e1 0.10698736280837039597e2 -0.12215760254444741754e2 0.72341822309820977868e1 -0.22216567686182446798e1 0.31883568286286899671e0 0 0 0; 0 0 0 0.67610846733109492706e0 -0.36661050678180486359e1 0.75771246506939812399e1 -0.79710602635488243814e1 0.49068751071627186058e1 -0.18841869963752696615e1 0.36124410255434790612e0 0 0; 0 0 0 0 0.55220578435115281189e0 -0.26873907470793164527e1 0.54681670852508283306e1 -0.66451562589578036671e1 0.55480476621608890820e1 -0.30776567832262484696e1 0.84178325750049836493e0 0; 0 0 0 0 0 0.44789845784655348861e0 -0.23565505827757565576e1 0.62553652361133177256e1 -0.11355022328644135767e2 0.15439628908176388584e2 -0.15481773512478815978e2 0.70504538217624485042e1; 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0; ];
284 DD_6 = sparse(DD_6);
285
286 DD_7 = (-diag(ones(m - 4, 1), -4) + 7 * diag(ones(m - 3, 1), -3) - 21 * diag(ones(m - 2, 1), -2) + 35 * diag(ones(m - 1, 1), -1) - 35 * diag(ones(m, 1), 0) + 21 * diag(ones(m - 1, 1), 1) - 7 * diag(ones(m - 2, 1), 2) + diag(ones(m - 3, 1), 3));
287 DD_7(1:10, 1:13) = [0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0; -0.70504538217624585762e1 0.16323556769979337662e2 -0.18517285691402663507e2 0.16903069990805048996e2 -0.12900521495071139822e2 0.78247176680265960663e1 -0.31352892049258744203e1 0.55220578435115360075e0 0 0 0 0 0; 0 -0.77136636008199086135e0 0.31512297676528003033e1 -0.68105129732006589582e1 0.10585680229488408490e2 -0.12315008394369015993e2 0.94058676147776075845e1 -0.38654404904580696832e1 0.61955060619091911792e0 0 0 0 0; 0 0 -0.32268062071305431283e0 0.19678458985349116625e1 -0.63675477999083230026e1 0.13582054493165302513e2 -0.17679957518285956226e2 0.12831367737363170226e2 -0.47327592713176644894e1 0.72167708116161363042e0 0 0 0; 0 0 0 -0.28975994984220671445e0 0.24321233335964448997e1 -0.99133841479577475282e1 0.21377580445278298070e2 -0.24963717988619759059e2 0.16177915963135193396e2 -0.56554717249372471666e1 0.83471406934702410357e0 0 0; 0 0 0 0 -0.43028213606032104517e0 0.42103703770684731501e1 -0.15829711232892153976e2 0.29347405933458541883e2 -0.30751033700284402671e2 0.18952643607153901924e2 -0.64293095477239701048e1 0.92991669927993083983e0 0; 0 0 0 0 0 -0.76177813656089336989e0 0.63167064935286161169e1 -0.19922469032198601202e2 0.33781909556822143066e2 -0.34078413772067976385e2 0.20555296711011785159e2 -0.68760337833423921963e1 0.98478196280731881078e0; ];
288 DD_7(m - 8:m, m - 12:m) = [-0.98478196280731881078e0 0.68760337833423921963e1 -0.20555296711011785159e2 0.34078413772067976385e2 -0.33781909556822143066e2 0.19922469032198601202e2 -0.63167064935286161169e1 0.76177813656089336989e0 0 0 0 0 0; 0 -0.92991669927993083983e0 0.64293095477239701048e1 -0.18952643607153901924e2 0.30751033700284402671e2 -0.29347405933458541883e2 0.15829711232892153976e2 -0.42103703770684731501e1 0.43028213606032104517e0 0 0 0 0; 0 0 -0.83471406934702410357e0 0.56554717249372471666e1 -0.16177915963135193396e2 0.24963717988619759059e2 -0.21377580445278298070e2 0.99133841479577475282e1 -0.24321233335964448997e1 0.28975994984220671445e0 0 0 0; 0 0 0 -0.72167708116161363042e0 0.47327592713176644894e1 -0.12831367737363170226e2 0.17679957518285956226e2 -0.13582054493165302513e2 0.63675477999083230026e1 -0.19678458985349116625e1 0.32268062071305431283e0 0 0; 0 0 0 0 -0.61955060619091911792e0 0.38654404904580696832e1 -0.94058676147776075845e1 0.12315008394369015993e2 -0.10585680229488408490e2 0.68105129732006589582e1 -0.31512297676528003033e1 0.77136636008199086135e0 0; 0 0 0 0 0 -0.55220578435115360075e0 0.31352892049258744203e1 -0.78247176680265960663e1 0.12900521495071139822e2 -0.16903069990805048996e2 0.18517285691402663507e2 -0.16323556769979337662e2 0.70504538217624585762e1; 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0; ];
289 DD_7 = sparse(DD_7);
290
291 DD_8 = (diag(ones(m - 4, 1), 4) - 8 * diag(ones(m - 3, 1), 3) + 28 * diag(ones(m - 2, 1), 2) - 56 * diag(ones(m - 1, 1), 1) + 70 * diag(ones(m, 1), 0) - 56 * diag(ones(m - 1, 1), -1) + 28 * diag(ones(m - 2, 1), -2) - 8 * diag(ones(m - 3, 1), -3) + diag(ones(m - 4, 1), -4));
292 DD_8(1:10, 1:14) = [0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0.70504538217624673892e1 -0.17094923130061349891e2 0.21668515459055490895e2 -0.23713582964005737596e2 0.23486201724559577669e2 -0.20139726062395637234e2 0.12541156819703497681e2 -0.44176462748092288060e1 0.61955060619091989235e0 0 0 0 0 0; 0 0.71430915910501867497e0 -0.32169486987428931518e1 0.81290324137612243914e1 -0.15699214646211457760e2 0.23981482952559874949e2 -0.25082313639406953559e2 0.15461761961832278733e2 -0.49564048495273529434e1 0.66829534663026066516e0 0 0 0 0; 0 0 0.29213206789956748018e0 -0.20438756549159061384e1 0.79665959467484077329e1 -0.21271097908734749058e2 0.35359915036571912453e2 -0.34216980632968453935e2 0.18931037085270657958e2 -0.57734166492929090433e1 0.75569070942147255138e0 0 0 0; 0 0 0 0.26637215914130374043e0 -0.26313682263734604148e1 0.12983764630211125301e2 -0.34204128712445276912e2 0.49927435977239518119e2 -0.43141109235027182388e2 0.22621886899748988667e2 -0.66777125547761928286e1 0.85485906228117671607e0 0 0; 0 0 0 0 0.41007336094698233166e0 -0.47386249189431794318e1 0.21106281643856205301e2 -0.46955849493533667013e2 0.61502067400568805342e2 -0.50540382952410405130e2 0.25717238190895880419e2 -0.74393335942394467187e1 0.93853036285882489957e0 0; 0 0 0 0 0 0.75161513192516571838e0 -0.72190931354612755621e1 0.26563292042931468269e2 -0.54051055290915428906e2 0.68156827544135952769e2 -0.54814124562698093757e2 0.27504135133369568785e2 -0.78782557024585504862e1 0.98665883917119316896e0; ];
293 DD_8(m - 9:m, m - 13:m) = [0.98665883917119316896e0 -0.78782557024585504862e1 0.27504135133369568785e2 -0.54814124562698093757e2 0.68156827544135952769e2 -0.54051055290915428906e2 0.26563292042931468269e2 -0.72190931354612755621e1 0.75161513192516571838e0 0 0 0 0 0; 0 0.93853036285882489957e0 -0.74393335942394467187e1 0.25717238190895880419e2 -0.50540382952410405130e2 0.61502067400568805342e2 -0.46955849493533667013e2 0.21106281643856205301e2 -0.47386249189431794318e1 0.41007336094698233166e0 0 0 0 0; 0 0 0.85485906228117671607e0 -0.66777125547761928286e1 0.22621886899748988667e2 -0.43141109235027182388e2 0.49927435977239518119e2 -0.34204128712445276912e2 0.12983764630211125301e2 -0.26313682263734604148e1 0.26637215914130374043e0 0 0 0; 0 0 0 0.75569070942147255138e0 -0.57734166492929090433e1 0.18931037085270657958e2 -0.34216980632968453935e2 0.35359915036571912453e2 -0.21271097908734749058e2 0.79665959467484077329e1 -0.20438756549159061384e1 0.29213206789956748018e0 0 0; 0 0 0 0 0.66829534663026066516e0 -0.49564048495273529434e1 0.15461761961832278733e2 -0.25082313639406953559e2 0.23981482952559874949e2 -0.15699214646211457760e2 0.81290324137612243914e1 -0.32169486987428931518e1 0.71430915910501867497e0 0; 0 0 0 0 0 0.61955060619091989235e0 -0.44176462748092288060e1 0.12541156819703497681e2 -0.20139726062395637234e2 0.23486201724559577669e2 -0.23713582964005737596e2 0.21668515459055490895e2 -0.17094923130061349891e2 0.70504538217624673892e1; 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0; ];
294 DD_8 = sparse(DD_8);
295
296 DD_9 = (-diag(ones(m - 5, 1), -5) + 9 * diag(ones(m - 4, 1), -4) - 36 * diag(ones(m - 3, 1), -3) + 84 * diag(ones(m - 2, 1), -2) - 126 * diag(ones(m - 1, 1), -1) + 126 * diag(ones(m, 1), 0) - 84 * diag(ones(m - 1, 1), 1) + 36 * diag(ones(m - 2, 1), 2) - 9 * diag(ones(m - 3, 1), 3) + diag(ones(m - 4, 1), 4));
297 DD_9(1:11, 1:15) = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; -0.70504538217624752230e1 0.17809232289166388354e2 -0.24885464157798411697e2 0.31842615377766997368e2 -0.39185416370771078968e2 0.44121209014955561206e2 -0.37623470459110493043e2 0.19879408236641529627e2 -0.55759554557182790312e1 0.66829534663026140771e0 0 0 0 0 0; 0 -0.66695396914459764542e0 0.32764459415493351213e1 -0.94984942131335544919e1 0.22096883550779531311e2 -0.42252556942280502185e2 0.56435205688665645507e2 -0.46385285885496836199e2 0.22303821822873088245e2 -0.60146581196723459865e1 0.70559212586023632254e0 0 0 0 0; 0 0 -0.26728718140337933088e0 0.21137687176984662199e1 -0.96966416347745772182e1 0.31341597391980823551e2 -0.63647847065829442415e2 0.76988206424179021354e2 -0.56793111255811973873e2 0.25980374921818090695e2 -0.68012163847932529624e1 0.78215606693622397877e0 0 0 0; 0 0 0 -0.24708804973573377055e0 0.28212553166921198174e1 -0.16439370707786854975e2 0.51306193068667915368e2 -0.89869384759031132614e2 0.97067495778811160373e2 -0.67865660699246966000e2 0.30049706496492867728e2 -0.76937315605305904447e1 0.87058511566721451662e0 0 0; 0 0 0 0 -0.39286387086790423439e0 0.52598319320161875843e1 -0.27136647827815121102e2 0.70433774240300500520e2 -0.11070372132102384962e3 0.11371586164292341154e3 -0.77151714572687641258e2 0.33477001174077510234e2 -0.84467732657294240961e1 0.94525186880633042479e0 0; 0 0 0 0 0 -0.74268863896636254047e0 0.81214797773939350074e1 -0.34152804055197602060e2 0.81076582936373143359e2 -0.12268228957944471498e3 0.12333178026607071095e3 -0.82512405400108706356e2 0.35452150661063477188e2 -0.88799295525407385207e1 0.98812358535685795524e0; ];
298 DD_9(m - 9:m, m - 14:m) = [-0.98812358535685795524e0 0.88799295525407385207e1 -0.35452150661063477188e2 0.82512405400108706356e2 -0.12333178026607071095e3 0.12268228957944471498e3 -0.81076582936373143359e2 0.34152804055197602060e2 -0.81214797773939350074e1 0.74268863896636254047e0 0 0 0 0 0; 0 -0.94525186880633042479e0 0.84467732657294240961e1 -0.33477001174077510234e2 0.77151714572687641258e2 -0.11371586164292341154e3 0.11070372132102384962e3 -0.70433774240300500520e2 0.27136647827815121102e2 -0.52598319320161875843e1 0.39286387086790423439e0 0 0 0 0; 0 0 -0.87058511566721451662e0 0.76937315605305904447e1 -0.30049706496492867728e2 0.67865660699246966000e2 -0.97067495778811160373e2 0.89869384759031132614e2 -0.51306193068667915368e2 0.16439370707786854975e2 -0.28212553166921198174e1 0.24708804973573377055e0 0 0 0; 0 0 0 -0.78215606693622397877e0 0.68012163847932529624e1 -0.25980374921818090695e2 0.56793111255811973873e2 -0.76988206424179021354e2 0.63647847065829442415e2 -0.31341597391980823551e2 0.96966416347745772182e1 -0.21137687176984662199e1 0.26728718140337933088e0 0 0; 0 0 0 0 -0.70559212586023632254e0 0.60146581196723459865e1 -0.22303821822873088245e2 0.46385285885496836199e2 -0.56435205688665645507e2 0.42252556942280502185e2 -0.22096883550779531311e2 0.94984942131335544919e1 -0.32764459415493351213e1 0.66695396914459764542e0 0; 0 0 0 0 0 -0.66829534663026140771e0 0.55759554557182790312e1 -0.19879408236641529627e2 0.37623470459110493043e2 -0.44121209014955561206e2 0.39185416370771078968e2 -0.31842615377766997368e2 0.24885464157798411697e2 -0.17809232289166388354e2 0.70504538217624752230e1; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; ];
299 DD_9 = sparse(DD_9);
300
301 DD_10 = (diag(ones(m - 5, 1), 5) - 10 * diag(ones(m - 4, 1), 4) + 45 * diag(ones(m - 3, 1), 3) - 120 * diag(ones(m - 2, 1), 2) + 210 * diag(ones(m - 1, 1), 1) - 252 * diag(ones(m, 1), 0) + 210 * diag(ones(m - 1, 1), -1) - 120 * diag(ones(m - 2, 1), -2) + 45 * diag(ones(m - 3, 1), -3) - 10 * diag(ones(m - 4, 1), -4) + diag(ones(m - 5, 1), -5));
302 DD_10(1:11, 1:16) = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0.70504538217624822734e1 -0.18476186258311004476e2 0.28161910099347774980e2 -0.41341109590900593201e2 0.61282299921550671561e2 -0.86373765957236149764e2 0.94058676147776232608e2 -0.66264694122138432090e2 0.27879777278591395156e2 -0.66829534663026140771e1 0.70559212586023702812e0 0 0 0 0 0; 0 0.62689419647750179604e0 -0.33308831364980034247e1 0.10914786591113557343e2 -0.29883886064802801254e2 0.69173811658980680918e2 -0.11287041137733129101e3 0.11596321471374209050e3 -0.74346072742910294151e2 0.30073290598361729932e2 -0.70559212586023632254e1 0.73517682146919258207e0 0 0 0 0; 0 0 0.24665297575614270036e0 -0.21786018467637256442e1 0.11551532389532584562e2 -0.44092342567416599454e2 0.10607974510971573736e3 -0.15397641284835804271e3 0.14198277813952993468e3 -0.86601249739393635650e2 0.34006081923966264812e2 -0.78215606693622397877e1 0.80337713279357912969e0 0 0 0; 0 0 0 0.23087142251463535911e0 -0.30031734426541638399e1 0.20275063024062368195e2 -0.73294561526668450525e2 0.14978230793171855436e3 -0.19413499155762232075e3 0.16966415174811741500e3 -0.10016568832164289243e3 0.38468657802652952223e2 -0.87058511566721451662e1 0.88321407619404757308e0 0 0; 0 0 0 0 0.37796280007363075810e0 -0.57748488744870271355e1 0.33920809784768901377e2 -0.10061967748614357217e3 0.18450620220170641603e3 -0.22743172328584682309e3 0.19287928643171910314e3 -0.11159000391359170078e3 0.42233866328647120480e2 -0.94525186880633042479e1 0.95064470121725563376e0 0; 0 0 0 0 0 0.73474076934244039806e0 -0.90238664193265944527e1 0.42691005068997002575e2 -0.11582368990910449051e3 0.20447048263240785831e3 -0.24666356053214142190e3 0.20628101350027176589e3 -0.11817383553687825729e3 0.44399647762703692603e2 -0.98812358535685795524e1 0.98929851729658394154e0; ];
303 DD_10(m - 10:m, m - 15:m) = [0.98929851729658394154e0 -0.98812358535685795524e1 0.44399647762703692603e2 -0.11817383553687825729e3 0.20628101350027176589e3 -0.24666356053214142190e3 0.20447048263240785831e3 -0.11582368990910449051e3 0.42691005068997002575e2 -0.90238664193265944527e1 0.73474076934244039806e0 0 0 0 0 0; 0 0.95064470121725563376e0 -0.94525186880633042479e1 0.42233866328647120480e2 -0.11159000391359170078e3 0.19287928643171910314e3 -0.22743172328584682309e3 0.18450620220170641603e3 -0.10061967748614357217e3 0.33920809784768901377e2 -0.57748488744870271355e1 0.37796280007363075810e0 0 0 0 0; 0 0 0.88321407619404757308e0 -0.87058511566721451662e1 0.38468657802652952223e2 -0.10016568832164289243e3 0.16966415174811741500e3 -0.19413499155762232075e3 0.14978230793171855436e3 -0.73294561526668450525e2 0.20275063024062368195e2 -0.30031734426541638399e1 0.23087142251463535911e0 0 0 0; 0 0 0 0.80337713279357912969e0 -0.78215606693622397877e1 0.34006081923966264812e2 -0.86601249739393635650e2 0.14198277813952993468e3 -0.15397641284835804271e3 0.10607974510971573736e3 -0.44092342567416599454e2 0.11551532389532584562e2 -0.21786018467637256442e1 0.24665297575614270036e0 0 0; 0 0 0 0 0.73517682146919258207e0 -0.70559212586023632254e1 0.30073290598361729932e2 -0.74346072742910294151e2 0.11596321471374209050e3 -0.11287041137733129101e3 0.69173811658980680918e2 -0.29883886064802801254e2 0.10914786591113557343e2 -0.33308831364980034247e1 0.62689419647750179604e0 0; 0 0 0 0 0 0.70559212586023702812e0 -0.66829534663026140771e1 0.27879777278591395156e2 -0.66264694122138432090e2 0.94058676147776232608e2 -0.86373765957236149764e2 0.61282299921550671561e2 -0.41341109590900593201e2 0.28161910099347774980e2 -0.18476186258311004476e2 0.70504538217624822734e1; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; ];
304 DD_10 = sparse(DD_10);
305
306 %[-1, 11, -55, 165, -330, 462, -462, 330, -165, 55, -11, 1, 0]
307 DD_11 = (-diag(ones(m - 6, 1), -6) + 11 * diag(ones(m - 5, 1), -5) - 55 * diag(ones(m - 4, 1), -4) + 165 * diag(ones(m - 3, 1), -3) - 330 * diag(ones(m - 2, 1), -2) + 462 * diag(ones(m - 1, 1), -1) - 462 * diag(ones(m, 1), 0) + 330 * diag(ones(m - 1, 1), 1) - 165 * diag(ones(m - 2, 1), 2) + 55 * diag(ones(m - 3, 1), 3) - 11 * diag(ones(m - 4, 1), 4) + diag(ones(m - 5, 1), 5));
308 DD_11(1:12, 1:17) = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; -0.70504538217624886828e1 0.19103080454788523638e2 -0.31492793235845807035e2 0.52255896182014198048e2 -0.91166185986353555693e2 0.15554757761621697209e3 -0.20692908752510771174e3 0.18222790883588068825e3 -0.10222585002150178224e3 0.36756244064664377424e2 -0.77615133844626073094e1 0.73517682146919325041e0 0 0 0 0 0; 0 -0.59247568548574727441e0 0.33811178542850964393e1 -0.12374519230919224800e2 0.39160480492663455349e2 -0.10704747746062038309e3 0.20692908752510736686e3 -0.25511907237023259909e3 0.20445170004300330891e3 -0.11026873219399300975e3 0.38807566922312997740e2 -0.80869450361611184028e1 0.75926914003985711330e0 0 0 0 0; 0 0 -0.22922038452398770506e0 0.22391799149842653766e1 -0.13526028855965613027e2 0.59818136859120542509e2 -0.16669674231526758728e3 0.28229009022198974497e3 -0.31236211190696585630e3 0.23815343678333249804e3 -0.12468896705454297098e3 0.43018583681492318832e2 -0.88371484607293704266e1 0.82079151707601599340e0 0 0 0; 0 0 0 -0.21701391114437616116e0 0.31781915325611402199e1 -0.24486327517266714070e2 0.10078002209916911947e3 -0.23537219817841487113e3 0.35591415118897425470e3 -0.37326113384585831300e3 0.27545564288451795418e3 -0.14105174527639415815e3 0.47882181361696798414e2 -0.97153548381345233039e1 0.89358450029368883150e0 0 0; 0 0 0 0 -0.36488508570871329747e0 0.62843543720560487741e1 -0.41458767514717546128e2 0.13835205654344741174e3 -0.28993831774553865375e3 0.41695815935738584232e3 -0.42433443014978202692e3 0.30687251076237717714e3 -0.15485750987170610843e3 0.51988852784348173364e2 -0.10457091713389811971e2 0.95506826122820715686e0 0; 0 0 0 0 0 -0.72758579432539352184e0 0.99262530612592538979e1 -0.52177895084329669814e2 0.15925757362501867446e3 -0.32131075842235520591e3 0.45221652764225927349e3 -0.45381822970059788496e3 0.32497804772641520756e3 -0.16279870846324687288e3 0.54346797194627187538e2 -0.10882283690262423357e2 0.99026190553785350209e0; ];
309 DD_11(m - 10:m, m - 16:m) = [-0.99026190553785350209e0 0.10882283690262423357e2 -0.54346797194627187538e2 0.16279870846324687288e3 -0.32497804772641520756e3 0.45381822970059788496e3 -0.45221652764225927349e3 0.32131075842235520591e3 -0.15925757362501867446e3 0.52177895084329669814e2 -0.99262530612592538979e1 0.72758579432539352184e0 0 0 0 0 0; 0 -0.95506826122820715686e0 0.10457091713389811971e2 -0.51988852784348173364e2 0.15485750987170610843e3 -0.30687251076237717714e3 0.42433443014978202692e3 -0.41695815935738584232e3 0.28993831774553865375e3 -0.13835205654344741174e3 0.41458767514717546128e2 -0.62843543720560487741e1 0.36488508570871329747e0 0 0 0 0; 0 0 -0.89358450029368883150e0 0.97153548381345233039e1 -0.47882181361696798414e2 0.14105174527639415815e3 -0.27545564288451795418e3 0.37326113384585831300e3 -0.35591415118897425470e3 0.23537219817841487113e3 -0.10078002209916911947e3 0.24486327517266714070e2 -0.31781915325611402199e1 0.21701391114437616116e0 0 0 0; 0 0 0 -0.82079151707601599340e0 0.88371484607293704266e1 -0.43018583681492318832e2 0.12468896705454297098e3 -0.23815343678333249804e3 0.31236211190696585630e3 -0.28229009022198974497e3 0.16669674231526758728e3 -0.59818136859120542509e2 0.13526028855965613027e2 -0.22391799149842653766e1 0.22922038452398770506e0 0 0; 0 0 0 0 -0.75926914003985711330e0 0.80869450361611184028e1 -0.38807566922312997740e2 0.11026873219399300975e3 -0.20445170004300330891e3 0.25511907237023259909e3 -0.20692908752510736686e3 0.10704747746062038309e3 -0.39160480492663455349e2 0.12374519230919224800e2 -0.33811178542850964393e1 0.59247568548574727441e0 0; 0 0 0 0 0 -0.73517682146919325041e0 0.77615133844626073094e1 -0.36756244064664377424e2 0.10222585002150178224e3 -0.18222790883588068825e3 0.20692908752510771174e3 -0.15554757761621697209e3 0.91166185986353555693e2 -0.52255896182014198048e2 0.31492793235845807035e2 -0.19103080454788523638e2 0.70504538217624886828e1; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; ];
310 DD_11 = sparse(DD_11);
311
312 %[1, -12, 66, -220, 495, -792, 924, -792, 495, -220, 66,-12, 1]
313 DD_12 = (diag(ones(m - 6, 1), 6) - 12 * diag(ones(m - 5, 1), 5) + 66 * diag(ones(m - 4, 1), 4) - 220 * diag(ones(m - 3, 1), 3) + 495 * diag(ones(m - 2, 1), 2) - 792 * diag(ones(m - 1, 1), 1) + 924 * diag(ones(m, 1), 0) - 792 * diag(ones(m - 1, 1), -1) + 495 * diag(ones(m - 2, 1), -2) - 220 * diag(ones(m - 3, 1), -3) + 66 * diag(ones(m - 4, 1), -4) - 12 * diag(ones(m - 5, 1), -5) + diag(ones(m - 6, 1), -6));
314 DD_12(1:12, 1:18) = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0.70504538217624945581e1 -0.19695556140274287325e2 0.34873911090130932535e2 -0.64630415412933476707e2 0.13032666647901711965e3 -0.26259505507683757401e3 0.41385817505021542347e3 -0.43734698120611365179e3 0.30667755006450534672e3 -0.14702497625865750970e3 0.46569080306775643856e2 -0.88221218576303190049e1 0.75926914003985774602e0 0 0 0 0 0; 0 0.56252054413792412864e0 -0.34278021535209417538e1 0.13874841106233553089e2 -0.50022717612825328606e2 0.15842900977125025845e3 -0.35473557861446977176e3 0.51023814474046519819e3 -0.49068408010320794140e3 0.33080619658197902926e3 -0.15523026768925199096e3 0.48521670216966710417e2 -0.91112296804782853596e1 0.77929289272158630803e0 0 0 0 0; 0 0 0.21428192906429919385e0 -0.22961219278627835050e1 0.15615594467315483385e2 -0.78810282483468544556e2 0.25004511347290138092e3 -0.48392586895198241994e3 0.62472422381393171260e3 -0.57156824827999799529e3 0.37406690116362891293e3 -0.17207433472596927533e3 0.53022890764376222560e2 -0.98494982049121919208e1 0.83534896297519895228e0 0 0 0; 0 0 0 0.20501361895561890151e0 -0.33471539032596363119e1 0.29069145008244604383e2 -0.13437336279889215930e3 0.35305829726762230670e3 -0.61013854489538443663e3 0.74652226769171662600e3 -0.66109354292284309003e3 0.42315523582918247446e3 -0.19152872544678719366e3 0.58292129028807139823e2 -0.10723014003524265978e2 0.90225552616201164306e0 0 0; 0 0 0 0 0.35327850260839005288e0 -0.67888982569607603271e1 0.49750521017661055353e2 -0.18446940872459654898e3 0.43490747661830798063e3 -0.71478541604123287255e3 0.84866886029956405384e3 -0.73649402582970522515e3 0.46457252961511832529e3 -0.20795541113739269345e3 0.62742550280338871828e2 -0.11460819134738485882e2 0.95876279102790935799e0 0; 0 0 0 0 0 0.72108566182178027310e0 -0.10828639703191913343e2 0.62613474101195603777e2 -0.21234343150002489927e3 0.48196613763353280887e3 -0.77522833310101589741e3 0.90763645940119576992e3 -0.77994731454339649814e3 0.48839612538974061864e3 -0.21738718877850875015e3 0.65293702141574540142e2 -0.11883142866454242025e2 0.99106616353107873319e0; ];
315 DD_12(m - 11:m, m - 17:m) = [0.99106616353107873319e0 -0.11883142866454242025e2 0.65293702141574540142e2 -0.21738718877850875015e3 0.48839612538974061864e3 -0.77994731454339649814e3 0.90763645940119576992e3 -0.77522833310101589741e3 0.48196613763353280887e3 -0.21234343150002489927e3 0.62613474101195603777e2 -0.10828639703191913343e2 0.72108566182178027310e0 0 0 0 0 0; 0 0.95876279102790935799e0 -0.11460819134738485882e2 0.62742550280338871828e2 -0.20795541113739269345e3 0.46457252961511832529e3 -0.73649402582970522515e3 0.84866886029956405384e3 -0.71478541604123287255e3 0.43490747661830798063e3 -0.18446940872459654898e3 0.49750521017661055353e2 -0.67888982569607603271e1 0.35327850260839005288e0 0 0 0 0; 0 0 0.90225552616201164306e0 -0.10723014003524265978e2 0.58292129028807139823e2 -0.19152872544678719366e3 0.42315523582918247446e3 -0.66109354292284309003e3 0.74652226769171662600e3 -0.61013854489538443663e3 0.35305829726762230670e3 -0.13437336279889215930e3 0.29069145008244604383e2 -0.33471539032596363119e1 0.20501361895561890151e0 0 0 0; 0 0 0 0.83534896297519895228e0 -0.98494982049121919208e1 0.53022890764376222560e2 -0.17207433472596927533e3 0.37406690116362891293e3 -0.57156824827999799529e3 0.62472422381393171260e3 -0.48392586895198241994e3 0.25004511347290138092e3 -0.78810282483468544556e2 0.15615594467315483385e2 -0.22961219278627835050e1 0.21428192906429919385e0 0 0; 0 0 0 0 0.77929289272158630803e0 -0.91112296804782853596e1 0.48521670216966710417e2 -0.15523026768925199096e3 0.33080619658197902926e3 -0.49068408010320794140e3 0.51023814474046519819e3 -0.35473557861446977176e3 0.15842900977125025845e3 -0.50022717612825328606e2 0.13874841106233553089e2 -0.34278021535209417538e1 0.56252054413792412864e0 0; 0 0 0 0 0 0.75926914003985774602e0 -0.88221218576303190049e1 0.46569080306775643856e2 -0.14702497625865750970e3 0.30667755006450534672e3 -0.43734698120611365179e3 0.41385817505021542347e3 -0.26259505507683757401e3 0.13032666647901711965e3 -0.64630415412933476707e2 0.34873911090130932535e2 -0.19695556140274287325e2 0.70504538217624945581e1; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; ];
316 DD_12 = sparse(DD_12);
317 %%%%%%%%%%%%%%%%%%%%%%%%%%%
318
319 %%%% Difference operators %%
320 D1 = H \ Q;
321
322 % Helper functions for constructing D2(c)
323 % TODO: Consider changing sparse(diag(...)) to spdiags(....)
324
325 % Minimal 13 point stencil width
326 function D2 = D2_fun_minimal(c)
327 % Here we add variable diffusion
328 C1 = sparse(diag(c));
329 C2 = 1/2 * diag(ones(m - 1, 1), -1) + 1/2 * diag(ones(m, 1), 0); C2(1, 2) = 1/2;
330 C3 = 5/18 * diag(ones(m - 1, 1), -1) + 5/18 * diag(ones(m - 1, 1), 1) + 4/9 * diag(ones(m, 1), 0); C3(1, 3) = 5/18; C3(m, m - 2) = 5/18;
331 C4 = 5/28 * diag(ones(m - 2, 1), -2) + 9/28 * diag(ones(m - 1, 1), -1) + 5/28 * diag(ones(m - 1, 1), 1) + 9/28 * diag(ones(m, 1), 0); C4(2, 4) = 5/28; C4(1, 3) = 5/28; C4(1, 4) = 9/28; C4(m, m - 2) = 5/28; C4(m, m - 3) = 5/28;
332 C5 = 1/7 * diag(ones(m - 2, 1), -2) + 1/7 * diag(ones(m - 2, 1), 2) + 8/35 * diag(ones(m - 1, 1), -1) + 8/35 * diag(ones(m - 1, 1), 1) + 9/35 * diag(ones(m, 1), 0); C5(2, 4) = 1/7; C5(1, 3) = 1/7; C5(1, 4) = 1/7; C5(1, 5) = 8/35; C5(2, 5) = 1/7; C5(m - 1, m - 4) = 1/7; C5(m, m - 4) = 8/35; C5(m, m - 3) = 1/7;
333 C6 = 1/6 * diag(ones(m - 3, 1), -3) + 1/6 * diag(ones(m - 2, 1), -2) + 1/6 * diag(ones(m - 2, 1), 2) + 1/6 * diag(ones(m - 1, 1), -1) + 1/6 * diag(ones(m - 1, 1), 1) + 1/6 * diag(ones(m, 1), 0);
334 C6(1, 4) = 1/6; C6(1, 5) = 1/6; C6(1, 6) = 1/6; C6(2, 5) = 1/6; C6(2, 6) = 1/6; C6(3, 6) = 1/6; C6(m - 1, m - 5) = 1/6; C6(m, m - 4) = 1/6; C6(m, m - 5) = 1/6;
335
336 C2 = sparse(diag(C2 * c));
337 C3 = sparse(diag(C3 * c));
338 C4 = sparse(diag(C4 * c));
339 C5 = sparse(diag(C5 * c));
340 C6 = sparse(diag(C6 * c));
341
342 % Remainder term added to wide second derivative operator
343 R = (1/30735936 / h) * transpose(DD_12) * C1 * DD_12 + (1/6403320 / h) * transpose(DD_11) * C2 * DD_11 + (1/1293600 / h) * transpose(DD_10) * C3 * DD_10 + (1/249480 / h) * transpose(DD_9) * C4 * DD_9 + (1/44352 / h) * transpose(DD_8) * C5 * DD_8 + (1/6468 / h) * transpose(DD_7) * C6 * DD_7;
344 D2 = D1 * C1 * D1 - H \ R;
345 end
346
347 % Few additional grid point in interior stencil cmp. to minimal
348 function D2 = D2_fun_nonminimal(c)
349 % Here we add variable diffusion
350 C1 = sparse(diag(c));
351 C2 = 1/2 * diag(ones(m - 1, 1), -1) + 1/2 * diag(ones(m, 1), 0); C2(1, 2) = 1/2;
352 C3 = 5/18 * diag(ones(m - 1, 1), -1) + 5/18 * diag(ones(m - 1, 1), 1) + 4/9 * diag(ones(m, 1), 0); C3(1, 3) = 5/18; C3(m, m - 2) = 5/18;
353 C4 = 5/28 * diag(ones(m - 2, 1), -2) + 9/28 * diag(ones(m - 1, 1), -1) + 5/28 * diag(ones(m - 1, 1), 1) + 9/28 * diag(ones(m, 1), 0); C4(2, 4) = 5/28; C4(1, 3) = 5/28; C4(1, 4) = 9/28; C4(m, m - 2) = 5/28; C4(m, m - 3) = 5/28;
354 C5 = 1/7 * diag(ones(m - 2, 1), -2) + 1/7 * diag(ones(m - 2, 1), 2) + 8/35 * diag(ones(m - 1, 1), -1) + 8/35 * diag(ones(m - 1, 1), 1) + 9/35 * diag(ones(m, 1), 0); C5(2, 4) = 1/7; C5(1, 3) = 1/7; C5(1, 4) = 1/7; C5(1, 5) = 8/35; C5(2, 5) = 1/7; C5(m - 1, m - 4) = 1/7; C5(m, m - 4) = 8/35; C5(m, m - 3) = 1/7;
355
356 C2 = sparse(diag(C2 * c));
357 C3 = sparse(diag(C3 * c));
358 C4 = sparse(diag(C4 * c));
359 C5 = sparse(diag(C5 * c));
360
361 % Remainder term added to wide second derivative operator
362 R = (1/30735936 / h) * transpose(DD_12) * C1 * DD_12 + (1/6403320 / h) * transpose(DD_11) * C2 * DD_11 + (1/1293600 / h) * transpose(DD_10) * C3 * DD_10 + (1/249480 / h) * transpose(DD_9) * C4 * DD_9 + (1/44352 / h) * transpose(DD_8) * C5 * DD_8;
363 D2 = D1 * C1 * D1 - H \ R;
364 end
365
366 % Wide stencil
367 function D2 = D2_wide(c)
368 % Here we add variable diffusion
369 C1 = sparse(diag(c));
370 D2 = D1 * C1 * D1;
371 end
372
373 switch options.stencil_width
374 case 'minimal'
375 D2 = @D2_fun_minimal;
376 case 'nonminimal'
377 D2 = @D2_fun_nonminimal;
378 case 'wide'
379 D2 = @D2_fun_wide;
380 otherwise
381 error('No option %s for stencil width', options.stencil_width)
382 end
383
384 %%%%%%%%%%%%%%%%%%%%%%%%%%%
385
386 %%%% Artificial dissipation operator %%%
387 switch options.AD
388 case 'upwind'
389 % This is the choice that yield 11th order Upwind
390 DI = H \ (transpose(DD_6) * DD_6) * (-1/5544);
391 case 'op'
392 % This choice will preserve the order of the underlying
393 % Non-dissipative D1 SBP operator
394 DI = H \ (transpose(DD_7) * DD_7) * (-1 / (5 * 5544));
395 % Notice that you can use any negative number instead of (-1/(5*5544))
396 otherwise
397 error("Artificial dissipation options '%s' not implemented.", option.AD)
398 end
399
400 %%%%%%%%%%%%%%%%%%%%%%%%%%%
401 end